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