00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*
00008   This is the header file for matrix.c.
00009 
00010   File:     matrix.h
00011   Author:   B. Douglas Ward
00012   Date:     23 April 1997
00013 
00014   Mod:      Added routines matrix_file_write and matrix_file_read.
00015   Date:     02 July 1999
00016 
00017   Mod:      Added routine for calculating square root of matrix.
00018   Date:     30 September 1999
00019 
00020   Mod:      Added routines matrix_sprint and vector_sprint.
00021   Date:     04 October 1999
00022 
00023   Mod:      Modified matrix_file_read to use mri_read_ascii routine.
00024   Date:     12 January 2000
00025 
00026   Mod:      Changed return type of vector_dot from float to double.
00027   Date:     13 April 2000
00028 
00029   Mod:      Added functions column_to_vector and matrix_extract_rows.
00030   Date:     21 April 2000
00031 
00032   Mod:      Added functions vector_dotself() and vector_multiply_subtract() -- RWCox.
00033   Date:     28 Dec 2001
00034 
00035   Mod:      Used "register" pointers in some places for speed -- RWCox.
00036   Date:     27 Feb 2003
00037 
00038   Mod:      Use one array instead of array of arrays for matrix -- RWCox.
00039   Date:     04 Mar 2005
00040 
00041 */
00042 
00043 /*---------------------------------------------------------------------------*/
00044 /*
00045   Define matrix and vector data structures.
00046 */
00047 
00048 #include "machdep.h"  /* 07 Mar 2005: to get DONT_USE_MATRIX_MAT */
00049 
00050 typedef struct matrix
00051 {
00052   int      rows;
00053   int      cols;
00054   float  ** elts;
00055 #ifndef DONT_USE_MATRIX_MAT
00056   float   * mat ;  /* 04 Mar 2005 */
00057 #endif
00058 }  matrix;
00059 
00060 
00061 typedef struct vector
00062 {
00063   int  dim;
00064   float  * elts;
00065 } vector;
00066 
00067 
00068 /*---------------------------------------------------------------------------*/
00069 /*
00070   Routine to print an error message and stop.
00071 */
00072 
00073 void matrix_error (char * message);
00074 
00075 
00076 /*---------------------------------------------------------------------------*/
00077 /*
00078   Initialize matrix data structure.
00079 */
00080 
00081 void matrix_initialize (matrix * m);
00082 
00083 
00084 /*---------------------------------------------------------------------------*/
00085 /*
00086   Destroy matrix data structure by deallocating memory.
00087 */
00088 
00089 void matrix_destroy (matrix * m);
00090 
00091 
00092 /*---------------------------------------------------------------------------*/
00093 /*
00094   Create matrix data structure by allocating memory and initializing values.
00095 */
00096 
00097 void matrix_create (int rows, int cols, matrix * m);
00098 
00099 
00100 /*---------------------------------------------------------------------------*/
00101 /*
00102   Print contents of matrix m.
00103 */
00104 
00105 void matrix_print (matrix m);
00106 
00107 
00108 /*---------------------------------------------------------------------------*/
00109 /*
00110   Print label and contents of matrix m.
00111 */
00112 
00113 void matrix_sprint (char * s, matrix m);
00114 
00115 
00116 /*---------------------------------------------------------------------------*/
00117 /*
00118   Print contents of matrix m to specified file.
00119 */
00120 
00121 void matrix_file_write (char * filename, matrix m);
00122 
00123 
00124 /*---------------------------------------------------------------------------*/
00125 /*
00126   Manual entry of matrix data.
00127 */
00128 
00129 void matrix_enter (matrix * m);
00130 
00131 
00132 /*---------------------------------------------------------------------------*/
00133 /*
00134   Read contents of matrix m from specified file.
00135   If unable to read matrix from file, or matrix has wrong dimensions:
00136      If error_exit flag is set, then print error message and exit.
00137      Otherwise, return null matrix.
00138 */
00139 
00140 void matrix_file_read (char * filename, int rows, int cols,  matrix * m,
00141                        int error_exit);
00142 
00143 
00144 /*---------------------------------------------------------------------------*/
00145 /*
00146   Convert simple array to matrix structure.
00147 */
00148 
00149 void array_to_matrix (int rows, int cols, float ** f, matrix * m);
00150 
00151 
00152 /*---------------------------------------------------------------------------*/
00153 /*
00154   Make a copy of the first matrix, return copy as the second matrix.
00155 */
00156 
00157 void matrix_equate (matrix a, matrix * b);
00158 
00159 
00160 /*---------------------------------------------------------------------------*/
00161 /*
00162   Extract p columns (specified by list) from matrix a.  Result is matrix b.
00163 */
00164 
00165 void matrix_extract (matrix a, int p, int * list, matrix * b);
00166 
00167 
00168 /*---------------------------------------------------------------------------*/
00169 /*
00170   Extract p rows (specified by list) from matrix a.  Result is matrix b.
00171 */
00172 
00173 void matrix_extract_rows (matrix a, int p, int * list, matrix * b);
00174 
00175 
00176 /*---------------------------------------------------------------------------*/
00177 /*
00178   Create n x n identity matrix.
00179 */
00180 
00181 void matrix_identity (int n, matrix * m);
00182 
00183 
00184 /*---------------------------------------------------------------------------*/
00185 /*
00186   Add matrix a to matrix b.  Result is matrix c.
00187 */
00188 
00189 void matrix_add (matrix a, matrix b, matrix * c);
00190 
00191 
00192 /*---------------------------------------------------------------------------*/
00193 /*
00194   Subtract matrix b from matrix a.  Result is matrix c.
00195 */
00196 
00197 void matrix_subtract (matrix a, matrix b, matrix * c);
00198 
00199 
00200 /*---------------------------------------------------------------------------*/
00201 /*
00202   Multiply matrix a by matrix b.  Result is matrix c.
00203 */
00204 
00205 void matrix_multiply (matrix a, matrix b, matrix * c);
00206 
00207 
00208 /*---------------------------------------------------------------------------*/
00209 /*
00210   Multiply matrix a by scalar constant k.  Result is matrix c.
00211 */
00212 
00213 void matrix_scale (float  k, matrix a, matrix * c);
00214 
00215 
00216 /*---------------------------------------------------------------------------*/
00217 /*
00218   Take transpose of matrix a.  Result is matrix t.
00219 */
00220 
00221 void matrix_transpose (matrix a, matrix * t);
00222 
00223  
00224 /*---------------------------------------------------------------------------*/
00225 /*
00226   Use Gaussian elimination to calculate inverse of matrix a.  Result is 
00227   matrix ainv.
00228 */
00229 
00230 int matrix_inverse (matrix a, matrix * ainv);
00231 int matrix_inverse_dsc (matrix a, matrix * ainv);
00232 
00233 
00234 /*---------------------------------------------------------------------------*/
00235 /*
00236   Calculate square root of symmetric positive definite matrix a.  
00237   Result is matrix s.
00238 */
00239 
00240 int matrix_sqrt (matrix a, matrix * s);
00241 
00242 
00243 /*---------------------------------------------------------------------------*/
00244 /*
00245   Initialize vector data structure.
00246 */
00247 
00248 void vector_initialize (vector * v);
00249 
00250 
00251 /*---------------------------------------------------------------------------*/
00252 /*
00253   Destroy vector data structure by deallocating memory.
00254 */
00255 
00256 void vector_destroy (vector * v);
00257 
00258 
00259 /*---------------------------------------------------------------------------*/
00260 /*
00261   Create vector v by allocating memory and initializing values.
00262 */
00263 
00264 void vector_create (int dim, vector * v);
00265 
00266 
00267 /*---------------------------------------------------------------------------*/
00268 /*
00269   Print contents of vector v.
00270 */
00271 
00272 void vector_print (vector v);
00273 
00274 
00275 /*---------------------------------------------------------------------------*/
00276 /*
00277   Print label and contents of vector v.
00278 */
00279 
00280 void vector_sprint (char * s, vector v);
00281 
00282 
00283 /*---------------------------------------------------------------------------*/
00284 /*
00285   Copy vector a.  Result is vector b.
00286 */
00287 
00288 void vector_equate (vector a, vector * b);
00289 
00290 
00291 /*---------------------------------------------------------------------------*/
00292 /*
00293   Convert simple array f into vector v.
00294 */
00295 
00296 void array_to_vector (int dim, float * f, vector * v);
00297 
00298 
00299 /*---------------------------------------------------------------------------*/
00300 /*
00301   Convert column c of matrix m into vector v.
00302 */
00303 
00304 void column_to_vector (matrix m, int c, vector * v);
00305 
00306 
00307 /*---------------------------------------------------------------------------*/
00308 /*
00309   Convert vector v into array f.
00310 */
00311 
00312 void vector_to_array (vector v, float * f);
00313 
00314 
00315 /*---------------------------------------------------------------------------*/
00316 /*
00317   Add vector a to vector b.  Result is vector c.
00318 */
00319 
00320 void vector_add (vector a, vector b, vector * c);
00321 
00322 
00323 /*---------------------------------------------------------------------------*/
00324 /*
00325   Subtract vector b from vector a.  Result is vector c.
00326 */
00327 
00328 void vector_subtract (vector a, vector b, vector * c);
00329 
00330 
00331 /*---------------------------------------------------------------------------*/
00332 /*
00333   Right multiply matrix a by vector b.  Result is vector c.
00334 */
00335 
00336 void vector_multiply (matrix a, vector b, vector * c);
00337 
00338 /*---------------------------------------------------------------------------*/
00339 /*
00340   Right multiply matrix a by vector b, then subtract c.  Result is vector d.
00341   Also returns sum of squares of elements of d.
00342 */
00343 
00344 float  vector_multiply_subtract (matrix a, vector b, vector c, vector * d) ;
00345 
00346 /*---------------------------------------------------------------------------*/
00347 /*
00348   Calculate dot product of vector a with vector b. 
00349 */
00350 
00351 float  vector_dot (vector a, vector b);
00352 
00353 float  vector_dotself (vector a);  /* 28 Dec 2001: RWCox */
00354 
00355 /*---------------------------------------------------------------------------*/
00356 
00357 float matrix_norm( matrix a ) ;    /* 03 Mar 2003: RWCox */
00358 
00359 int * matrix_check_columns( matrix a , double eps ) ; /* 14 Jul 2004: RWCox */
00360 
00361 double * matrix_singvals( matrix X ) ; /* 14 Jul 2004 */
00362 
00363 void matrix_psinv( matrix X , matrix *XtXinv , matrix *XtXinvXt ) ;  /* 19 Jul 2004 */