source: issm/oecreview/Archive/20545-21336/ISSM-20728-20729.diff@ 21337

../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp

@@ -21 +21 @@
         /*TripleMultiply    Perform triple matrix product a*b*c+d.*/
 
         int         idima,idimb,idimc,idimd;
-        IssmDouble dtemp_static[600];
+        IssmDouble  dtemp_static[600];
+        IssmDouble* dtemp_dynamic = NULL;
+        IssmDouble* dtemp         = NULL;
         _assert_(a && b && c && d);
 
         /*  set up dimensions for triple product  */
@@ -51 +53 @@
                 if (ncolc != idimc) _error_("Matrix B and C inner vectors not equal size.");
                 idimd=nrowc;
         }
+
+        /*Depending on the size of incoming matrices, we might need to use a dynamic allocation*/
+        if(idima*idimc>600){
+                IssmDouble *dtemp_dynamic=xNew<IssmDouble>(idima*idimc);
+                dtemp = dtemp_dynamic;
+        }
+        else{
+                dtemp = &dtemp_static[0];
+        }
         
         /*  perform the matrix triple product in the order that minimizes the
                  number of multiplies and the temporary space used, noting that
@@ -60 +71 @@
 
         /*  multiply (a*b)*c+d  */
         if (idima*idimc*(idimb+idimd) <= idimb*idimd*(idima+idimc)) {
-                if(idima*idimc>600) {
-                        /*Use dynamic allocation instead of dtemp_static*/
-                        IssmDouble *dtemp;
-                        dtemp=xNew<IssmDouble>(idima*idimc);
-
-                        MatrixMultiply(a,nrowa,ncola,itrna,b,nrowb,ncolb,itrnb,dtemp,0);
-                        MatrixMultiply(dtemp,idima,idimc,0,c,nrowc,ncolc,itrnc,d,iaddd);
-                        xDelete<IssmDouble>(dtemp);
-                        return 1;
-                }
-                MatrixMultiply(a,nrowa,ncola,itrna,b,nrowb,ncolb,itrnb,&dtemp_static[0],0);
-                MatrixMultiply(&dtemp_static[0],idima,idimc,0,c,nrowc,ncolc,itrnc,d,iaddd);
+                MatrixMultiply(a,nrowa,ncola,itrna,b,nrowb,ncolb,itrnb,dtemp,0);
+                MatrixMultiply(dtemp,idima,idimc,0,c,nrowc,ncolc,itrnc,d,iaddd);
         }
 
         /*  multiply a*(b*c)+d  */
         else{
-                if(idimb*idimd>600) {
-                        /*Use dynamic allocation instead of dtemp_static*/
-                        IssmDouble *dtemp;
-                        dtemp=xNew<IssmDouble>(idimb*idimd);
-
-                        MatrixMultiply(b,nrowb,ncolb,itrnb,c,nrowc,ncolc,itrnc,dtemp,0);
-                        MatrixMultiply(a,nrowa,ncola,itrna,dtemp,idimb,idimd,0,d,iaddd);
-                        xDelete<IssmDouble>(dtemp);
-                        return 1;
-                }
-                MatrixMultiply(b,nrowb,ncolb,itrnb,c,nrowc,ncolc,itrnc,&dtemp_static[0],0);
-                MatrixMultiply(a,nrowa,ncola,itrna,&dtemp_static[0],idimb,idimd,0,d,iaddd);
+                MatrixMultiply(b,nrowb,ncolb,itrnb,c,nrowc,ncolc,itrnc,dtemp,0);
+                MatrixMultiply(a,nrowa,ncola,itrna,dtemp,idimb,idimd,0,d,iaddd);
         }
 
+        /*Cleanup and return*/
+        xDelete<IssmDouble>(dtemp);
         return 1;
 }/*}}}*/
 int MatrixMultiply(IssmDouble* a, int nrowa, int ncola, int itrna, IssmDouble* b, int nrowb, int ncolb, int itrnb, IssmDouble* c, int iaddc ){/*{{{*/
@@ -336 +329 @@
 void Matrix2x2Invert(IssmDouble* Ainv,IssmDouble* A){/*{{{*/
 
         /*Intermediaries*/
-        IssmDouble det;
-        IssmDouble det_reciprocal;
+        IssmDouble det,det_reciprocal;
 
         /*Compute determinant*/
         Matrix2x2Determinant(&det,A);
         if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
         
         /*Multiplication is faster than divsion, so we multiply by the reciprocal*/
-        det_reciprocal = 1/det;  
+        det_reciprocal = 1./det;  
 
         /*Compute invert*/
         Ainv[0]=   A[3]*det_reciprocal; /* =  d/det */
@@ -457 +449 @@
 void Matrix3x3Invert(IssmDouble* Ainv,IssmDouble* A){/*{{{*/
 
         /*Intermediaries*/
-        IssmDouble det;
-        IssmDouble det_reciprocal;
+        IssmDouble det,det_reciprocal;
 
         /*Compute determinant*/
         Matrix3x3Determinant(&det,A);
         if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
 
         /*Multiplication is faster than divsion, so we multiply by the reciprocal*/
-        det_reciprocal = 1/det;  
+        det_reciprocal = 1./det;  
 
         /*Compute invert*/
         Ainv[0]=(A[4]*A[8]-A[5]*A[7])*det_reciprocal; /* = (e*i-f*h)/det */
@@ -562 +553 @@
 void Matrix4x4Invert(IssmDouble* Ainv,IssmDouble* A){/*{{{*/
 
         /*Intermediaries*/
-        IssmDouble det;
-        IssmDouble det_reciprocal;
+        IssmDouble det,det_reciprocal;
 
         /*Compute determinant*/
         Matrix4x4Determinant(&det,A);
-        //if(fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
+        if(fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
 
         /*Multiplication is faster than division, so we multiply by the reciprocal*/
-        det_reciprocal = 1/det;
+        det_reciprocal = 1./det;
 
         /*Compute adjoint matrix*/
         Matrix4x4Adjoint(Ainv,A);