/***********************************************************************Calculates the pseudo inverse for a given matrix using Recipes svdroutineM.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]);}