Doxygen Source Code Documentation
matrix_f.c File Reference
#include "mri_image.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "matrix_f.h"
#include <string.h>
Go to the source code of this file.
Define Documentation
|
|
|
|
|
|
|
Definition at line 990 of file matrix_f.c. |
|
|
Function Documentation
|
Convert simple array to matrix structure. Definition at line 378 of file matrix_f.c. References cols, matrix::elts, i, matrix_create(), and rows. Referenced by calc_reduced_model(), and read_glt_matrix().
|
|
Convert simple array f into vector v. Definition at line 898 of file matrix_f.c. References vector::elts, i, v, and vector_create_noinit(). Referenced by calc_reduced_model(), form_clusters(), and print_cluster().
|
|
Convert column c of matrix m into vector v. Definition at line 915 of file matrix_f.c. References c, matrix::elts, vector::elts, i, matrix::rows, v, and vector_create_noinit(). Referenced by init_delay(), and init_regression_analysis().
|
|
Add matrix a to matrix b. Result is matrix c. Definition at line 487 of file matrix_f.c. References a, c, matrix::cols, cols, matrix::elts, flops, i, matrix_create(), matrix_error(), matrix::rows, and rows. Referenced by ComputeDeltaTau(), and matrix_sqrt().
00488 { 00489 register int rows, cols; 00490 register int i, j; 00491 00492 if ((a.rows != b.rows) || (a.cols != b.cols)) 00493 matrix_error ("Incompatible dimensions for matrix addition"); 00494 00495 rows = a.rows; 00496 cols = a.cols; 00497 00498 matrix_create (rows, cols, c); 00499 00500 for (i = 0; i < rows; i++) 00501 for (j = 0; j < cols; j++) 00502 c->elts[i][j] = a.elts[i][j] + b.elts[i][j]; 00503 } |
|
Search a matrix for nearly identical column pairs, where "nearly identical" means they are correlated closer than 1-eps. Return is a pointer to an int array of the form [ i1 j1 i2 j2 ... -1 -1 ] where columns (i1,j1) are nearly the same, (i2,j2) also, etc. In addition:
Definition at line 1225 of file matrix_f.c. References a, matrix::cols, cols, matrix::elts, i, realloc, matrix::rows, and rows. Referenced by calculate_results().
01226 { 01227 int i,j,k , rows=a.rows , cols=a.cols ; 01228 int *iar=NULL , nar=0 ; 01229 double sumi,sumj,sumd ; 01230 01231 if( eps <= 0.0f ) eps = 1.e-5 ; 01232 01233 for( i=0 ; i < cols ; i++ ){ 01234 sumi = 0.0 ; 01235 for( k=0 ; k < rows ; k++ ) sumi += a.elts[k][i] * a.elts[k][i] ; 01236 if( sumi <= 0.0 ){ 01237 iar = (int *)realloc( (void *)iar , sizeof(int)*2*(nar+1) ) ; 01238 iar[2*nar] = i ; iar[2*nar+1] = -1 ; nar++ ; 01239 continue ; /* skip to next column i */ 01240 } 01241 for( j=i+1 ; j < cols ; j++ ){ 01242 sumj = sumd = 0.0 ; 01243 for( k=0 ; k < rows ; k++ ){ 01244 sumj += a.elts[k][j] * a.elts[k][j] ; 01245 sumd += a.elts[k][j] * a.elts[k][i] ; 01246 } 01247 if( sumj > 0.0 ){ 01248 sumd = fabs(sumd) / sqrt(sumi*sumj) ; 01249 if( sumd >= 1.0-eps ){ 01250 iar = (int *)realloc( (void *)iar , sizeof(int)*2*(nar+1) ) ; 01251 iar[2*nar] = i ; iar[2*nar+1] = j ; nar++ ; 01252 } 01253 } 01254 } 01255 } 01256 01257 if( iar != NULL ){ 01258 iar = (int *)realloc( (void *)iar , sizeof(int)*2*(nar+1) ) ; 01259 iar[2*nar] = iar[2*nar+1] = -1 ; 01260 } 01261 01262 return iar ; 01263 } |
|
Create matrix data structure by allocating memory and initializing values. Definition at line 161 of file matrix_f.c. References calloc, matrix::cols, cols, matrix::elts, i, malloc, matrix::mat, matrix_destroy(), matrix_error(), matrix::rows, and rows. Referenced by analyze_results(), array_to_matrix(), calc_covariance(), ComputeD0(), ComputeNewD(), cubic_spline(), Form_R_Matrix(), get_inputs(), init_indep_var_matrix(), InitGlobals(), matrix_add(), matrix_enter(), matrix_equate(), matrix_extract(), matrix_extract_rows(), matrix_file_read(), matrix_identity(), matrix_multiply(), matrix_psinv(), matrix_scale(), matrix_subtract(), matrix_transpose(), niml_to_matrix(), and poly_field().
00162 { 00163 register int i, j; 00164 00165 matrix_destroy (m); 00166 00167 if ((rows < 0) || (cols < 0)) 00168 matrix_error ("Illegal dimensions for new matrix"); 00169 00170 m->rows = rows; 00171 m->cols = cols; 00172 if ((rows < 1) || (cols < 1)) return; 00173 00174 m->elts = (float **) malloc (sizeof(float *) * rows); 00175 if (m->elts == NULL) 00176 matrix_error ("Memory allocation error"); 00177 00178 #ifdef DONT_USE_MATRIX_MAT 00179 for (i = 0; i < rows; i++){ 00180 m->elts[i] = (float *) calloc (sizeof(float) , cols); 00181 if (m->elts[i] == NULL) matrix_error ("Memory allocation error"); 00182 } 00183 #else 00184 m->mat = (float *) calloc( sizeof(float) , rows*cols ) ; 00185 if( m->mat == NULL ) 00186 matrix_error ("Memory allocation error"); 00187 for (i = 0; i < rows; i++) 00188 m->elts[i] = m->mat + (i*cols) ; /* 04 Mar 2005: offsets into mat */ 00189 #endif 00190 } |
|
|
Manual entry of matrix data. Definition at line 278 of file matrix_f.c. References cols, matrix::elts, i, matrix_create(), and rows.
00279 { 00280 int rows, cols; 00281 int i, j; 00282 float fval; 00283 00284 printf ("Enter number of rows: "); fflush(stdout); 00285 scanf ("%d", &rows); 00286 printf ("Enter number of cols: "); fflush(stdout); 00287 scanf ("%d", &cols); 00288 00289 matrix_create (rows, cols, m); 00290 00291 for (i = 0; i < rows; i++) 00292 for (j = 0; j < cols; j++) 00293 { 00294 printf ("elts[%d][%d] = ", i, j); fflush(stdout); 00295 scanf ("%f", &fval); 00296 m->elts[i][j] = fval; 00297 } 00298 } |
|
Make a copy of the first matrix, return copy as the second matrix. Definition at line 395 of file matrix_f.c. References a, matrix::cols, cols, matrix::elts, i, matrix_create(), matrix::rows, and rows. Referenced by init_indep_var_matrix(), matrix_inverse(), matrix_inverse_dsc(), and matrix_sqrt().
00396 { 00397 register int i, j; 00398 register int rows, cols; 00399 00400 rows = a.rows; 00401 cols = a.cols; 00402 00403 matrix_create (rows, cols, b); 00404 00405 for (i = 0; i < rows; i++){ 00406 #if 0 00407 for (j = 0; j < cols; j++) 00408 b->elts[i][j] = a.elts[i][j]; 00409 #else 00410 if( cols > 0 ) 00411 memcpy( b->elts[i] , a.elts[i] , sizeof(float )*cols ) ; /* RWCox */ 00412 #endif 00413 } 00414 } |
|
Routine to print and error message and stop. Definition at line 112 of file matrix_f.c. Referenced by matrix_add(), matrix_create(), matrix_file_read(), matrix_file_write(), matrix_identity(), matrix_inverse(), matrix_inverse_dsc(), matrix_multiply(), matrix_sqrt(), matrix_subtract(), vector_add(), vector_create(), vector_create_noinit(), vector_dot(), vector_multiply(), vector_multiply_subtract(), and vector_subtract().
00113 {
00114 printf ("Matrix error: %s \n", message);
00115 exit (1);
00116 }
|
|
Extract p columns (specified by list) from matrix a. Result is matrix b. Definition at line 422 of file matrix_f.c. References a, cols, matrix::elts, i, matrix_create(), p, matrix::rows, and rows. Referenced by calc_matrices().
|
|
Extract p rows (specified by list) from matrix a. Result is matrix b. Definition at line 443 of file matrix_f.c. References a, matrix::cols, cols, matrix::elts, i, matrix_create(), p, and rows. Referenced by init_indep_var_matrix().
|
|
Read contents of matrix m from specified file. If unable to read matrix from file, or matrix has wrong dimensions: If error_exit flag is set, then print error message and exit. Otherwise, return null matrix. Definition at line 309 of file matrix_f.c. References cols, matrix::elts, error_exit(), far, i, matrix_create(), matrix_destroy(), matrix_error(), MRI_FLOAT_PTR, mri_free(), mri_read_1D(), MRI_IMAGE::nx, MRI_IMAGE::ny, and rows. Referenced by DC_main(), do_xrestore_stuff(), markov_array(), read_glt_matrix(), and read_input_data().
00311 { 00312 int i, j; 00313 00314 MRI_IMAGE *im, *flim; /* pointers to image structures 00315 -- used to read ASCII file */ 00316 float * far; /* pointer to MRI_IMAGE floating point data */ 00317 char message [80]; /* error message */ 00318 00319 00320 /*----- First, check for empty file name -----*/ 00321 if (filename == NULL) matrix_error ("Missing matrix file name"); 00322 00323 00324 /*----- Read the matrix file -----*/ 00325 flim = mri_read_1D (filename); 00326 if (flim == NULL) 00327 if (error_exit) 00328 { 00329 sprintf (message, "Unable to read matrix from file: %s", filename); 00330 matrix_error (message); 00331 } 00332 else 00333 { 00334 matrix_destroy (m); 00335 return; 00336 } 00337 00338 00339 /*----- Set pointer to data -----*/ 00340 far = MRI_FLOAT_PTR(flim); 00341 00342 00343 /*----- Test for correct dimensions -----*/ 00344 if ( (rows != flim->nx) || (cols != flim->ny) ) 00345 if (error_exit) 00346 { 00347 sprintf (message, 00348 "In matrix file: %s Expected: %d x %d Actual: %d x %d", 00349 filename, rows, cols, flim->nx, flim->ny); 00350 matrix_error (message); 00351 } 00352 else 00353 { 00354 matrix_destroy (m); 00355 return; 00356 } 00357 00358 00359 matrix_create (rows, cols, m); 00360 00361 00362 /*----- Copy data from image structure to matrix structure -----*/ 00363 for (i = 0; i < rows; i++) 00364 for (j = 0; j < cols; j++) 00365 m->elts[i][j] = far[i + j*rows]; 00366 00367 00368 mri_free (flim); flim = NULL; 00369 00370 } |
|
Print contents of matrix m to specified file. Definition at line 246 of file matrix_f.c. References matrix::cols, cols, matrix::elts, i, matrix_error(), matrix::rows, and rows.
00247 { 00248 int i, j; 00249 int rows, cols; 00250 FILE * outfile = NULL; 00251 00252 00253 /*----- First, check for empty file name -----*/ 00254 if (filename == NULL) matrix_error ("Missing matrix file name"); 00255 00256 00257 outfile = fopen (filename, "w"); 00258 00259 rows = m.rows; 00260 cols = m.cols; 00261 00262 for (i = 0; i < rows; i++) 00263 { 00264 for (j = 0; j < cols; j++) 00265 fprintf (outfile, " %g", m.elts[i][j]); 00266 fprintf (outfile, " \n"); 00267 } 00268 fprintf (outfile, " \n"); 00269 00270 fclose (outfile); 00271 } |
|
Create n x n identity matrix. Definition at line 464 of file matrix_f.c. References matrix::elts, i, matrix_create(), and matrix_error(). Referenced by calc_glt_matrix(), form_clusters(), matrix_inverse(), and matrix_sqrt().
00465 { 00466 register int i, j; 00467 00468 if (n < 0) 00469 matrix_error ("Illegal dimensions for identity matrix"); 00470 00471 matrix_create (n, n, m); 00472 00473 for (i = 0; i < n; i++) 00474 for (j = 0; j < n; j++) 00475 if (i == j) 00476 m->elts[i][j] = 1.0; 00477 else 00478 m->elts[i][j] = 0.0; 00479 } |
|
Initialize matrix data structure. Definition at line 124 of file matrix_f.c. References matrix::cols, matrix::elts, matrix::mat, and matrix::rows. Referenced by analyze_results(), calc_covariance(), calc_glt_matrix(), calc_linear_regression(), calc_matrices(), calc_reduced_model(), calculate_results(), ComputeD0(), ComputeDeltaTau(), ComputeHpHm(), ComputeNewD(), cubic_spline(), do_xrestore_stuff(), form_clusters(), Form_R_Matrix(), init_delay(), init_indep_var_matrix(), init_regression_analysis(), InitGlobals(), initialize_options(), initialize_program(), markov_array(), matrix_destroy(), matrix_inverse(), matrix_inverse_dsc(), matrix_sqrt(), poly_field(), and read_input_data().
|
|
Use Gaussian elimination to calculate inverse of matrix a. Result is matrix ainv. Definition at line 607 of file matrix_f.c. References a, matrix::cols, matrix::elts, flops, i, matrix_destroy(), matrix_equate(), matrix_error(), matrix_identity(), matrix_initialize(), matrix_sprint(), p, and matrix::rows. Referenced by analyze_results(), calc_covariance(), calc_linear_regression(), cubic_spline(), matrix_inverse_dsc(), matrix_sqrt(), and poly_field().
00608 { 00609 const float epsilon = 1.0e-10; 00610 matrix tmp; 00611 register int i, j, ii, n; 00612 register float fval; 00613 register float fmax; 00614 register float * p; 00615 00616 matrix_initialize (&tmp); 00617 00618 00619 if (a.rows != a.cols) 00620 matrix_error ("Illegal dimensions for matrix inversion"); 00621 00622 00623 n = a.rows; 00624 matrix_identity (n, ainv); 00625 matrix_equate (a, &tmp); 00626 00627 for (i = 0; i < n; i++) 00628 { 00629 fmax = fabs(tmp.elts[i][i]); 00630 for (j = i+1; j < n; j++) 00631 if (fabs(tmp.elts[j][i]) > fmax) 00632 { 00633 fmax = fabs(tmp.elts[j][i]); 00634 p = tmp.elts[i]; 00635 tmp.elts[i] = tmp.elts[j]; 00636 tmp.elts[j] = p; 00637 p = ainv->elts[i]; 00638 ainv->elts[i] = ainv->elts[j]; 00639 ainv->elts[j] = p; 00640 } 00641 00642 if (fmax < epsilon) 00643 { 00644 matrix_destroy (&tmp); 00645 return (0); 00646 } 00647 00648 00649 fval = 1.0 / tmp.elts[i][i]; /* RWCox: change division by this to */ 00650 for (j = 0; j < n; j++) /* multiplication by 1.0/this */ 00651 { 00652 tmp.elts[i][j] *= fval; 00653 ainv->elts[i][j] *= fval; 00654 } 00655 for (ii = 0; ii < n; ii++) 00656 if (ii != i) 00657 { 00658 fval = tmp.elts[ii][i]; 00659 for (j = 0; j < n; j++) 00660 { 00661 tmp.elts[ii][j] -= fval*tmp.elts[i][j]; 00662 ainv->elts[ii][j] -= fval*ainv->elts[i][j]; 00663 } 00664 } 00665 00666 } 00667 matrix_destroy (&tmp); 00668 return (1); 00669 } |
|
Use Gaussian elimination to calculate inverse of matrix a, with diagonal scaling applied for stability. Result is matrix ainv. Definition at line 677 of file matrix_f.c. References a, matrix::cols, matrix::elts, flops, free, i, malloc, matrix_destroy(), matrix_equate(), matrix_error(), matrix_initialize(), matrix_inverse(), and matrix::rows. Referenced by calc_glt_matrix(), and calc_matrices().
00678 { 00679 matrix atmp; 00680 register int i, j, n; 00681 register double *diag ; 00682 int mir ; 00683 00684 if (a.rows != a.cols) 00685 matrix_error ("Illegal dimensions for matrix inversion"); 00686 00687 matrix_initialize (&atmp); 00688 00689 n = a.rows; 00690 matrix_equate (a, &atmp); 00691 diag = (double *)malloc( sizeof(double)*n ) ; 00692 for( i=0 ; i < n ; i++ ){ 00693 diag[i] = fabs(atmp.elts[i][i]) ; 00694 if( diag[i] == 0.0 ) diag[i] = 1.0 ; /* shouldn't happen? */ 00695 diag[i] = 1.0 / sqrt(diag[i]) ; 00696 } 00697 00698 for( i=0 ; i < n ; i++ ) /* scale a */ 00699 for( j=0 ; j < n ; j++ ) 00700 atmp.elts[i][j] *= diag[i]*diag[j] ; 00701 00702 mir = matrix_inverse( atmp , ainv ) ; /* invert */ 00703 00704 for( i=0 ; i < n ; i++ ) /* scale inverse */ 00705 for( j=0 ; j < n ; j++ ) 00706 ainv->elts[i][j] *= diag[i]*diag[j] ; 00707 00708 matrix_destroy (&atmp); free((void *)diag) ; 00709 return (mir); 00710 } |
|
Multiply matrix a by matrix b. Result is matrix c. Definition at line 535 of file matrix_f.c. References a, c, matrix::cols, cols, matrix::elts, flops, i, matrix_create(), matrix_error(), matrix::rows, and rows. Referenced by analyze_results(), calc_glt_matrix(), calc_linear_regression(), calc_matrices(), ComputeD0(), ComputeDeltaTau(), ComputeHpHm(), ComputeNewD(), matrix_sqrt(), and poly_field().
00536 { 00537 int rows, cols; 00538 register int i, j, k; 00539 register float sum ; 00540 00541 if (a.cols != b.rows) 00542 matrix_error ("Incompatible dimensions for matrix multiplication"); 00543 00544 rows = a.rows; 00545 cols = b.cols; 00546 00547 matrix_create (rows, cols, c); 00548 00549 for (i = 0; i < rows; i++) 00550 for (j = 0; j < cols; j++) 00551 { 00552 sum = 0.0 ; 00553 for (k = 0; k < a.cols; k++) 00554 sum += a.elts[i][k] * b.elts[k][j]; 00555 c->elts[i][j] = sum; 00556 } 00557 } |
|
Compute the L_infinity norm of a matrix: the max absolute row sum. Definition at line 1197 of file matrix_f.c. References a, matrix::cols, cols, matrix::elts, flops, i, matrix::rows, and rows.
|
|
Print contents of matrix m. Definition at line 198 of file matrix_f.c. References matrix::cols, cols, matrix::elts, i, matrix::rows, rows, and zork. Referenced by calculate_results(), matrix_sprint(), and read_glt_matrix().
00199 { 00200 int i, j; 00201 int rows, cols; 00202 float val ; 00203 int ipr ; 00204 00205 rows = m.rows; 00206 cols = m.cols; 00207 00208 for( i=0 ; i < rows ; i++ ){ 00209 for( j=0 ; j < cols ; j++ ){ 00210 val = (int)m.elts[i][j] ; 00211 if( val != m.elts[i][j] || fabs(val) > 9.0f ) goto zork ; 00212 } 00213 } 00214 zork: 00215 ipr = (i==rows && j==cols) ; 00216 00217 for (i = 0; i < rows; i++) 00218 { 00219 for (j = 0; j < cols; j++) 00220 if( ipr ) printf (" %2d" , (int)m.elts[i][j]); 00221 else printf (" %10.4g", m.elts[i][j]); 00222 printf (" \n"); 00223 } 00224 printf (" \n"); fflush(stdout); 00225 } |
|
Given MxN matrix X, return the NxN matrix [ T ]-1 [ T ]-1 T [ X X ] and the NxM matrix [ X X ] X ----------------------------------------------------------------------------- Definition at line 1313 of file matrix_f.c. References calloc, matrix::cols, matrix::elts, flops, free, matrix_create(), matrix::rows, and svd_double(). Referenced by calc_matrices(), and Form_R_Matrix().
01314 { 01315 int m = X.rows , n = X.cols , ii,jj,kk ; 01316 double *amat , *umat , *vmat , *sval , *xfac , smax,del,sum ; 01317 01318 if( m < 1 || n < 1 || m < n || (XtXinv == NULL && XtXinvXt == NULL) ) return; 01319 01320 amat = (double *)calloc( sizeof(double),m*n ) ; /* input matrix */ 01321 umat = (double *)calloc( sizeof(double),m*n ) ; /* left singular vectors */ 01322 vmat = (double *)calloc( sizeof(double),n*n ) ; /* right singular vectors */ 01323 sval = (double *)calloc( sizeof(double),n ) ; /* singular values */ 01324 xfac = (double *)calloc( sizeof(double),n ) ; /* column norms of [a] */ 01325 01326 #define A(i,j) amat[(i)+(j)*m] 01327 #define U(i,j) umat[(i)+(j)*m] 01328 #define V(i,j) vmat[(i)+(j)*n] 01329 01330 /* copy input matrix into amat */ 01331 01332 for( ii=0 ; ii < m ; ii++ ) 01333 for( jj=0 ; jj < n ; jj++ ) A(ii,jj) = X.elts[ii][jj] ; 01334 01335 /* scale each column to have norm 1 */ 01336 01337 for( jj=0 ; jj < n ; jj++ ){ 01338 sum = 0.0 ; 01339 for( ii=0 ; ii < m ; ii++ ) sum += A(ii,jj)*A(ii,jj) ; 01340 if( sum > 0.0 ) sum = 1.0/sqrt(sum) ; 01341 xfac[jj] = sum ; 01342 for( ii=0 ; ii < m ; ii++ ) A(ii,jj) *= sum ; 01343 } 01344 01345 /* compute SVD of scaled matrix */ 01346 01347 svd_double( m , n , amat , sval , umat , vmat ) ; 01348 01349 free((void *)amat) ; /* done with this */ 01350 01351 /* find largest singular value */ 01352 01353 smax = sval[0] ; 01354 for( ii=1 ; ii < n ; ii++ ) 01355 if( sval[ii] > smax ) smax = sval[ii] ; 01356 01357 if( smax <= 0.0 ){ /* this is bad */ 01358 free((void *)xfac); free((void *)sval); 01359 free((void *)vmat); free((void *)umat); return; 01360 } 01361 01362 for( ii=0 ; ii < n ; ii++ ) 01363 if( sval[ii] < 0.0 ) sval[ii] = 0.0 ; 01364 01365 #ifdef FLOATIZE 01366 #define PSINV_EPS 1.e-8 01367 #else 01368 #define PSINV_EPS 1.e-16 01369 #endif 01370 01371 /* "reciprocals" of singular values: 1/s is actually s/(s^2+del) */ 01372 01373 del = PSINV_EPS * smax*smax ; 01374 for( ii=0 ; ii < n ; ii++ ) 01375 sval[ii] = sval[ii] / ( sval[ii]*sval[ii] + del ) ; 01376 01377 /* create pseudo-inverse */ 01378 01379 if( XtXinvXt != NULL ){ 01380 matrix_create( n , m , XtXinvXt ) ; 01381 for( ii=0 ; ii < n ; ii++ ){ 01382 for( jj=0 ; jj < m ; jj++ ){ 01383 sum = 0.0 ; 01384 for( kk=0 ; kk < n ; kk++ ) 01385 sum += sval[kk] * V(ii,kk) * U(jj,kk) ; 01386 XtXinvXt->elts[ii][jj] = sum * xfac[ii] ; 01387 } 01388 } 01389 } 01390 01391 if( XtXinv != NULL ){ 01392 matrix_create( n , n , XtXinv ) ; 01393 for( ii=0 ; ii < n ; ii++ ) sval[ii] = sval[ii] * sval[ii] ; 01394 matrix_create( n , n , XtXinv ) ; 01395 for( ii=0 ; ii < n ; ii++ ){ 01396 for( jj=0 ; jj < n ; jj++ ){ 01397 sum = 0.0 ; 01398 for( kk=0 ; kk < n ; kk++ ) 01399 sum += sval[kk] * V(ii,kk) * V(jj,kk) ; 01400 XtXinv->elts[ii][jj] = sum * xfac[ii] * xfac[jj] ; 01401 } 01402 } 01403 } 01404 01405 free((void *)xfac); free((void *)sval); 01406 free((void *)vmat); free((void *)umat); return; 01407 } |
|
Multiply matrix a by scalar constant k. Result is matrix c. Definition at line 565 of file matrix_f.c. References a, c, cols, matrix::cols, matrix::elts, i, matrix_create(), rows, and matrix::rows. Referenced by matrix_sqrt().
|
|
Return the eigenvalues of matrix X-transpose X. The output points to a vector of doubles, of length X.cols. This should be free()-ed when you are done with it. ----------------------------------------------------------------------------- Definition at line 1271 of file matrix_f.c. References a, matrix::cols, matrix::elts, flops, free, i, malloc, matrix::rows, and symeigval_double(). Referenced by calculate_results().
01272 { 01273 int i,j,k , M=X.rows , N=X.cols ; 01274 double *a , *e , sum ; 01275 01276 a = (double *) malloc( sizeof(double)*N*N ) ; 01277 e = (double *) malloc( sizeof(double)*N ) ; 01278 01279 for( i=0 ; i < N ; i++ ){ 01280 for( j=0 ; j <= i ; j++ ){ 01281 sum = 0.0 ; 01282 for( k=0 ; k < M ; k++ ) sum += X.elts[k][i] * X.elts[k][j] ; 01283 a[j+N*i] = sum ; 01284 if( j < i ) a[i+N*j] = sum ; 01285 } 01286 } 01287 01288 for( i=0 ; i < N ; i++ ){ 01289 if( a[i+N*i] > 0.0 ) e[i] = 1.0 / sqrt(a[i+N*i]) ; 01290 else e[i] = 1.0 ; 01291 } 01292 01293 for( i=0 ; i < N ; i++ ){ 01294 for( j=0 ; j < N ; j++ ) a[j+N*i] *= sqrt(e[i]*e[j]) ; 01295 } 01296 01297 symeigval_double( N , a , e ) ; 01298 free( (void *)a ) ; 01299 return e ; 01300 } |
|
Print label and contents of matrix m. Definition at line 233 of file matrix_f.c. References matrix_print(). Referenced by calc_covariance(), calculate_results(), markov_array(), matrix_inverse(), and poly_field().
00234 { 00235 printf ("%s \n", s); 00236 00237 matrix_print (m); 00238 } |
|
Calculate square root of symmetric positive definite matrix a. Result is matrix s. Definition at line 719 of file matrix_f.c. References a, matrix::cols, matrix::elts, i, matrix_add(), matrix_destroy(), matrix_equate(), matrix_error(), matrix_identity(), matrix_initialize(), matrix_inverse(), matrix_multiply(), matrix_scale(), matrix_subtract(), and matrix::rows. Referenced by calc_covariance().
00720 { 00721 const int MAX_ITER = 100; 00722 int n; 00723 int ok; 00724 int iter; 00725 register float sse, psse; 00726 register int i, j; 00727 matrix x, xinv, axinv, xtemp, error; 00728 00729 matrix_initialize (&x); 00730 matrix_initialize (&xinv); 00731 matrix_initialize (&axinv); 00732 matrix_initialize (&xtemp); 00733 matrix_initialize (&error); 00734 00735 00736 if (a.rows != a.cols) 00737 matrix_error ("Illegal dimensions for matrix square root"); 00738 00739 00740 n = a.rows; 00741 matrix_identity (n, &x); 00742 00743 00744 psse = 1.0e+30; 00745 for (iter = 0; iter < MAX_ITER; iter++) 00746 { 00747 ok = matrix_inverse (x, &xinv); 00748 if (! ok) return (0); 00749 matrix_multiply (a, xinv, &axinv); 00750 matrix_add (x, axinv, &xtemp); 00751 matrix_scale (0.5, xtemp, &x); 00752 00753 matrix_multiply (x, x, &xtemp); 00754 matrix_subtract (a, xtemp, &error); 00755 sse = 0.0; 00756 for (i = 0; i < n; i++) 00757 for (j = 0; j < n; j++) 00758 sse += error.elts[i][j] * error.elts[i][j] ; 00759 00760 if (sse >= psse) break; 00761 00762 psse = sse; 00763 } 00764 00765 if (iter == MAX_ITER) return (0); 00766 00767 matrix_equate (x, s); 00768 00769 matrix_destroy (&x); 00770 matrix_destroy (&xinv); 00771 matrix_destroy (&axinv); 00772 matrix_destroy (&xtemp); 00773 00774 return (1); 00775 00776 } |
|
Subtract matrix b from matrix a. Result is matrix c. Definition at line 511 of file matrix_f.c. References a, c, matrix::cols, cols, matrix::elts, flops, i, matrix_create(), matrix_error(), matrix::rows, and rows. Referenced by calc_glt_matrix(), and matrix_sqrt().
00512 { 00513 register int rows, cols; 00514 register int i, j; 00515 00516 if ((a.rows != b.rows) || (a.cols != b.cols)) 00517 matrix_error ("Incompatible dimensions for matrix subtraction"); 00518 00519 rows = a.rows; 00520 cols = a.cols; 00521 00522 matrix_create (rows, cols, c); 00523 00524 for (i = 0; i < rows; i++) 00525 for (j = 0; j < cols; j++) 00526 c->elts[i][j] = a.elts[i][j] - b.elts[i][j]; 00527 } |
|
Take transpose of matrix a. Result is matrix t. Definition at line 586 of file matrix_f.c. References a, matrix::cols, cols, matrix::elts, i, matrix_create(), matrix::rows, and rows. Referenced by analyze_results(), calc_glt_matrix(), calc_linear_regression(), calc_matrices(), calculate_results(), ComputeNewD(), and poly_field().
|
|
Read an ASCII file as columns, transpose to rows, allow column selectors.
Definition at line 2079 of file mri_read.c.
02080 { 02081 MRI_IMAGE *inim , *outim , *flim ; 02082 char dname[512] , *cpt , *dpt ; 02083 int ii,jj,nx,ny,nts , *ivlist , *ivl , *sslist ; 02084 float *far , *oar ; 02085 02086 ENTRY("mri_read_1D") ; 02087 02088 if( fname == NULL || fname[0] == '\0' || strlen(fname) > 511 ) RETURN(NULL) ; 02089 02090 if( strncmp(fname,"1D:",3) == 0 ){ /* 28 Apr 2003 */ 02091 return mri_1D_fromstring( fname+3 ) ; 02092 } 02093 02094 /*-- split filename and subvector list --*/ 02095 02096 cpt = strstr(fname,"[") ; 02097 dpt = strstr(fname,"{") ; /* 30 Apr 2003: subsampling list */ 02098 02099 if( cpt == fname || dpt == fname ){ /* can't be at start of filename! */ 02100 fprintf(stderr,"*** Illegal filename in mri_read_1D: %s\n",fname) ; 02101 RETURN(NULL) ; 02102 } else { /* got a subvector list */ 02103 strcpy( dname , fname ) ; 02104 if( cpt != NULL ){ ii = cpt-fname; dname[ii] = '\0'; } 02105 if( dpt != NULL ){ ii = dpt-fname; dname[ii] = '\0'; } 02106 } 02107 02108 /*-- read file in --*/ 02109 02110 inim = mri_read_ascii(dname) ; 02111 if( inim == NULL ) RETURN(NULL) ; 02112 flim = mri_transpose(inim) ; mri_free(inim) ; 02113 02114 /*-- get the subvector and subsampling lists, if any --*/ 02115 02116 nx = flim->nx ; ny = flim->ny ; 02117 02118 ivlist = MCW_get_intlist( ny , cpt ) ; /* subvector list */ 02119 sslist = MCW_get_intlist( nx , dpt ) ; /* subsampling list */ 02120 02121 /* if have subvector list, extract those rows into a new image */ 02122 02123 if( ivlist != NULL && ivlist[0] > 0 ){ 02124 nts = ivlist[0] ; /* number of subvectors */ 02125 ivl = ivlist + 1 ; /* start of array of subvectors */ 02126 02127 for( ii=0 ; ii < nts ; ii++ ){ /* check them out */ 02128 if( ivl[ii] < 0 || ivl[ii] >= ny ){ 02129 fprintf(stderr,"*** Out-of-range subvector [list] in mri_read_1D: %s\n",fname) ; 02130 mri_free(flim) ; free(ivlist) ; RETURN(NULL) ; 02131 } 02132 } 02133 02134 outim = mri_new( nx , nts , MRI_float ) ; /* make output image */ 02135 far = MRI_FLOAT_PTR( flim ) ; 02136 oar = MRI_FLOAT_PTR( outim ) ; 02137 02138 for( ii=0 ; ii < nts ; ii++ ) /* copy desired rows */ 02139 memcpy( oar + ii*nx , far + ivl[ii]*nx , sizeof(float)*nx ) ; 02140 02141 mri_free(flim); free(ivlist); flim = outim; ny = nts; 02142 } 02143 02144 /* if have subsampling list, extract those columns into a new image */ 02145 02146 if( sslist != NULL && sslist[0] > 0 ){ 02147 nts = sslist[0] ; /* number of columns to get */ 02148 ivl = sslist + 1 ; /* start of array of column indexes */ 02149 02150 for( ii=0 ; ii < nts ; ii++ ){ /* check them out */ 02151 if( ivl[ii] < 0 || ivl[ii] >= nx ){ 02152 fprintf(stderr,"*** Out-of-range subsampling {list} in mri_read_1D: %s\n",fname) ; 02153 mri_free(flim) ; free(sslist) ; RETURN(NULL) ; 02154 } 02155 } 02156 02157 outim = mri_new( nts , ny , MRI_float ) ; /* make output image */ 02158 far = MRI_FLOAT_PTR( flim ) ; 02159 oar = MRI_FLOAT_PTR( outim ) ; 02160 02161 for( ii=0 ; ii < nts ; ii++ ) /* copy desired columns */ 02162 for( jj=0 ; jj < ny ; jj++ ) 02163 oar[ii+jj*nts] = far[ivl[ii]+jj*nx] ; 02164 02165 mri_free(flim); free(sslist); flim = outim; 02166 } 02167 02168 RETURN(flim) ; 02169 } |
|
Compute SVD of double precision matrix: T [a] = [u] diag[s] [v]
Definition at line 81 of file cs_symeig.c. References a, free, malloc, svd_(), and v. Referenced by DMAT_svd(), main(), matrix_psinv(), and mri_psinv().
00082 { 00083 integer mm,nn , lda,ldu,ldv , ierr ; 00084 doublereal *aa, *ww , *uu , *vv , *rv1 ; 00085 logical matu , matv ; 00086 00087 if( a == NULL || s == NULL || m < 1 || n < 1 ) return ; 00088 00089 mm = m ; 00090 nn = n ; 00091 aa = a ; 00092 lda = m ; 00093 ww = s ; 00094 00095 if( u == NULL ){ 00096 matu = (logical) 0 ; 00097 uu = (doublereal *)malloc(sizeof(double)*m*n) ; 00098 } else { 00099 matu = (logical) 1 ; 00100 uu = u ; 00101 } 00102 ldu = m ; 00103 00104 if( v == NULL ){ 00105 matv = (logical) 0 ; 00106 vv = NULL ; 00107 } else { 00108 matv = (logical) 1 ; 00109 vv = v ; 00110 } 00111 ldv = n ; 00112 00113 rv1 = (double *) malloc(sizeof(double)*n) ; /* workspace */ 00114 00115 (void) svd_( &mm , &nn , &lda , aa , ww , 00116 &matu , &ldu , uu , &matv , &ldv , vv , &ierr , rv1 ) ; 00117 00118 free((void *)rv1) ; 00119 00120 if( u == NULL ) free((void *)uu) ; 00121 return ; 00122 } |
|
Add vector a to vector b. Result is vector c. Definition at line 948 of file matrix_f.c. References a, c, vector::dim, vector::elts, flops, i, matrix_error(), and vector_create_noinit().
|
|
Create vector v by allocating memory and initializing values. Definition at line 808 of file matrix_f.c. References calloc, vector::dim, vector::elts, i, matrix_error(), v, and vector_destroy(). Referenced by calc_covariance(), calc_tcoef(), calculate_results(), cubic_spline(), do_xrestore_stuff(), glt_analysis(), initialize_state_history(), poly_field(), and regression_analysis().
00809 { 00810 register int i; 00811 00812 vector_destroy (v); 00813 00814 if (dim < 0) matrix_error ("Illegal dimensions for new vector"); 00815 00816 v->dim = dim; 00817 if (dim < 1) return; 00818 00819 v->elts = (float *) calloc (sizeof(float) , dim); 00820 if (v->elts == NULL) 00821 matrix_error ("Memory allocation error"); 00822 } |
|
Definition at line 825 of file matrix_f.c. References vector::dim, vector::elts, i, malloc, matrix_error(), v, and vector_destroy(). Referenced by array_to_vector(), column_to_vector(), DWItoDT_tsfunc(), vector_add(), vector_equate(), vector_multiply(), vector_multiply_subtract(), and vector_subtract().
00825 : RWCox */ 00826 { 00827 register int i; 00828 00829 vector_destroy (v); 00830 00831 if (dim < 0) matrix_error ("Illegal dimensions for new vector"); 00832 00833 v->dim = dim; 00834 if (dim < 1) return; 00835 00836 v->elts = (float *) malloc (sizeof(float ) * dim); 00837 if (v->elts == NULL) 00838 matrix_error ("Memory allocation error"); 00839 } |
|
Destroy vector data structure by deallocating memory. Definition at line 796 of file matrix_f.c. References vector::elts, free, v, and vector_initialize(). Referenced by calc_covariance(), calc_linear_regression(), calc_reduced_model(), calc_resids(), calc_response(), calc_sse(), calc_sse_fit(), calculate_results(), cubic_spline(), DWItoDT_tsfunc(), form_clusters(), FreeGlobals(), glt_analysis(), init_delay(), init_regression_analysis(), poly_field(), print_cluster(), regression_analysis(), reset_options(), terminate_program(), vector_create(), and vector_create_noinit().
00797 { 00798 if (v->elts != NULL) free (v->elts); 00799 vector_initialize (v); 00800 } |
|
Calculate dot product of vector a with vector b. Definition at line 1145 of file matrix_f.c. References a, vector::dim, vector::elts, flops, i, and matrix_error(). Referenced by calc_linear_regression(), calc_resids(), and calc_sse().
01146 { 01147 register int i, dim; 01148 register float sum; 01149 register float *aa , *bb ; 01150 01151 if (a.dim != b.dim) 01152 matrix_error ("Incompatible dimensions for vector dot product"); 01153 01154 dim = a.dim; 01155 01156 sum = 0.0f; 01157 aa = a.elts ; bb = b.elts ; 01158 for (i = 0; i < dim; i++) 01159 #if 0 01160 sum += a.elts[i] * b.elts[i]; 01161 #else 01162 sum += aa[i] * bb[i] ; 01163 #endif 01164 01165 return (sum); 01166 } |
|
Calculate dot product of vector a with itself -- 28 Dec 2001, RWCox. Definition at line 1173 of file matrix_f.c. References a, vector::dim, vector::elts, flops, and i. Referenced by calc_sse_fit().
01174 { 01175 register int i, dim; 01176 register float sum; 01177 register float *aa ; 01178 01179 dim = a.dim; 01180 sum = 0.0f; 01181 aa = a.elts ; 01182 for (i = 0; i < dim; i++) 01183 #if 0 01184 sum += a.elts[i] * a.elts[i]; 01185 #else 01186 sum += aa[i] * aa[i] ; 01187 #endif 01188 01189 return (sum); 01190 } |
|
Copy vector a. Result is vector b. Definition at line 875 of file matrix_f.c. References a, vector::dim, vector::elts, i, and vector_create_noinit(). Referenced by regression_analysis().
00876 { 00877 register int i, dim; 00878 00879 dim = a.dim; 00880 00881 vector_create_noinit (dim, b); 00882 00883 #if 0 00884 for (i = 0; i < dim; i++) 00885 b->elts[i] = a.elts[i]; 00886 #else 00887 if( dim > 0 ) 00888 memcpy( b->elts , a.elts , sizeof(float )*dim ) ; /* RWCox */ 00889 #endif 00890 } |
|
Initialize vector data structure. Definition at line 784 of file matrix_f.c. References vector::dim, vector::elts, and v. Referenced by calc_covariance(), calc_linear_regression(), calc_reduced_model(), calc_resids(), calc_response(), calc_sse(), calc_sse_fit(), calculate_results(), cubic_spline(), do_xrestore_stuff(), DWItoDT_tsfunc(), form_clusters(), glt_analysis(), init_delay(), init_regression_analysis(), InitGlobals(), initialize_options(), initialize_state_history(), poly_field(), print_cluster(), regression_analysis(), and vector_destroy().
|
|
Right multiply matrix a by vector b. Result is vector c. Definition at line 996 of file matrix_f.c. References a, c, matrix::cols, cols, vector::dim, dotnum, dotsum, matrix::elts, vector::elts, flops, i, matrix_error(), MATVEC, matrix::rows, rows, and vector_create_noinit(). Referenced by calc_coef(), calc_lcoef(), calc_linear_regression(), calc_rcoef(), calc_resids(), calc_response(), calc_sse(), calc_sse_fit(), cubic_spline(), do_xrestore_stuff(), DWItoDT_tsfunc(), form_clusters(), poly_field(), and print_cluster().
00997 { 00998 register int rows, cols; 00999 register int i, j; 01000 register float sum ; 01001 register float *bb ; 01002 #ifdef DOTP 01003 register float **aa , *cc ; 01004 #else 01005 register float *aa ; 01006 #endif 01007 01008 if (a.cols != b.dim) 01009 matrix_error ("Incompatible dimensions for vector multiplication"); 01010 01011 rows = a.rows; 01012 cols = a.cols; 01013 01014 vector_create_noinit (rows, c); 01015 01016 if( cols <= 0 ){ 01017 for( i=0 ; i < rows ; i++ ) c->elts[i] = 0.0f ; 01018 return ; 01019 } 01020 01021 bb = b.elts ; 01022 01023 #ifdef DOTP 01024 aa = a.elts ; cc = c->elts ; 01025 i = rows%2 ; 01026 if( i == 1 ) DOTP(cols,aa[0],bb,cc) ; 01027 for( ; i < rows ; i+=2 ){ 01028 DOTP(cols,aa[i] ,bb,cc+i ) ; 01029 DOTP(cols,aa[i+1],bb,cc+(i+1)) ; 01030 } 01031 #else 01032 01033 #ifdef UNROLL_VECMUL 01034 if( cols%2 == 0 ){ /* even number of cols */ 01035 for (i = 0; i < rows; i++){ 01036 sum = 0.0f ; aa = a.elts[i] ; 01037 for (j = 0; j < cols; j+=2 ) 01038 sum += aa[j]*bb[j] + aa[j+1]*bb[j+1]; 01039 c->elts[i] = sum ; 01040 } 01041 } else { /* odd number of cols */ 01042 for (i = 0; i < rows; i++){ 01043 aa = a.elts[i] ; sum = aa[0]*bb[0] ; 01044 for (j = 1; j < cols; j+=2 ) 01045 sum += aa[j]*bb[j] + aa[j+1]*bb[j+1]; 01046 c->elts[i] = sum ; 01047 } 01048 } 01049 #else 01050 for (i = 0; i < rows; i++){ 01051 sum = 0.0f ; aa = a.elts[i] ; 01052 for (j = 0; j < cols; j++ ) sum += aa[j]*bb[j] ; 01053 c->elts[i] = sum ; 01054 } 01055 #endif /* UNROLL_VECMUL */ 01056 01057 #endif /* DOTP */ 01058 01059 } |
|
Compute d = c-a*b: a is a matrix; b,c,d are vectors -- RWCox 26 Feb 2002: return value is sum of squares of d vector Definition at line 1067 of file matrix_f.c. References a, c, matrix::cols, cols, vector::dim, dotnum, dotsum, matrix::elts, vector::elts, flops, free, i, malloc, matrix_error(), matrix::rows, rows, and vector_create_noinit(). Referenced by calc_resids(), and calc_sse().
01068 { 01069 register int rows, cols; 01070 register int i, j; 01071 register float *bb ; 01072 #ifdef DOTP 01073 float qsum,sum , **aa , *dd,*cc,*ee ; 01074 #else 01075 register float qsum,sum, *aa ; 01076 #endif 01077 01078 if (a.cols != b.dim || a.rows != c.dim ) 01079 matrix_error ("Incompatible dimensions for vector multiplication-subtraction"); 01080 01081 rows = a.rows; 01082 cols = a.cols; 01083 01084 vector_create_noinit (rows, d); 01085 01086 if( cols <= 0 ){ 01087 qsum = 0.0f ; 01088 for( i=0 ; i < rows ; i++ ){ 01089 d->elts[i] = c.elts[i] ; 01090 qsum += d->elts[i] * d->elts[i] ; 01091 } 01092 return qsum ; 01093 } 01094 01095 qsum = 0.0f ; bb = b.elts ; 01096 01097 #ifdef DOTP 01098 aa = a.elts ; dd = d->elts ; cc = c.elts ; 01099 ee = (float *)malloc(sizeof(float)*rows) ; 01100 i = rows%2 ; 01101 if( i == 1 ) DOTP(cols,aa[0],bb,ee) ; 01102 for( ; i < rows ; i+=2 ){ 01103 DOTP(cols,aa[i] ,bb,ee+i ) ; 01104 DOTP(cols,aa[i+1],bb,ee+(i+1)) ; 01105 } 01106 VSUB(rows,cc,ee,dd) ; 01107 DOTP(rows,dd,dd,&qsum) ; 01108 free((void *)ee) ; 01109 #else 01110 01111 #ifdef UNROLL_VECMUL 01112 if( cols%2 == 0 ){ /* even number */ 01113 for (i = 0; i < rows; i++){ 01114 aa = a.elts[i] ; sum = c.elts[i] ; 01115 for (j = 0; j < cols; j+=2) 01116 sum -= aa[j]*bb[j] + aa[j+1]*bb[j+1]; 01117 d->elts[i] = sum ; qsum += sum*sum ; 01118 } 01119 } else { /* odd number */ 01120 for (i = 0; i < rows; i++){ 01121 aa = a.elts[i] ; sum = c.elts[i] - aa[0]*bb[0] ; 01122 for (j = 1; j < cols; j+=2) 01123 sum -= aa[j]*bb[j] + aa[j+1]*bb[j+1]; 01124 d->elts[i] = sum ; qsum += sum*sum ; 01125 } 01126 } 01127 #else 01128 for (i = 0; i < rows; i++){ 01129 aa = a.elts[i] ; sum = c.elts[i] ; 01130 for (j = 0; j < cols; j++) sum -= aa[j] * bb[j] ; 01131 d->elts[i] = sum ; qsum += sum*sum ; 01132 } 01133 #endif /* UNROLL_VECMUL */ 01134 01135 #endif /* DOTP */ 01136 01137 return qsum ; /* 26 Feb 2003 */ 01138 } |
|
Print contents of vector v. Definition at line 846 of file matrix_f.c. References vector::dim, vector::elts, i, and v. Referenced by print_cluster(), and vector_sprint().
|
|
Print label and contents of vector v. Definition at line 862 of file matrix_f.c. References v, and vector_print(). Referenced by calc_covariance().
00863 { 00864 printf ("%s \n", s); 00865 00866 vector_print (v); 00867 } |
|
Subtract vector b from vector a. Result is vector c. Definition at line 969 of file matrix_f.c. References a, c, vector::dim, vector::elts, flops, i, matrix_error(), and vector_create_noinit(). Referenced by calc_linear_regression(), calc_resids(), calc_sse(), and calc_sse_fit().
00970 { 00971 register int i, dim; 00972 register float *aa,*bb,*cc ; 00973 00974 if (a.dim != b.dim) 00975 matrix_error ("Incompatible dimensions for vector subtraction"); 00976 00977 dim = a.dim; 00978 00979 vector_create_noinit (dim, c); 00980 00981 aa = a.elts ; bb = b.elts ; cc = c->elts ; 00982 for (i = 0; i < dim; i++) 00983 #if 0 00984 c->elts[i] = a.elts[i] - b.elts[i]; 00985 #else 00986 cc[i] = aa[i] - bb[i] ; 00987 #endif 00988 } |
|
Convert vector v into array f. Definition at line 934 of file matrix_f.c. References vector::dim, vector::elts, i, and v. Referenced by calc_reduced_model(), calculate_results(), DWItoDT_tsfunc(), and form_clusters().
|