1 /*** These are from numerical_recipes in C. */
2 /*** PROBLEM: Thier counts went from 1 to N. **/
3 /*** I Changed that. ***/
8 #define MAX_MATRIX_SIZE MAX_NUM_SCORE_DIM
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;}
14 double a[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],b[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE];
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();
25 for (j=0;j<n;j++) ipiv[j]=0;
32 if (fabs(a[j][k]) >= big) {
37 } else if (ipiv[k] > 1) nrerror("GAUSSJ: Singular Matrix-1");
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])
46 if (a[icol][icol] == 0.0) nrerror("GAUSSJ: Singular Matrix-2");
47 pivinv=1.0/a[icol][icol];
49 for (l=0;l<n;l++) a[icol][l] *= pivinv;
50 for (l=0;l<m;l++) b[icol][l] *= pivinv;
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;
59 for (l=n-1;l>=0;l--) {
60 if (indxr[l] != indxc[l])
62 SWAP(a[k][indxr[l]],a[k][indxc[l]]);
64 free_ivector(ipiv,1,n);
65 free_ivector(indxr,1,n);
66 free_ivector(indxc,1,n);
73 /***************************************************************/
78 void ludcmp(a,n,indx,d)
80 double a[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],*d;
83 double big,dum,sum,temp;
85 void nrerror(),free_vector();
87 /* fprintf(stderr,"\n Declaring space for vv"); */
90 /* fprintf(stderr,"\nSuccessfully declared vv"); */
96 if ((temp=fabs(a[i][j])) > big) big=temp;
98 if (big == 0.0) nrerror("Singular matrix in routine LUDCMP");
105 for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
112 sum -= a[i][k]*a[k][j];
114 if ( (dum=vv[i]*fabs(sum)) >= big) {
129 if (a[j][j] == 0.0) a[j][j]=TINY;
132 for (i=j+1;i<n;i++) a[i][j] *= dum;
136 /* free_vector(vv,1,n); */
137 free_vector(vv,0,n-1);
146 /***********************************/
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];
160 for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
164 for (i=n-1;i>=0;i--) {
166 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
173 /*****************************************************/
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
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],
187 double temp_matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE];
188 double col[MAX_MATRIX_SIZE];
189 int indx[MAX_MATRIX_SIZE];
193 for (j=0; j<n; j++) /** So matrix doesn't get destroyed. **/
194 temp_matrix[i][j] = matrix[i][j];
196 /** fprintf(stderr,"\n Matrix to invert + find det of is:\n");
197 print_matrix(matrix,n,n,stderr);
200 /* fprintf(stderr,"\n Into ludcmp"); */
202 ludcmp(temp_matrix,n,indx,det);
206 *det *= temp_matrix[j][j];
208 /* fprintf(stderr,"\n did det=%lf",*det); */
210 /*** Now convert current form of "inverse" into the real matrix inverse. */
211 /** only if want to compute inverse though. ****/
213 for (j=0; j<n; j++) {
218 lubksb(temp_matrix,n,indx,col);
220 inverse[i][j] = col[i];
227 /*************************************************************************/
230 double determ(double matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE], int n)
234 matrix_inverse_and_det (matrix, NULL, &det,n);
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])
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];
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])
273 double x_trans[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE];
275 double product1[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE];
276 int all_x_non_infinity =1;
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;
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);
291 matrix_mult(num_x_column, num_M_row_column, num_M_row_column, x_trans, M,
293 matrix_mult(num_x_column, num_M_row_column, num_x_column, product1, x,
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);
309 if (all_x_non_infinity) return 1;
319 print_matrix(double M[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE], int num_rows,
320 int num_columns, FILE *f_out)
325 for (i=0; i<num_rows; i++){
327 for (j=0;j<num_columns;j++)
328 fprintf(f_out," %lf",M[i][j]);
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)
341 for (i=0; i<num_rows; i++){
343 for (j=0;j<num_columns;j++)
344 fprintf(f_out," %lf",M[i][j]);