recipes_pseudo.c~

混沌分析的C语言程序的

/***********************************************************************

Calculates the pseudo inverse for a given matrix using Recipes svd
routine

M.Banbrook 12/5/95



**********************************************************************/




#include <stdio.h>
#include <math.h>
#include "/home/mb/Recipes/nrutil.c"
#include "/home/mb/Recipes/svdcmp.c"
#include "/home/mb/Recipes/jacobi.c"
#include "/home/mb/Recipes/eigsrt.c"


/***************** print a matrix  ******************************/
void print_matrix(array,row,col)
     double  **array;
     int    row,col;

{
int i,j;

printf("\n");
for (i=1;i<=row;i++)
  {
    for (j=1;j<=col;j++)
      printf(" %lf ",array[i][j]);
    printf("\n");
  }
}

/***************** transpose a matrix  ******************************/
void transpose(double **in_matrix,double **tran,int num_row,int num_col)

{
  int i,j;

  
  for (i=1;i<=num_col;i++)
    for (j=1;j<=num_row;j++)
	tran[i][j]=in_matrix[j][i];
}



/************** multiply two matrices   *****************************/  
void multiply_matrix(out_matrix,matrix1,matrix2,row1,col1,row2,col2)
     int row1,col1,row2,col2;
     double **out_matrix,**matrix1,**matrix2;

{
  int i,j,k;
 
  for (i=1;i<=row1;i++)
    for (j=1;j<=col2;j++)
      out_matrix[i][j]=0.0;

  for (i=1;i<=row1;i++)
    {
      for (j=1;j<=col2;j++)
	{
	  for (k=1;k<=row2;k++)
	    {
	      out_matrix[i][j]=out_matrix[i][j]+ matrix1[i][k] * matrix2[k][j];
	    }
	}
    }
}



void pseudo(matrix_inv,matrix_array,num_row,num_col)
     double **matrix_array,**matrix_inv;
     int   num_row,num_col;

{
  double **s;
  double **corr,*E;
  int   i,j,k;

  s=dmatrix(1,num_row,1,num_col);
  E=dvector(1,num_col);
  corr=dmatrix(1,num_col,1,num_col);

  for (i=1;i<=num_col;i++)
    for (j=1;j<=num_row;j++)
      matrix_inv[i][j]=0.0;
  
   
  for(i=1;i<=num_row;i++)
    for(j=1;j<=num_col;j++)
      s[i][j]=matrix_array[i][j];
 
  svdcmp(s,num_row,num_col,E,corr);

  
  for (i=1;i<=num_col;i++)
    for (j=1;j<=num_row;j++)
      for (k=1;k<=num_col;k++)
	matrix_inv[i][j] += s[j][k] * corr[i][k] * (1/E[k]);

}