1 #ifndef MATRIX_F_INCLUDED 2 #define MATRIX_F_INCLUDED 3 /***************************************************************************** 4 Major portions of this software are copyrighted by the Medical College 5 of Wisconsin, 1994-2000, and are released under the Gnu General Public 6 License, Version 2. See the file README.Copyright for details. 7 ******************************************************************************/ 8 9 /* 10 This is the header file for matrix.c. 11 12 File: matrix.h 13 Author: B. Douglas Ward 14 Date: 23 April 1997 15 16 Mod: Added routines matrix_file_write and matrix_file_read. 17 Date: 02 July 1999 18 19 Mod: Added routine for calculating square root of matrix. 20 Date: 30 September 1999 21 22 Mod: Added routines matrix_sprint and vector_sprint. 23 Date: 04 October 1999 24 25 Mod: Modified matrix_file_read to use mri_read_ascii routine. 26 Date: 12 January 2000 27 28 Mod: Changed return type of vector_dot from float to double. 29 Date: 13 April 2000 30 31 Mod: Added functions column_to_vector and matrix_extract_rows. 32 Date: 21 April 2000 33 34 Mod: Added functions vector_dotself() and vector_multiply_subtract() -- RWCox. 35 Date: 28 Dec 2001 36 37 Mod: Used "register" pointers in some places for speed -- RWCox. 38 Date: 27 Feb 2003 39 40 Mod: Use one array instead of array of arrays for matrix -- RWCox. 41 Date: 04 Mar 2005 42 43 */ 44 45 /*---------------------------------------------------------------------------*/ 46 /* 47 Define matrix and vector data structures. 48 */ 49 50 #include "machdep.h" /* 07 Mar 2005: to get DONT_USE_MATRIX_MAT */ 51 52 typedef struct matrix 53 { 54 int rows; 55 int cols; 56 float ** elts; 57 #ifndef DONT_USE_MATRIX_MAT 58 float * mat ; /* 04 Mar 2005 */ 59 #endif 60 } matrix; 61 62 63 typedef struct vector 64 { 65 int dim; 66 float * elts; 67 } vector; 68 69 70 #undef ISVALID_MATRIX 71 #define ISVALID_MATRIX(m) ((m).rows > 0 && (m).cols > 0 && (m).elts != NULL) 72 73 #undef ISVALID_VECTOR 74 #define ISVALID_VECTOR(v) ((v).dim > 0 && (v).elts != NULL) 75 76 /*---------------------------------------------------------------------------*/ 77 /* 78 Routine to print an error message and stop. 79 */ 80 81 void matrix_error (char * message); 82 83 84 /*---------------------------------------------------------------------------*/ 85 /* 86 Initialize matrix data structure. 87 */ 88 89 void matrix_initialize (matrix * m); 90 91 92 /*---------------------------------------------------------------------------*/ 93 /* 94 Destroy matrix data structure by deallocating memory. 95 */ 96 97 void matrix_destroy (matrix * m); 98 99 100 /*---------------------------------------------------------------------------*/ 101 /* 102 Create matrix data structure by allocating memory and initializing values. 103 */ 104 105 void matrix_create (int rows, int cols, matrix * m); 106 107 108 /*---------------------------------------------------------------------------*/ 109 /* 110 Print contents of matrix m. 111 */ 112 113 void matrix_print (matrix m); 114 115 116 /*---------------------------------------------------------------------------*/ 117 /* 118 Print label and contents of matrix m. 119 */ 120 121 void matrix_sprint (char * s, matrix m); 122 123 124 /*---------------------------------------------------------------------------*/ 125 /* 126 Print contents of matrix m to specified file. 127 */ 128 129 void matrix_file_write (char * filename, matrix m); 130 131 132 /*---------------------------------------------------------------------------*/ 133 /* 134 Manual entry of matrix data. 135 */ 136 137 void matrix_enter (matrix * m); 138 139 140 /*---------------------------------------------------------------------------*/ 141 /* 142 Read contents of matrix m from specified file. 143 If unable to read matrix from file, or matrix has wrong dimensions: 144 If error_exit flag is set, then print error message and exit. 145 Otherwise, return null matrix. 146 */ 147 148 void matrix_file_read (char * filename, int rows, int cols, matrix * m, 149 int error_exit); 150 151 152 /*---------------------------------------------------------------------------*/ 153 /* 154 Convert simple array to matrix structure. 155 */ 156 157 void array_to_matrix (int rows, int cols, float ** f, matrix * m); 158 159 160 /*---------------------------------------------------------------------------*/ 161 /* 162 Make a copy of the first matrix, return copy as the second matrix. 163 */ 164 165 void matrix_equate (matrix a, matrix * b); 166 167 168 /*---------------------------------------------------------------------------*/ 169 /* 170 Extract p columns (specified by list) from matrix a. Result is matrix b. 171 */ 172 173 void matrix_extract (matrix a, int p, int * list, matrix * b); 174 175 #define matrix_extract_cols matrix_extract 176 177 /*---------------------------------------------------------------------------*/ 178 /* 179 Extract p rows (specified by list) from matrix a. Result is matrix b. 180 */ 181 182 void matrix_extract_rows (matrix a, int p, int * list, matrix * b); 183 184 185 /*---------------------------------------------------------------------------*/ 186 /* 187 Add columns that are all 0 except for a single 1, as indicated. 188 */ 189 190 void matrix_augment_01_columns( matrix a, int nadd, int *addlist, matrix *b ); 191 192 193 /*---------------------------------------------------------------------------*/ 194 /* 195 Create n x n identity matrix. 196 */ 197 198 void matrix_identity (int n, matrix * m); 199 200 201 /*---------------------------------------------------------------------------*/ 202 /* 203 Add matrix a to matrix b. Result is matrix c. 204 */ 205 206 void matrix_add (matrix a, matrix b, matrix * c); 207 208 209 /*---------------------------------------------------------------------------*/ 210 /* 211 Subtract matrix b from matrix a. Result is matrix c. 212 */ 213 214 void matrix_subtract (matrix a, matrix b, matrix * c); 215 216 217 /*---------------------------------------------------------------------------*/ 218 /* 219 Multiply matrix a by matrix b. Result is matrix c. 220 */ 221 222 void matrix_multiply (matrix a, matrix b, matrix * c); 223 224 225 /*---------------------------------------------------------------------------*/ 226 /* 227 Multiply matrix a by scalar constant k. Result is matrix c. 228 */ 229 230 void matrix_scale (float k, matrix a, matrix * c); 231 232 233 /*---------------------------------------------------------------------------*/ 234 /* 235 Take transpose of matrix a. Result is matrix t. 236 */ 237 238 void matrix_transpose (matrix a, matrix * t); 239 240 241 /*---------------------------------------------------------------------------*/ 242 /* 243 Use Gaussian elimination to calculate inverse of matrix a. Result is 244 matrix ainv. 245 */ 246 247 int matrix_inverse (matrix a, matrix * ainv); 248 int matrix_inverse_dsc (matrix a, matrix * ainv); 249 250 251 /*---------------------------------------------------------------------------*/ 252 /* 253 Calculate square root of symmetric positive definite matrix a. 254 Result is matrix s. 255 */ 256 257 int matrix_sqrt (matrix a, matrix * s); 258 259 260 /*---------------------------------------------------------------------------*/ 261 /* 262 Initialize vector data structure. 263 */ 264 265 void vector_initialize (vector * v); 266 267 268 /*---------------------------------------------------------------------------*/ 269 /* 270 Destroy vector data structure by deallocating memory. 271 */ 272 273 void vector_destroy (vector * v); 274 275 276 /*---------------------------------------------------------------------------*/ 277 /* 278 Create vector v by allocating memory and initializing values. 279 */ 280 281 void vector_create (int dim, vector * v); 282 283 284 /*---------------------------------------------------------------------------*/ 285 /* 286 Print contents of vector v. 287 */ 288 289 void vector_print (vector v); 290 291 292 /*---------------------------------------------------------------------------*/ 293 /* 294 Print label and contents of vector v. 295 */ 296 297 void vector_sprint (char * s, vector v); 298 299 300 /*---------------------------------------------------------------------------*/ 301 /* 302 Copy vector a. Result is vector b. 303 */ 304 305 void vector_equate (vector a, vector * b); 306 307 308 /*---------------------------------------------------------------------------*/ 309 /* 310 Convert simple array f into vector v. 311 */ 312 313 void array_to_vector (int dim, float * f, vector * v); 314 315 316 /*---------------------------------------------------------------------------*/ 317 /* 318 Convert column c of matrix m into vector v. 319 */ 320 321 void column_to_vector (matrix m, int c, vector * v); 322 323 324 /*---------------------------------------------------------------------------*/ 325 /* 326 Convert vector v into array f. 327 */ 328 329 void vector_to_array (vector v, float * f); 330 331 332 /*---------------------------------------------------------------------------*/ 333 /* 334 Add vector a to vector b. Result is vector c. 335 */ 336 337 void vector_add (vector a, vector b, vector * c); 338 339 340 /*---------------------------------------------------------------------------*/ 341 /* 342 Subtract vector b from vector a. Result is vector c. 343 */ 344 345 void vector_subtract (vector a, vector b, vector * c); 346 347 348 /*---------------------------------------------------------------------------*/ 349 /* 350 Right multiply matrix a by vector b. Result is vector c. 351 */ 352 353 void vector_multiply (matrix a, vector b, vector * c); 354 void vector_multiply_transpose (matrix a, vector b, vector * c); 355 356 /*---------------------------------------------------------------------------*/ 357 /* 358 Right multiply matrix a by vector b, then subtract c. Result is vector d. 359 Also returns sum of squares of elements of d. 360 */ 361 362 float vector_multiply_subtract (matrix a, vector b, vector c, vector * d) ; 363 364 /*---------------------------------------------------------------------------*/ 365 /* 366 Calculate dot product of vector a with vector b. 367 */ 368 369 float vector_dot (vector a, vector b); 370 371 float vector_dotself (vector a); /* 28 Dec 2001: RWCox */ 372 373 /*---------------------------------------------------------------------------*/ 374 375 float matrix_norm( matrix a ) ; /* 03 Mar 2003: RWCox */ 376 377 int * matrix_check_columns( matrix a , double eps ) ; /* 14 Jul 2004: RWCox */ 378 379 double * matrix_singvals( matrix X ) ; /* 14 Jul 2004 */ 380 381 void matrix_psinv( matrix X , matrix *XtXinv , matrix *XtXinvXt ) ; /* 19 Jul 2004 */ 382 383 int matrix_qrr( matrix X , matrix *R ) ; /* 03 Jul 2008 */ 384 void vector_rr_solve( matrix R , vector b , vector *x ) ; 385 void vector_rrtran_solve( matrix R , vector b , vector *x ) ; 386 387 #endif 388