sources/multicoil/matrix_operations.c

   1 /*** These are from numerical_recipes  in  C. */
   2 /*** PROBLEM: Thier counts went from 1 to N. **/
   3 /*** I Changed that. ***/
   4 #include <stdio.h>
   5 #include <math.h>
   6 #include "scconst.h"
   7
   8 #define MAX_MATRIX_SIZE MAX_NUM_SCORE_DIM
   9
  10 /*** gaussj replaces a by its inverse?? If b starts as identity matrix?? ***/
  11 #define SWAP(a,b) {double temp=(a);(a)=(b);(b)=temp;}
  12
  13 void gaussj(a,n,b,m)
  14 double a[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],b[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE];
  15 int n,m;
  16 {
  17         int *indxc,*indxr,*ipiv;
  18         int i,icol,irow,j,k,l,ll,*ivector();
  19         double big,dum,pivinv;
  20         void nrerror(),free_ivector();
  21
  22         indxc=ivector(1,n);
  23         indxr=ivector(1,n);
  24         ipiv=ivector(1,n);
  25         for (j=0;j<n;j++) ipiv[j]=0;
  26         for (i=0;i<n;i++) {
  27                 big=0.0;
  28                 for (j=0;j<n;j++)
  29                         if (ipiv[j] != 1)
  30                                 for (k=0;k<n;k++) {
  31                                         if (ipiv[k] == 0) {
  32                                                 if (fabs(a[j][k]) >= big) {
  33                                                         big=fabs(a[j][k]);
  34                                                         irow=j;
  35                                                         icol=k;
  36                                                 }
  37                                         } else if (ipiv[k] > 1) nrerror("GAUSSJ: Singular Matrix-1");
  38                                 }
  39                 ++(ipiv[icol]);
  40                 if (irow != icol) {
  41                         for (l=0;l<n;l++) SWAP(a[irow][l],a[icol][l])
  42                         for (l=0;l<m;l++) SWAP(b[irow][l],b[icol][l])
  43                 }
  44                 indxr[i]=irow;
  45                 indxc[i]=icol;
  46                 if (a[icol][icol] == 0.0) nrerror("GAUSSJ: Singular Matrix-2");
  47                 pivinv=1.0/a[icol][icol];
  48                 a[icol][icol]=1.0;
  49                 for (l=0;l<n;l++) a[icol][l] *= pivinv;
  50                 for (l=0;l<m;l++) b[icol][l] *= pivinv;
  51                 for (ll=0;ll<n;ll++)
  52                         if (ll != icol) {
  53                                 dum=a[ll][icol];
  54                                 a[ll][icol]=0.0;
  55                                 for (l=0;l<n;l++) a[ll][l] -= a[icol][l]*dum;
  56                                 for (l=0;l<m;l++) b[ll][l] -= b[icol][l]*dum;
  57                         }
  58         }
  59         for (l=n-1;l>=0;l--) {
  60                 if (indxr[l] != indxc[l])
  61                         for (k=0;k<n;k++)
  62                                 SWAP(a[k][indxr[l]],a[k][indxc[l]]);
  63         }
  64         free_ivector(ipiv,1,n);
  65         free_ivector(indxr,1,n);
  66         free_ivector(indxc,1,n);
  67 }
  68
  69 #undef SWAP
  70
  71
  72
  73 /***************************************************************/
  74
  75
  76 #define TINY 1.0e-20;
  77
  78 void ludcmp(a,n,indx,d)
  79 int n,*indx;
  80 double a[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],*d;
  81 {
  82         int i,imax,j,k;
  83         double big,dum,sum,temp;
  84         double *vv,*vector();
  85         void nrerror(),free_vector();
  86
  87 /*      fprintf(stderr,"\n Declaring space for vv"); */
  88 /*      vv=vector(1,n); */
  89         vv=vector(0,n-1);
  90 /*      fprintf(stderr,"\nSuccessfully declared vv"); */
  91
  92         *d=1.0;
  93         for (i=0;i<n;i++) {
  94                 big=0.0;
  95                 for (j=0;j<n;j++) {
  96                         if ((temp=fabs(a[i][j])) > big) big=temp;
  97                       }
  98                 if (big == 0.0) nrerror("Singular matrix in routine LUDCMP");
  99                 vv[i]=1.0/big;
 100         }
 101
 102         for (j=0;j<n;j++) {
 103                 for (i=0;i<j;i++) {
 104                         sum=a[i][j];
 105                         for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
 106                         a[i][j]=sum;
 107                 }
 108                 big=0.0;
 109                 for (i=j;i<n;i++) {
 110                         sum=a[i][j];
 111                         for (k=0;k<j;k++)
 112                                 sum -= a[i][k]*a[k][j];
 113                         a[i][j]=sum;
 114                         if ( (dum=vv[i]*fabs(sum)) >= big) {
 115                                 big=dum;
 116                                 imax=i;
 117                         }
 118                 }
 119                 if (j != imax) {
 120                         for (k=0;k<n;k++) {
 121                                 dum=a[imax][k];
 122                                 a[imax][k]=a[j][k];
 123                                 a[j][k]=dum;
 124                         }
 125                         *d = -(*d);
 126                         vv[imax]=vv[j];
 127                 }
 128                 indx[j]=imax;
 129                 if (a[j][j] == 0.0) a[j][j]=TINY;
 130                 if (j != n-1) {
 131                         dum=1.0/(a[j][j]);
 132                         for (i=j+1;i<n;i++) a[i][j] *= dum;
 133                 }
 134         }
 135
 136 /*      free_vector(vv,1,n);  */
 137         free_vector(vv,0,n-1);
 138
 139 }
 140
 141 #undef TINY
 142
 143
 144
 145
 146 /***********************************/
 147
 148 void lubksb(a,n,indx,b)
 149 double a[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],b[MAX_MATRIX_SIZE];
 150 int n,indx[MAX_MATRIX_SIZE];
 151 {
 152         int i,ii=-1,ip,j;
 153         double sum;
 154
 155         for (i=0;i<n;i++) {
 156                 ip=indx[i];
 157                 sum=b[ip];
 158                 b[ip]=b[i];
 159                 if (ii != -1)
 160                         for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
 161                 else if (sum) ii=i;
 162                 b[i]=sum;
 163         }
 164         for (i=n-1;i>=0;i--) {
 165                 sum=b[i];
 166                 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
 167                 b[i]=sum/a[i][i];
 168         }
 169 }
 170
 171
 172
 173 /*****************************************************/
 174
 175
 176 /***
 177   matrix and inverse should be n X n arrays.
 178   indx should be an array of n
 179   col should be an array of n
 180 ****/
 181 /*** If just want det then set inverse to NULL ***/
 182 void matrix_inverse_and_det (double matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 183                              double inverse[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 184                              double *det,int n)
 185 {
 186   int i,j;
 187   double  temp_matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE];
 188   double col[MAX_MATRIX_SIZE];
 189   int indx[MAX_MATRIX_SIZE];
 190
 191
 192   for (i=0; i<n; i++)
 193     for (j=0; j<n; j++)      /** So matrix doesn't get destroyed. **/
 194       temp_matrix[i][j] = matrix[i][j];
 195
 196 /**  fprintf(stderr,"\n Matrix to invert + find det of is:\n");
 197   print_matrix(matrix,n,n,stderr);
 198 ***/
 199
 200 /*  fprintf(stderr,"\n Into ludcmp"); */
 201
 202   ludcmp(temp_matrix,n,indx,det);
 203
 204
 205   for (j=0; j<n; j++)
 206     *det *= temp_matrix[j][j];
 207
 208 /*  fprintf(stderr,"\n did det=%lf",*det); */
 209
 210 /*** Now convert current form of "inverse" into the real matrix inverse.  */
 211 /** only if want to compute inverse though. ****/
 212   if (inverse)
 213     for (j=0; j<n; j++) {
 214       for (i=0; i<n; i++)
 215         col[i] = 0.0;
 216       col[j]=1.0;
 217
 218       lubksb(temp_matrix,n,indx,col);
 219       for (i=0; i<n; i++)
 220         inverse[i][j] = col[i];
 221     }
 222
 223 }
 224
 225
 226
 227 /*************************************************************************/
 228
 229
 230 double determ(double matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE], int n)
 231 {
 232   double det;
 233
 234   matrix_inverse_and_det (matrix, NULL,   &det,n);
 235
 236   return (det);
 237 }
 238
 239
 240
 241 void  matrix_mult(int num_rows1, int num_columns1_rows2, int num_columns2,
 242                   double matrix1[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 243                   double matrix2[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 244                   double matrix_product[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE])
 245 {
 246
 247   int i, j, k;
 248
 249   for (i=0; i<num_rows1; i++)
 250     for(k=0; k<num_columns2; k++) {
 251       matrix_product[i][k]=0;
 252       for (j=0; j<num_columns1_rows2; j++)
 253         matrix_product[i][k] += matrix1[i][j] * matrix2[j][k];
 254     }
 255
 256 }
 257
 258
 259
 260
 261
 262
 263 /** Return 0 if x had some -HUGE_VAL values in it, since that means **/
 264 /** the normalized matrix is grabage. ***?
 265 /*** compute x(transpose)  * M * x ******/
 266 int normalize (int num_M_row_column,
 267                 int num_x_column, int num_x_row,
 268                 double out_matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 269                  double M[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 270                 double x[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE])
 271
 272 {
 273   double x_trans[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE];
 274   int i,k;
 275   double product1[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE];
 276   int all_x_non_infinity =1;
 277
 278   for(i=0; i<num_x_column; i++)
 279     for (k=0; k < num_x_row; k++) {
 280       x_trans[i][k] = x[k][i];
 281       if (x[k][i] <= -HUGE_VAL)
 282         all_x_non_infinity=0;
 283     }
 284
 285
 286   if (num_x_row != num_M_row_column) {   /** num_x_row = num_xtrans_column. **/
 287     fprintf(stderr, "\n Error:  Can't multiply transpose by matrix because number_x_row = %d, number_M_row_column=%d.", num_x_row, num_M_row_column);
 288     exit(-1);
 289   }
 290
 291   matrix_mult(num_x_column, num_M_row_column, num_M_row_column, x_trans, M,
 292               product1);
 293   matrix_mult(num_x_column, num_M_row_column, num_x_column, product1, x,
 294               out_matrix);
 295
 296 /***********
 297   if (all_x_non_infinity) {
 298     fprintf(stderr,"\n x_trans=");
 299     print_matrix(x_trans, num_x_column, num_M_row_column,stderr);
 300     fprintf(stderr,"\n M=");
 301     print_matrix(M, num_M_row_column, num_M_row_column,stderr);
 302     fprintf(stderr,"\n x_trans*M=");
 303     print_matrix(product1, num_x_column, num_M_row_column,stderr);
 304     fprintf(stderr,"\n x_trans*M*x=");
 305     print_matrix(out_matrix, num_x_column,num_x_column,stderr);
 306   }
 307 *************/
 308
 309   if (all_x_non_infinity) return 1;
 310   else return 0;
 311
 312 }
 313
 314
 315
 316
 317
 318
 319 print_matrix(double M[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE], int num_rows,
 320                   int num_columns, FILE *f_out)
 321 {
 322   int i,j;
 323
 324
 325   for (i=0; i<num_rows; i++){
 326     fprintf(f_out,"\n");
 327     for (j=0;j<num_columns;j++)
 328       fprintf(f_out," %lf",M[i][j]);
 329   }
 330   fprintf(f_out,"\n");
 331 }
 332
 333
 334 /*** For printing a matrix of size used for computing ALL dimensions. ***/
 335 print_matrix2(double M[NUM_DIM_IN_ORIG_MATRIX][NUM_DIM_IN_ORIG_MATRIX],
 336               int num_rows, int num_columns, FILE *f_out)
 337 {
 338   int i,j;
 339
 340
 341   for (i=0; i<num_rows; i++){
 342     fprintf(f_out,"\n");
 343     for (j=0;j<num_columns;j++)
 344       fprintf(f_out," %lf",M[i][j]);
 345   }
 346   fprintf(f_out,"\n");
 347 }