pseudo.c

混沌分析的C语言程序的

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

Calculates the pseudo inverse for a given matrix 

Efficient version for prediction progranm cast_pred
M.Banbrook    6/10/94


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




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










/************** 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];
	    }
	}
    }
}

/************** multiply two matrices   *****************************/  
void multiply_matrixE(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 (k=1;k<=row2;k++)
	{
	  out_matrix[i][i]=out_matrix[i][i]+ matrix1[i][k] * matrix2[k][i];
	}
    }
}

/************** copy matrix ***************************************/
void copy_matrix(in_matrix,copy,num_row,num_col)
     int num_row,num_col;
     double **in_matrix,**copy;

{
  int i,j;

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






/***************** rearrange diagonal elements of E *******************
****************** matrix and reciprocate          *******************/
void  order_matrix(original,num)
     double **original;
     int num;

{
int i,j;
double tmp;

for (i=1; i<=num;i++)
  {
    for (j=i;j<=num;j++)
      {
	if (fabs(original[i][i])<fabs(original[j][j]))
	  {
	    tmp=original[j][j];
	    original[j][j]=original[i][i];
	    original[i][i]=tmp;
	  }
      }
  }

for (i=1; i<=num;i++)
  for (j=1;j<=num;j++)
    {
      if (i!=j)
	original[i][j]=0.0;
      else
	{
	  if (original[i][i]==0.0);
	  else
	    original[i][i]=1.0/original[i][i];
	}
    }
}

/****************** order the diagonal components   *******************/
void  order_E(original,num)
     double **original;
     int num;

{
int i,j;
double tmp;

for (i=1; i<=num;i++)
  {
    for (j=i;j<=num;j++)
      {
	if (fabs(original[i][i])<fabs(original[j][j]))
	  {
	    tmp=original[j][j];
	    original[j][j]=original[i][i];
	    original[i][i]=tmp;
	  }
      }
  }
}

/***************** 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];
}







/*****************  print out matrix to screen ********************/
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");
  }
}
/*****************************************************************/




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

{
double **stran,**E,**matrix_transpose,*values;
double **s,**corr,**tmp1,**tmp2,**tmp3,**tmp4;
int   i,j,nrot;

/************ initialise all matrices  ****************************/
matrix_transpose=dmatrix(1,num_col,1,num_row);
corr=dmatrix(1,num_col,1,num_col);
tmp1=dmatrix(1,num_col,1,num_col);
E=dmatrix(1,num_col,1,num_col);
values=dvector(1,num_row);

/*************  calculate C matrix ***********************************/
transpose(matrix_array,matrix_transpose,num_row,num_col);
multiply_matrix(tmp1,matrix_transpose,matrix_array,num_col,num_row,num_row,num_col);
jacobi(tmp1,num_col,values,corr,&nrot);
eigsrt(values,corr,num_col);

free_dmatrix(tmp1,1,num_col,1,num_col);

/*************** calculate S matrix  ********************************/
tmp2=dmatrix(1,num_row,1,num_row);

multiply_matrix(tmp2,matrix_array,matrix_transpose,num_row,num_col,num_col,num_row);
free_dmatrix(matrix_transpose,1,num_col,1,num_row);
s=dmatrix(1,num_row,1,num_row);
jacobi(tmp2,num_row,values,s,&nrot);
eigsrt(values,s,num_row);


free_dmatrix(tmp2,1,num_row,1,num_row);



/************** calculate E matrix  **********************************/
stran=dmatrix(1,num_col,1,num_row);
tmp3=dmatrix(1,num_row,1,num_col);

transpose(s,stran,num_row,num_col);
multiply_matrix(tmp3,matrix_array,corr,num_row,num_col,num_col,num_col);
multiply_matrixE(E,stran,tmp3,num_col,num_row,num_row,num_col);
order_matrix(E,num_col);


free_dmatrix(tmp3,1,num_row,1,num_col);
free_dmatrix(s,1,num_row,1,num_row);

/*********** calculate B+ matrix  ***********************************/
tmp4=dmatrix(1,num_col,1,num_row);

multiply_matrix(tmp4,E,stran,num_col,num_col,num_col,num_row);
multiply_matrix(matrix_inv,corr,tmp4,num_col,num_col,num_col,num_row);

free_dmatrix(E,1,num_col,1,num_col);
free_dmatrix(corr,1,num_col,1,num_col);
free_dmatrix(stran,1,num_col,1,num_row);
free_dmatrix(tmp4,1,num_col,1,num_row);
free_dvector(values,1,num_row);
}