Ice Sheet System Model  4.18
Code documentation
Functions
matrix.h File Reference

prototypes for matrix.h More...

#include "../Numerics/types.h"

Go to the source code of this file.

Functions

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)
 
int MatrixMultiply (IssmDouble *a, int nrowa, int ncola, int itrna, IssmDouble *b, int nrowb, int ncolb, int itrnb, IssmDouble *c, int iaddc)
 
int MatrixInverse (IssmDouble *a, int ndim, int nrow, IssmDouble *b, int nvec, IssmDouble *pdet)
 
void Matrix2x2Invert (IssmDouble *Ainv, IssmDouble *A)
 
void Matrix2x2Determinant (IssmDouble *Adet, IssmDouble *A)
 
void Matrix2x2Eigen (IssmDouble *plambda1, IssmDouble *plambda2, IssmDouble *pvx, IssmDouble *pvy, IssmDouble a11, IssmDouble a21, IssmDouble a22)
 
void Matrix3x3Invert (IssmDouble *Ainv, IssmDouble *A)
 
void Matrix3x3Determinant (IssmDouble *Adet, IssmDouble *A)
 
IssmDouble Matrix3x3Determinant (IssmDouble a1, IssmDouble a2, IssmDouble a3, IssmDouble b1, IssmDouble b2, IssmDouble b3, IssmDouble c1, IssmDouble c2, IssmDouble c3)
 
void Matrix3x3Solve (IssmDouble *X, IssmDouble *A, IssmDouble *B)
 
void Matrix4x4Adjoint (IssmDouble *Aadj, IssmDouble *A)
 
void Matrix4x4Invert (IssmDouble *Ainv, IssmDouble *A)
 
void Matrix4x4Determinant (IssmDouble *Adet, IssmDouble *A)
 
void Matrix4x4Solve (IssmDouble *X, IssmDouble *A, IssmDouble *B)
 
void newcell (IssmDouble **pcell, IssmDouble newvalue, bool top, int m)
 
IssmDouble cellsum (IssmDouble *cell, int m)
 
void celldelete (IssmDouble **pcell, int m, int *indices, int nind)
 
void cellsplit (IssmDouble **pcell, int m, int i, IssmDouble scale)
 
void cellecho (int numcells, int m,...)
 

Detailed Description

prototypes for matrix.h

Definition in file matrix.h.

Function Documentation

◆ TripleMultiply()

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 
)

Definition at line 20 of file MatrixUtils.cpp.

20  {/*{{{*/
21  /*TripleMultiply Perform triple matrix product a*b*c+d.*/
22 
23  int idima,idimb,idimc,idimd;
24  IssmDouble dtemp_static[600];
25  IssmDouble* dtemp_dynamic = NULL;
26  IssmDouble* dtemp = NULL;
27  _assert_(a && b && c && d);
28 
29  /* set up dimensions for triple product */
30  if (!itrna){
31  idima=nrowa;
32  idimb=ncola;
33  }
34  else{
35  idima=ncola;
36  idimb=nrowa;
37  }
38 
39  if (!itrnb){
40  if (nrowb != idimb) _error_("Matrix A and B inner vectors not equal size.");
41  idimc=ncolb;
42  }
43  else{
44  if (ncolb != idimb) _error_("Matrix A and B inner vectors not equal size.");
45  idimc=nrowb;
46  }
47 
48  if (!itrnc) {
49  if (nrowc != idimc) _error_("Matrix B and C inner vectors not equal size.");
50  idimd=ncolc;
51  }
52  else{
53  if (ncolc != idimc) _error_("Matrix B and C inner vectors not equal size.");
54  idimd=nrowc;
55  }
56 
57  /*Depending on the size of incoming matrices, we might need to use a dynamic allocation*/
58  if(idima*idimc>600){
59  dtemp_dynamic = xNew<IssmDouble>(idima*idimc);
60  dtemp = dtemp_dynamic;
61  }
62  else{
63  dtemp = &dtemp_static[0];
64  }
65 
66  /* perform the matrix triple product in the order that minimizes the
67  number of multiplies and the temporary space used, noting that
68  (a*b)*c requires ac(b+d) multiplies and ac IssmDoubles, and a*(b*c)
69  requires bd(a+c) multiplies and bd IssmDoubles (both are the same for
70  a symmetric triple product) */
71 
72  /* multiply (a*b)*c+d */
73  if (idima*idimc*(idimb+idimd) <= idimb*idimd*(idima+idimc)) {
74  MatrixMultiply(a,nrowa,ncola,itrna,b,nrowb,ncolb,itrnb,dtemp,0);
75  MatrixMultiply(dtemp,idima,idimc,0,c,nrowc,ncolc,itrnc,d,iaddd);
76  }
77 
78  /* multiply a*(b*c)+d */
79  else{
80  MatrixMultiply(b,nrowb,ncolb,itrnb,c,nrowc,ncolc,itrnc,dtemp,0);
81  MatrixMultiply(a,nrowa,ncola,itrna,dtemp,idimb,idimd,0,d,iaddd);
82  }
83 
84  /*Cleanup and return*/
85  xDelete<IssmDouble>(dtemp_dynamic);
86  return 1;
87 }/*}}}*/

◆ MatrixMultiply()

int MatrixMultiply ( IssmDouble *  a,
int  nrowa,
int  ncola,
int  itrna,
IssmDouble *  b,
int  nrowb,
int  ncolb,
int  itrnb,
IssmDouble *  c,
int  iaddc 
)

Definition at line 88 of file MatrixUtils.cpp.

88  {/*{{{*/
89  /*MatrixMultiply Perform matrix multiplication a*b+c.*/
90  int noerr=1;
91  int i,j,k,ipta,iptb,iptc;
92  int nrowc,ncolc,iinca,jinca,iincb,jincb,ntrma,ntrmb,nterm;
93 
94  _assert_(a && b && c);
95 
96  /* set up dimensions and increments for matrix a */
97  if (!itrna) {
98  nrowc=nrowa;
99  ntrma=ncola;
100  iinca=ncola;
101  jinca=1;
102  }
103  else {
104  nrowc=ncola;
105  ntrma=nrowa;
106  iinca=1;
107  jinca=ncola;
108  }
109 
110  /* set up dimensions and increments for matrix b */
111  if (!itrnb) {
112  ncolc=ncolb;
113  ntrmb=nrowb;
114  iincb=ncolb;
115  jincb=1;
116  }
117  else {
118  ncolc=nrowb;
119  ntrmb=ncolb;
120  iincb=1;
121  jincb=ncolb;
122  }
123 
124  if (ntrma != ntrmb) _error_("Matrix A and B inner vectors not equal size");
125 
126  nterm=ntrma;
127 
128  /* zero matrix c, if not being added to product */
129  if (!iaddc) for (i=0;i<nrowc*ncolc;i++) c[i]=0.;
130 
131  /* perform the matrix multiplication */
132  iptc=0;
133  for (i=0; i<nrowc; i++){
134  for (j=0; j<ncolc; j++){
135  ipta=i*iinca;
136  iptb=j*jincb;
137 
138  for (k=0; k<nterm; k++){
139  c[iptc]+=a[ipta]*b[iptb];
140  ipta+=jinca;
141  iptb+=iincb;
142  }
143  iptc++;
144  }
145  }
146 
147  return noerr;
148 }/*}}}*/

◆ MatrixInverse()

int MatrixInverse ( IssmDouble *  a,
int  ndim,
int  nrow,
IssmDouble *  b,
int  nvec,
IssmDouble *  pdet 
)

Definition at line 149 of file MatrixUtils.cpp.

149  {/*{{{*/
150  /* MatrixInverse Perform matrix inversion and linear equation solution.
151 
152  This function uses Gaussian elimination on the original matrix
153  augmented by an identity matrix of the same size to calculate
154  the inverse (see for example, "Modern Methods of Engineering
155  Computation", Sec. 6.4). By noting how the matrices are
156  unpopulated and repopulated, the calculation may be done in place.
157 
158  Gaussian elimination is inherently inefficient, and so this is
159  intended for small matrices. */
160  int noerr=1;
161  int i,j,k,ipt,jpt,irow,icol,ipiv,ncol;
162  int *pivrc1,*pivrc2,*pindx;
163  IssmDouble pivot,det,dtemp;
164 
165  if (!b && nvec) {
166  _error_("No right-hand side for nvec=" << nvec << ".");
167  noerr=0;
168  return noerr;
169  }
170 
171  /*In debugging mode, check that we are not dealing with simple matrices*/
172  _assert_(!(ndim==2 && nrow==2));
173  _assert_(!(ndim==3 && nrow==3));
174 
175  /* initialize local variables and arrays */
176 
177  ncol=nrow;
178  det=1.;
179  pivrc1 =xNew<int>(nrow);
180  pivrc2 =xNew<int>(nrow);
181  pindx =xNew<int>(nrow);
182 
183  /* loop over the rows/columns of the matrix */
184 
185  for (i=0; i<nrow; i++) {
186 
187  /* search for pivot, finding the term with the greatest magnitude
188  in the rows/columns not yet used */
189 
190  pivot=0.;
191  for (j=0; j<nrow; j++)
192  if (!pindx[j])
193  for (k=0; k<ncol; k++)
194  if (!pindx[k])
195  if (fabs(a[j*ndim+k]) > fabs(pivot)) {
196  irow=j;
197  icol=k;
198  pivot=a[j*ndim+k];
199  }
200 
201  if (fabs(pivot) < DBL_EPSILON) {
202  xDelete<int>(pivrc1);
203  xDelete<int>(pivrc2);
204  xDelete<int>(pindx);
205  _error_("Pivot " << pivot << " less than machine epsilon");
206  noerr=0;
207  return noerr;
208  }
209 
210  pivrc1[i]=irow;
211  pivrc2[i]=icol;
212 
213  ipiv=icol;
214  pindx[ipiv]++;
215 
216  // _printf_("pivot for i=" << i << ": irow=" << irow << ", icol=" << icol << ", pindx[" << ipiv << "]=" << pindx[ipiv] << "\n\n");
217 
218  /* switch rows to put pivot element on diagonal, noting that the
219  column stays the same and the determinant changes sign */
220 
221  if (irow != icol) {
222  // _printf_("row switch for i=" << i << ": irow=" << irow << ", icol=" << icol << "\n\n");
223 
224  ipt=irow*ndim;
225  jpt=icol*ndim;
226  for (k=0; k<ncol; k++) {
227  dtemp =a[ipt+k];
228  a[ipt+k]=a[jpt+k];
229  a[jpt+k]=dtemp;
230  }
231 
232  ipt=irow*nvec;
233  jpt=icol*nvec;
234  for (k=0; k<nvec; k++) {
235  dtemp =b[ipt+k];
236  b[ipt+k]=b[jpt+k];
237  b[jpt+k]=dtemp;
238  }
239 
240  det=-det;
241  }
242 
243  /* divide pivot row by pivot element, noting that the original
244  matrix will have 1 on the diagonal, which will be discarded,
245  and the augmented matrix will start with 1 from the identity
246  matrix and then have 1/pivot, which is part of the inverse. */
247 
248  a[ipiv*ndim+ipiv]=1.;
249 
250  ipt=ipiv*ndim;
251  for (k=0; k<ncol; k++)
252  a[ipt+k]/=pivot;
253 
254  ipt=ipiv*nvec;
255  for (k=0; k<nvec; k++)
256  b[ipt+k]/=pivot;
257 
258  /* reduce non-pivot rows such that they will have 0 in the pivot
259  column, which will be discarded, and the augmented matrix will
260  start with 0 from the identity matrix and then have non-zero
261  in the corresponding column, which is part of the inverse.
262  only one column of the augmented matrix is populated at a time,
263  which corresponds to the only column of the original matrix
264  being zeroed, so that the inverse may be done in place. */
265 
266  for (j=0; j<nrow; j++) {
267  if (j == ipiv) continue;
268 
269  dtemp=a[j*ndim+ipiv];
270  a[j*ndim+ipiv]=0.;
271 
272  if (fabs(dtemp) > DBL_EPSILON) {
273  ipt=j *ndim;
274  jpt=ipiv*ndim;
275  for (k=0; k<ncol; k++)
276  a[ipt+k]-=dtemp*a[jpt+k];
277 
278  ipt=j *nvec;
279  jpt=ipiv*nvec;
280  for (k=0; k<nvec; k++)
281  b[ipt+k]-=dtemp*b[jpt+k];
282  }
283  }
284 
285  /* for a diagonal matrix, the determinant is the product of the
286  diagonal terms, and so it may be accumulated from the pivots,
287  noting that switching rows changes the sign as above */
288 
289  det*=pivot;
290  }
291 
292  /* switch columns back in reverse order, noting that a row switch
293  in the original matrix corresponds to a column switch in the
294  inverse matrix */
295 
296  for (i=0; i<nrow; i++) {
297  j=(nrow-1)-i;
298 
299  if (pivrc1[j] != pivrc2[j]) {
300  irow=pivrc1[j];
301  icol=pivrc2[j];
302 
303  // _printf_("column switch back for j=" << j << ": irow=" << irow << ", icol=" << icol << "\n\n");
304 
305  ipt=0;
306  for (k=0; k<nrow; k++) {
307  dtemp =a[ipt+irow];
308  a[ipt+irow]=a[ipt+icol];
309  a[ipt+icol]=dtemp;
310  ipt+=ndim;
311  }
312  }
313  }
314 
315  if (pdet) *pdet=det;
316  xDelete<int>(pivrc1);
317  xDelete<int>(pivrc2);
318  xDelete<int>(pindx);
319  return noerr;
320 }/*}}}*/

◆ Matrix2x2Invert()

void Matrix2x2Invert ( IssmDouble *  Ainv,
IssmDouble *  A 
)

Definition at line 329 of file MatrixUtils.cpp.

329  {/*{{{*/
330 
331  /*Intermediaries*/
332  IssmDouble det,det_reciprocal;
333 
334  /*Compute determinant*/
335  Matrix2x2Determinant(&det,A);
336  if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
337 
338  /*Multiplication is faster than divsion, so we multiply by the reciprocal*/
339  det_reciprocal = 1./det;
340 
341  /*Compute invert*/
342  Ainv[0]= A[3]*det_reciprocal; /* = d/det */
343  Ainv[1]= - A[1]*det_reciprocal; /* = -b/det */
344  Ainv[2]= - A[2]*det_reciprocal; /* = -c/det */
345  Ainv[3]= A[0]*det_reciprocal; /* = a/det */
346 
347 }/*}}}*/

◆ Matrix2x2Determinant()

void Matrix2x2Determinant ( IssmDouble *  Adet,
IssmDouble *  A 
)

Definition at line 322 of file MatrixUtils.cpp.

322  {/*{{{*/
323  /*Compute determinant of a 2x2 matrix*/
324 
325  /*det = a*d - c*b*/
326  *Adet= A[0]*A[3]-A[2]*A[1];
327 }

◆ Matrix2x2Eigen()

void Matrix2x2Eigen ( IssmDouble *  plambda1,
IssmDouble *  plambda2,
IssmDouble *  pvx,
IssmDouble *  pvy,
IssmDouble  a11,
IssmDouble  a21,
IssmDouble  a22 
)

Definition at line 348 of file MatrixUtils.cpp.

348  {/*{{{*/
349  /*From symetric matrix (a11,a21;a21,a22), get eigen values lambda1 and lambda2 and one eigen vector v*/
350 
351  /*Output*/
352  IssmDouble lambda1,lambda2;
353  IssmDouble vx,vy;
354 
355  /*To get the eigen values, we must solve the following equation:
356  * | a11 - lambda a21 |
357  * det | | = 0
358  * | a21 a22-lambda |
359  *
360  * We have to solve the following polynom:
361  * lamda^2 + ( -a11 -a22)*lambda + (a11*a22-a21*a21) = 0*/
362 
363  /*Compute polynom determinant*/
364  IssmDouble b=-a11-a22;
365  IssmDouble delta=b*b - 4*(a11*a22-a21*a21);
366 
367  /*Compute norm of M to avoid round off errors*/
368  IssmDouble normM=a11*a11 + a22*a22 + a21*a21;
369 
370  /*1: normM too small: eigen values = 0*/
371  if(normM<1.e-30){
372  lambda1=0.;
373  lambda2=0.;
374  vx=1.;
375  vy=0.;
376  }
377  /*2: delta is small -> double root*/
378  else if (delta < 1.e-5*normM){
379  lambda1=-b/2.;
380  lambda2=-b/2.;
381  vx=1.;
382  vy=0.;
383  }
384  /*3: general case -> two roots*/
385  else{
386  delta = sqrt(delta);
387  lambda1 = (-b-delta)/2.;
388  lambda2 = (-b+delta)/2.;
389 
390  /*Now, one must find the eigen vectors. For that we use the following property of the inner product
391  * <Ax,y> = <x,tAy>
392  * Here, M'(M-lambda*Id) is symmetrical, which gives:
393  * \forall (x,y)\in R²xR² <M'x,y> = <M'y,x>
394  * And we have the following:
395  * if y\in Ker(M'), \forall x\in R² <M'x,y> = <x,M'y> = 0
396  * We have shown that
397  * Im(M') \perp Ker(M')
398  *
399  * To find the eigen vectors of M, we only have to find two vectors
400  * of the image of M' and take their perpendicular as long as they are
401  * not 0.
402  * To do that, we take the images (1,0) and (0,1):
403  * x1 = (a11 - lambda) x2 = a21
404  * y1 = a21 y2 = (a22-lambda)
405  *
406  * We take the vector that has the larger norm and take its perpendicular.*/
407 
408  IssmDouble norm1 = (a11-lambda1)*(a11-lambda1) + a21*a21;
409  IssmDouble norm2 = a21*a21 + (a22-lambda1)*(a22-lambda1);
410 
411  if(norm2<norm1){
412  norm1=sqrt(norm1);
413  vx = - a21/norm1;
414  vy = (a11-lambda1)/norm1;
415  }
416  else{
417  norm2=sqrt(norm2);
418  vx = - (a22-lambda1)/norm2;
419  vy = a21/norm2;
420  }
421  }
422 
423  /*Assign output*/
424  *plambda1 = lambda1;
425  *plambda2 = lambda2;
426  if(pvx) *pvx = vx;
427  if(pvy) *pvy = vy;
428 
429 }/*}}}*/

◆ Matrix3x3Invert()

void Matrix3x3Invert ( IssmDouble *  Ainv,
IssmDouble *  A 
)

Definition at line 448 of file MatrixUtils.cpp.

448  {/*{{{*/
449 
450  /*Intermediaries*/
451  IssmDouble det,det_reciprocal;
452 
453  /*Compute determinant*/
454  Matrix3x3Determinant(&det,A);
455  if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
456 
457  /*Multiplication is faster than divsion, so we multiply by the reciprocal*/
458  det_reciprocal = 1./det;
459 
460  /*Compute invert*/
461  Ainv[0]=(A[4]*A[8]-A[5]*A[7])*det_reciprocal; /* = (e*i-f*h)/det */
462  Ainv[1]=(A[2]*A[7]-A[1]*A[8])*det_reciprocal; /* = (c*h-b*i)/det */
463  Ainv[2]=(A[1]*A[5]-A[2]*A[4])*det_reciprocal; /* = (b*f-c*e)/det */
464  Ainv[3]=(A[5]*A[6]-A[3]*A[8])*det_reciprocal; /* = (f*g-d*i)/det */
465  Ainv[4]=(A[0]*A[8]-A[2]*A[6])*det_reciprocal; /* = (a*i-c*g)/det */
466  Ainv[5]=(A[2]*A[3]-A[0]*A[5])*det_reciprocal; /* = (c*d-a*f)/det */
467  Ainv[6]=(A[3]*A[7]-A[4]*A[6])*det_reciprocal; /* = (d*h-e*g)/det */
468  Ainv[7]=(A[1]*A[6]-A[0]*A[7])*det_reciprocal; /* = (b*g-a*h)/det */
469  Ainv[8]=(A[0]*A[4]-A[1]*A[3])*det_reciprocal; /* = (a*e-b*d)/det */
470 }/*}}}*/

◆ Matrix3x3Determinant() [1/2]

void Matrix3x3Determinant ( IssmDouble *  Adet,
IssmDouble *  A 
)

Definition at line 431 of file MatrixUtils.cpp.

431  {/*{{{*/
432  /*Compute determinant of a 3x3 matrix*/
433 
434  /*det = a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)*/
435  *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];
436 }

◆ Matrix3x3Determinant() [2/2]

IssmDouble Matrix3x3Determinant ( IssmDouble  a1,
IssmDouble  a2,
IssmDouble  a3,
IssmDouble  b1,
IssmDouble  b2,
IssmDouble  b3,
IssmDouble  c1,
IssmDouble  c2,
IssmDouble  c3 
)

Definition at line 438 of file MatrixUtils.cpp.

438  {/*{{{*/
439  /*Compute determinant of a 3x3 matrix*/
440 
441  /*det = a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)
442  * a b c a1 a2 a3
443  * d e f b1 b2 b3
444  * g h i c1 c2 c3 */
445  return a1*b2*c3-a1*b3*c2-b1*a2*c3+b1*a3*c2+c1*a2*b3-c1*a3*b2;
446 }

◆ Matrix3x3Solve()

void Matrix3x3Solve ( IssmDouble *  X,
IssmDouble *  A,
IssmDouble *  B 
)

Definition at line 471 of file MatrixUtils.cpp.

471  {/*{{{*/
472 
473  IssmDouble Ainv[3][3];
474 
475  Matrix3x3Invert(&Ainv[0][0],A);
476  for(int i=0;i<3;i++) X[i]=Ainv[i][0]*B[0] + Ainv[i][1]*B[1] + Ainv[i][2]*B[2];
477 
478 }/*}}}*/

◆ Matrix4x4Adjoint()

void Matrix4x4Adjoint ( IssmDouble *  Aadj,
IssmDouble *  A 
)

Definition at line 509 of file MatrixUtils.cpp.

509  {/*{{{*/
510 
511  IssmDouble a1 = A[0*4+0];
512  IssmDouble b1 = A[0*4+1];
513  IssmDouble c1 = A[0*4+2];
514  IssmDouble d1 = A[0*4+3];
515 
516  IssmDouble a2 = A[1*4+0];
517  IssmDouble b2 = A[1*4+1];
518  IssmDouble c2 = A[1*4+2];
519  IssmDouble d2 = A[1*4+3];
520 
521  IssmDouble a3 = A[2*4+0];
522  IssmDouble b3 = A[2*4+1];
523  IssmDouble c3 = A[2*4+2];
524  IssmDouble d3 = A[2*4+3];
525 
526  IssmDouble a4 = A[3*4+0];
527  IssmDouble b4 = A[3*4+1];
528  IssmDouble c4 = A[3*4+2];
529  IssmDouble d4 = A[3*4+3];
530 
531  /* Row column labeling reversed since we transpose rows & columns*/
532  Aadj[0*4+0] = Matrix3x3Determinant(b2, b3, b4, c2, c3, c4, d2, d3, d4);
533  Aadj[1*4+0] = - Matrix3x3Determinant(a2, a3, a4, c2, c3, c4, d2, d3, d4);
534  Aadj[2*4+0] = Matrix3x3Determinant(a2, a3, a4, b2, b3, b4, d2, d3, d4);
535  Aadj[3*4+0] = - Matrix3x3Determinant(a2, a3, a4, b2, b3, b4, c2, c3, c4);
536 
537  Aadj[0*4+1] = - Matrix3x3Determinant(b1, b3, b4, c1, c3, c4, d1, d3, d4);
538  Aadj[1*4+1] = Matrix3x3Determinant(a1, a3, a4, c1, c3, c4, d1, d3, d4);
539  Aadj[2*4+1] = - Matrix3x3Determinant(a1, a3, a4, b1, b3, b4, d1, d3, d4);
540  Aadj[3*4+1] = Matrix3x3Determinant(a1, a3, a4, b1, b3, b4, c1, c3, c4);
541 
542  Aadj[0*4+2] = Matrix3x3Determinant(b1, b2, b4, c1, c2, c4, d1, d2, d4);
543  Aadj[1*4+2] = - Matrix3x3Determinant(a1, a2, a4, c1, c2, c4, d1, d2, d4);
544  Aadj[2*4+2] = Matrix3x3Determinant(a1, a2, a4, b1, b2, b4, d1, d2, d4);
545  Aadj[3*4+2] = - Matrix3x3Determinant(a1, a2, a4, b1, b2, b4, c1, c2, c4);
546 
547  Aadj[0*4+3] = - Matrix3x3Determinant(b1, b2, b3, c1, c2, c3, d1, d2, d3);
548  Aadj[1*4+3] = Matrix3x3Determinant(a1, a2, a3, c1, c2, c3, d1, d2, d3);
549  Aadj[2*4+3] = - Matrix3x3Determinant(a1, a2, a3, b1, b2, b3, d1, d2, d3);
550  Aadj[3*4+3] = Matrix3x3Determinant(a1, a2, a3, b1, b2, b3, c1, c2, c3);
551 }/*}}}*/

◆ Matrix4x4Invert()

void Matrix4x4Invert ( IssmDouble *  Ainv,
IssmDouble *  A 
)

Definition at line 552 of file MatrixUtils.cpp.

552  {/*{{{*/
553 
554  /*Intermediaries*/
555  IssmDouble det,det_reciprocal;
556 
557  /*Compute determinant*/
558  Matrix4x4Determinant(&det,A);
559  if(fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
560 
561  /*Multiplication is faster than division, so we multiply by the reciprocal*/
562  det_reciprocal = 1./det;
563 
564  /*Compute adjoint matrix*/
565  Matrix4x4Adjoint(Ainv,A);
566 
567  /*Scalte adjoint matrix to get inverse*/
568  for(int i=0;i<4*4;i++) Ainv[i] = Ainv[i]*det_reciprocal;
569 }/*}}}*/

◆ Matrix4x4Determinant()

void Matrix4x4Determinant ( IssmDouble *  Adet,
IssmDouble *  A 
)

Definition at line 480 of file MatrixUtils.cpp.

480  {/*{{{*/
481  /*Compute determinant of a 4x4 matrix*/
482 
483  IssmDouble a1 = A[0*4+0];
484  IssmDouble b1 = A[0*4+1];
485  IssmDouble c1 = A[0*4+2];
486  IssmDouble d1 = A[0*4+3];
487 
488  IssmDouble a2 = A[1*4+0];
489  IssmDouble b2 = A[1*4+1];
490  IssmDouble c2 = A[1*4+2];
491  IssmDouble d2 = A[1*4+3];
492 
493  IssmDouble a3 = A[2*4+0];
494  IssmDouble b3 = A[2*4+1];
495  IssmDouble c3 = A[2*4+2];
496  IssmDouble d3 = A[2*4+3];
497 
498  IssmDouble a4 = A[3*4+0];
499  IssmDouble b4 = A[3*4+1];
500  IssmDouble c4 = A[3*4+2];
501  IssmDouble d4 = A[3*4+3];
502 
503  *Adet= a1 * Matrix3x3Determinant(b2, b3, b4, c2, c3, c4, d2, d3, d4)
504  - b1 * Matrix3x3Determinant(a2, a3, a4, c2, c3, c4, d2, d3, d4)
505  + c1 * Matrix3x3Determinant(a2, a3, a4, b2, b3, b4, d2, d3, d4)
506  - d1 * Matrix3x3Determinant(a2, a3, a4, b2, b3, b4, c2, c3, c4);
507 }

◆ Matrix4x4Solve()

void Matrix4x4Solve ( IssmDouble *  X,
IssmDouble *  A,
IssmDouble *  B 
)

Definition at line 570 of file MatrixUtils.cpp.

570  {/*{{{*/
571  IssmDouble Ainv[4][4];
572 
573  Matrix4x4Invert(&Ainv[0][0],A);
574  for(int i=0;i<4;i++) X[i]=Ainv[i][0]*B[0] + Ainv[i][1]*B[1] + Ainv[i][2]*B[2] + Ainv[i][3]*B[3];
575 }/*}}}*/

◆ newcell()

void newcell ( IssmDouble **  pcell,
IssmDouble  newvalue,
bool  top,
int  m 
)

Definition at line 577 of file MatrixUtils.cpp.

577  { /*{{{*/
578 
579  IssmDouble* cell=NULL;
580  IssmDouble* dummy=NULL;
581 
582  /*recover pointer: */
583  cell=*pcell;
584 
585  /*reallocate:*/
586  dummy=xNew<IssmDouble>(m+1);
587 
588  /*copy data:*/
589  if(top){
590  dummy[0]=newvalue;
591  for(int i=0;i<m;i++)dummy[i+1]=cell[i];
592  }
593  else{
594  dummy[m]=newvalue;
595  for(int i=0;i<m;i++)dummy[i]=cell[i];
596  }
597 
598  /*delete and reassign: */
599  xDelete<IssmDouble>(cell); cell=dummy;
600 
601  /*assign output pointer:*/
602  *pcell=cell;
603 } /*}}}*/

◆ cellsum()

IssmDouble cellsum ( IssmDouble *  cell,
int  m 
)

Definition at line 604 of file MatrixUtils.cpp.

604  { /*{{{*/
605 
606  IssmDouble sum=0;
607 
608  for(int i=0;i<m;i++)sum+=cell[i];
609 
610  return sum;
611 } /*}}}*/

◆ celldelete()

void celldelete ( IssmDouble **  pcell,
int  m,
int *  indices,
int  nind 
)

Definition at line 612 of file MatrixUtils.cpp.

612  { /*{{{*/
613 
614  /*input: */
615  IssmDouble* cell=*pcell;
616 
617  /*output: */
618  IssmDouble* newcell=xNew<IssmDouble>(nind);
619 
620  for(int i=0;i<nind;i++){
621  newcell[i]=cell[indices[i]];
622  }
623 
624  /*free allocation:*/
625  xDelete<IssmDouble>(cell);
626 
627  /*assign output pointers: */
628  *pcell=newcell;
629 } /*}}}*/

◆ cellsplit()

void cellsplit ( IssmDouble **  pcell,
int  m,
int  i,
IssmDouble  scale 
)

Definition at line 630 of file MatrixUtils.cpp.

630  { /*{{{*/
631 
632  /*input: */
633  IssmDouble* cell=*pcell;
634 
635  /*output: */
636  IssmDouble* newcell=xNew<IssmDouble>(m+1);
637 
638  for(int j=0;j<i;j++)newcell[j]=cell[j];
639  newcell[i]=scale*cell[i];
640  newcell[i+1]=scale* cell[i];
641  for(int j=i+2;j<m+1;j++)newcell[j]=cell[j-1];
642 
643  /*free allocation:*/
644  xDelete<IssmDouble>(cell);
645 
646  /*assign output pointers: */
647  *pcell=newcell;
648 } /*}}}*/

◆ cellecho()

void cellecho ( int  numcells,
int  m,
  ... 
)

Definition at line 649 of file MatrixUtils.cpp.

649  { /*{{{*/
650 
651  va_list arguments;
652  IssmDouble** celllist= NULL;
653 
654  /*allocate variable length array: */
655  celllist=xNew<IssmDouble*>(numcells);
656 
657  va_start(arguments,m);
658 
659  for ( int x = 0; x < numcells; x++ ){
660  celllist[x]= va_arg ( arguments, IssmDouble*);
661  }
662  va_end ( arguments );
663 
664  _printf_("Echo of cell: \n");
665  for(int i=0;i<m;i++){
666  _printf_(i << ": ");
667  for (int j=0;j<numcells;j++)_printf_(setprecision(10) << celllist[j][i] << " ");
668  _printf_("\n");
669  }
670 
671  /*deallocate:*/
672  xDelete<IssmDouble*>(celllist);
673 
674 } /*}}}*/
Matrix4x4Adjoint
void Matrix4x4Adjoint(IssmDouble *Aadj, IssmDouble *A)
Definition: MatrixUtils.cpp:509
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
Matrix3x3Invert
void Matrix3x3Invert(IssmDouble *Ainv, IssmDouble *A)
Definition: MatrixUtils.cpp:448
Matrix4x4Invert
void Matrix4x4Invert(IssmDouble *Ainv, IssmDouble *A)
Definition: MatrixUtils.cpp:552
bamg::det
long long det(const I2 &a, const I2 &b, const I2 &c)
Definition: det.h:8
newcell
void newcell(IssmDouble **pcell, IssmDouble newvalue, bool top, int m)
Definition: MatrixUtils.cpp:577
Matrix4x4Determinant
void Matrix4x4Determinant(IssmDouble *Adet, IssmDouble *A)
Definition: MatrixUtils.cpp:480
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
MatrixMultiply
int MatrixMultiply(IssmDouble *a, int nrowa, int ncola, int itrna, IssmDouble *b, int nrowb, int ncolb, int itrnb, IssmDouble *c, int iaddc)
Definition: MatrixUtils.cpp:88
Matrix3x3Determinant
void Matrix3x3Determinant(IssmDouble *Adet, IssmDouble *A)
Definition: MatrixUtils.cpp:431
Matrix2x2Determinant
void Matrix2x2Determinant(IssmDouble *Adet, IssmDouble *A)
Definition: MatrixUtils.cpp:322
Generated on Thu Jul 2 2020 08:09:25 for Ice Sheet System Model by   doxygen 1.8.19