![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 13142
- Timestamp:
- 08/22/12 17:21:15 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
r13056 r13142 17 17 18 18 /*FUNCTION TripleMultiply {{{*/ 19 int TripleMultiply( 19 int TripleMultiply(IssmDouble* a, int nrowa, int ncola, int itrna, IssmDouble* b, int nrowb, int ncolb, int itrnb, IssmDouble* c, int nrowc, int ncolc, int itrnc, IssmDouble* d, int iaddd){ 20 20 /*TripleMultiply Perform triple matrix product a*b*c+d.*/ 21 22 int idima,idimb,idimc,idimd;23 IssmDouble *dtemp;24 25 /* set up dimensions for triple product */26 27 if (!itrna) 21 22 int idima,idimb,idimc,idimd; 23 IssmDouble *dtemp; 24 25 /* set up dimensions for triple product */ 26 27 if (!itrna){ 28 28 idima=nrowa; 29 29 idimb=ncola; 30 30 } 31 else 31 else{ 32 32 idima=ncola; 33 33 idimb=nrowa; 34 34 } 35 35 36 if (!itrnb) { 37 if (nrowb != idimb) { 38 _error_("Matrix A and B inner vectors not equal size."); 39 } 36 if (!itrnb){ 37 if (nrowb != idimb) _error_("Matrix A and B inner vectors not equal size."); 40 38 idimc=ncolb; 41 39 } 42 else { 43 if (ncolb != idimb) { 44 _error_("Matrix A and B inner vectors not equal size."); 45 } 40 else{ 41 if (ncolb != idimb) _error_("Matrix A and B inner vectors not equal size."); 46 42 idimc=nrowb; 47 43 } 48 44 49 45 if (!itrnc) { 50 if (nrowc != idimc) { 51 _error_("Matrix B and C inner vectors not equal size."); 52 } 46 if (nrowc != idimc) _error_("Matrix B and C inner vectors not equal size."); 53 47 idimd=ncolc; 54 48 } 55 else { 56 if (ncolc != idimc) { 57 _error_("Matrix B and C inner vectors not equal size."); 58 } 49 else{ 50 if (ncolc != idimc) _error_("Matrix B and C inner vectors not equal size."); 59 51 idimd=nrowc; 60 52 } 61 53 62 /* perform the matrix triple product in the order that minimizes the 63 number of multiplies and the temporary space used, noting that 64 (a*b)*c requires ac(b+d) multiplies and ac IssmDoubles, and a*(b*c) 65 requires bd(a+c) multiplies and bd IssmDoubles (both are the same for 66 a symmetric triple product) */ 67 68 /* multiply (a*b)*c+d */ 69 54 /* perform the matrix triple product in the order that minimizes the 55 number of multiplies and the temporary space used, noting that 56 (a*b)*c requires ac(b+d) multiplies and ac IssmDoubles, and a*(b*c) 57 requires bd(a+c) multiplies and bd IssmDoubles (both are the same for 58 a symmetric triple product) */ 59 60 /* multiply (a*b)*c+d */ 70 61 if (idima*idimc*(idimb+idimd) <= idimb*idimd*(idima+idimc)) { 71 62 dtemp=xNew<IssmDouble>(idima*idimc); … … 76 67 } 77 68 78 /* multiply a*(b*c)+d */ 79 80 else { 69 /* multiply a*(b*c)+d */ 70 else{ 81 71 dtemp=xNew<IssmDouble>(idimb*idimd); 82 72 … … 89 79 }/*}}}*/ 90 80 /*FUNCTION MatrixMuliply {{{*/ 91 int MatrixMultiply( 81 int MatrixMultiply(IssmDouble* a, int nrowa, int ncola, int itrna, IssmDouble* b, int nrowb, int ncolb, int itrnb, IssmDouble* c, int iaddc ){ 92 82 /*MatrixMultiply Perform matrix multiplication a*b+c.*/ 93 83 int noerr=1; … … 95 85 int nrowc,ncolc,iinca,jinca,iincb,jincb,ntrma,ntrmb,nterm; 96 86 97 /* set up dimensions and increments for matrix a */87 /* set up dimensions and increments for matrix a */ 98 88 if (!itrna) { 99 89 nrowc=nrowa; … … 109 99 } 110 100 111 /* set up dimensions and increments for matrix b */101 /* set up dimensions and increments for matrix b */ 112 102 if (!itrnb) { 113 103 ncolc=ncolb; … … 123 113 } 124 114 125 if (ntrma != ntrmb) { 126 _error_("Matrix A and B inner vectors not equal size"); 127 noerr=0; 128 return noerr; 129 } 130 else 131 nterm=ntrma; 132 133 /* zero matrix c, if not being added to product */ 134 135 if (!iaddc) 136 for (i=0; i<nrowc*ncolc; i++) 137 *(c+i)=0.; 138 139 /* perform the matrix multiplication */ 140 115 if (ntrma != ntrmb) _error_("Matrix A and B inner vectors not equal size"); 116 117 nterm=ntrma; 118 119 /* zero matrix c, if not being added to product */ 120 if (!iaddc) for (i=0;i<nrowc*ncolc;i++) c[i]=0.; 121 122 /* perform the matrix multiplication */ 141 123 iptc=0; 142 for (i=0; i<nrowc; i++) 143 for (j=0; j<ncolc; j++) 124 for (i=0; i<nrowc; i++){ 125 for (j=0; j<ncolc; j++){ 144 126 ipta=i*iinca; 145 127 iptb=j*jincb; 146 128 147 for (k=0; k<nterm; k++) 148 *(c+iptc)+=*(a+ipta)**(b+iptb);129 for (k=0; k<nterm; k++){ 130 c[iptc]+=a[ipta]*b[iptb]; 149 131 ipta+=jinca; 150 132 iptb+=iincb; 151 133 } 152 153 134 iptc++; 154 135 } … … 159 140 /*FUNCTION MatrixInverse {{{*/ 160 141 int MatrixInverse( IssmDouble* a, int ndim, int nrow, IssmDouble* b, int nvec, IssmDouble* pdet ){ 161 /* MatrixInverse Perform matrix inversion and linear equation solution.162 163 This function uses Gaussian elimination on the original matrix164 augmented by an identity matrix of the same size to calculate165 the inverse (see for example, "Modern Methods of Engineering166 Computation", Sec. 6.4). By noting how the matrices are167 unpopulated and repopulated, the calculation may be done in place.168 169 Gaussian elimination is inherently inefficient, and so this is170 intended for small matrices. */142 /* MatrixInverse Perform matrix inversion and linear equation solution. 143 144 This function uses Gaussian elimination on the original matrix 145 augmented by an identity matrix of the same size to calculate 146 the inverse (see for example, "Modern Methods of Engineering 147 Computation", Sec. 6.4). By noting how the matrices are 148 unpopulated and repopulated, the calculation may be done in place. 149 150 Gaussian elimination is inherently inefficient, and so this is 151 intended for small matrices. */ 171 152 int noerr=1; 172 153 int i,j,k,ipt,jpt,irow,icol,ipiv,ncol; … … 184 165 _assert_(!(ndim==3 && nrow==3)); 185 166 186 /* initialize local variables and arrays */167 /* initialize local variables and arrays */ 187 168 188 169 ncol=nrow; … … 192 173 pindx =xNew<int>(nrow); 193 174 194 /* loop over the rows/columns of the matrix */175 /* loop over the rows/columns of the matrix */ 195 176 196 177 for (i=0; i<nrow; i++) { 197 178 198 /* search for pivot, finding the term with the greatest magnitude199 in the rows/columns not yet used */179 /* search for pivot, finding the term with the greatest magnitude 180 in the rows/columns not yet used */ 200 181 201 182 pivot=0.; 202 183 for (j=0; j<nrow; j++) 203 204 205 206 207 208 209 210 184 if (!pindx[j]) 185 for (k=0; k<ncol; k++) 186 if (!pindx[k]) 187 if (fabs(a[j*ndim+k]) > fabs(pivot)) { 188 irow=j; 189 icol=k; 190 pivot=a[j*ndim+k]; 191 } 211 192 212 193 if (fabs(pivot) < DBL_EPSILON) { … … 225 206 pindx[ipiv]++; 226 207 227 // _printf_(true,"pivot for i=%d: irow=%d, icol=%d, pindx[%d]=%d\n",228 // i,irow,icol,ipiv,pindx[ipiv]);229 230 /* switch rows to put pivot element on diagonal, noting that the231 column stays the same and the determinant changes sign */208 // _printf_(true,"pivot for i=%d: irow=%d, icol=%d, pindx[%d]=%d\n", 209 // i,irow,icol,ipiv,pindx[ipiv]); 210 211 /* switch rows to put pivot element on diagonal, noting that the 212 column stays the same and the determinant changes sign */ 232 213 233 214 if (irow != icol) { 234 // _printf_(true,"row switch for i=%d: irow=%d, icol=%d\n",235 // i,irow,icol);215 // _printf_(true,"row switch for i=%d: irow=%d, icol=%d\n", 216 // i,irow,icol); 236 217 237 218 ipt=irow*ndim; … … 254 235 } 255 236 256 /* divide pivot row by pivot element, noting that the original257 matrix will have 1 on the diagonal, which will be discarded,258 and the augmented matrix will start with 1 from the identity259 matrix and then have 1/pivot, which is part of the inverse. */237 /* divide pivot row by pivot element, noting that the original 238 matrix will have 1 on the diagonal, which will be discarded, 239 and the augmented matrix will start with 1 from the identity 240 matrix and then have 1/pivot, which is part of the inverse. */ 260 241 261 242 a[ipiv*ndim+ipiv]=1.; … … 263 244 ipt=ipiv*ndim; 264 245 for (k=0; k<ncol; k++) 265 246 a[ipt+k]/=pivot; 266 247 267 248 ipt=ipiv*nvec; 268 249 for (k=0; k<nvec; k++) 269 270 271 /* reduce non-pivot rows such that they will have 0 in the pivot272 column, which will be discarded, and the augmented matrix will273 start with 0 from the identity matrix and then have non-zero274 in the corresponding column, which is part of the inverse.275 only one column of the augmented matrix is populated at a time,276 which corresponds to the only column of the original matrix277 being zeroed, so that the inverse may be done in place. */250 b[ipt+k]/=pivot; 251 252 /* reduce non-pivot rows such that they will have 0 in the pivot 253 column, which will be discarded, and the augmented matrix will 254 start with 0 from the identity matrix and then have non-zero 255 in the corresponding column, which is part of the inverse. 256 only one column of the augmented matrix is populated at a time, 257 which corresponds to the only column of the original matrix 258 being zeroed, so that the inverse may be done in place. */ 278 259 279 260 for (j=0; j<nrow; j++) { … … 287 268 jpt=ipiv*ndim; 288 269 for (k=0; k<ncol; k++) 289 270 a[ipt+k]-=dtemp*a[jpt+k]; 290 271 291 272 ipt=j *nvec; 292 273 jpt=ipiv*nvec; 293 274 for (k=0; k<nvec; k++) 294 295 } 296 } 297 298 /* for a diagonal matrix, the determinant is the product of the299 diagonal terms, and so it may be accumulated from the pivots,300 noting that switching rows changes the sign as above */275 b[ipt+k]-=dtemp*b[jpt+k]; 276 } 277 } 278 279 /* for a diagonal matrix, the determinant is the product of the 280 diagonal terms, and so it may be accumulated from the pivots, 281 noting that switching rows changes the sign as above */ 301 282 302 283 det*=pivot; 303 284 } 304 285 305 /* switch columns back in reverse order, noting that a row switch306 in the original matrix corresponds to a column switch in the307 inverse matrix */286 /* switch columns back in reverse order, noting that a row switch 287 in the original matrix corresponds to a column switch in the 288 inverse matrix */ 308 289 309 290 for (i=0; i<nrow; i++) { … … 314 295 icol=pivrc2[j]; 315 296 316 // _printf_(true,"column switch back for j=%d: irow=%d, icol=%d\n",317 // j,irow,icol);297 // _printf_(true,"column switch back for j=%d: irow=%d, icol=%d\n", 298 // j,irow,icol); 318 299 319 300 ipt=0; … … 363 344 364 345 /*det = a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)*/ 365 346 *Adet= A[0]*A[4]*A[8]-A[0]*A[5]*A[7]-A[3]*A[1]*A[8]+A[3]*A[2]*A[7]+A[6]*A[1]*A[5]-A[6]*A[2]*A[4]; 366 347 } 367 348 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.
