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

Index: ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp    (revision 20728)
+++ ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp    (revision 20729)
@@ -21,7 +21,9 @@
        /*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,6 +53,15 @@
                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,36 +71,18 @@
 
        /*  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,15 +329,14 @@
 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,15 +449,14 @@
 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,15 +553,14 @@
 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);