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().
|
