Changeset 20728

Timestamp:

06/14/16 00:11:40 (9 years ago)

Author:

agscott1

Message:

CHG TripleMultiply now uses static allocation in most cases, and MatrixInvert functions now multiply by the reciprocal of a determinant instead of dividing by the determinant.

File:

• issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (modified) (8 diffs)

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

@@ -22 +22 @@
 
         int         idima,idimb,idimc,idimd;
-        IssmDouble *dtemp;
+        IssmDouble dtemp_static[600];
         _assert_(a && b && c && d);
 
@@ -52 +52 @@
                 idimd=nrowc;
         }
-
+        
         /*  perform the matrix triple product in the order that minimizes the
                  number of multiplies and the temporary space used, noting that
@@ -61 +61 @@
         /*  multiply (a*b)*c+d  */
         if (idima*idimc*(idimb+idimd) <= idimb*idimd*(idima+idimc)) {
-                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);
+                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);
         }
 
         /*  multiply a*(b*c)+d  */
         else{
-                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);
+                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);
         }
 
@@ -324 +338 @@
         /*Intermediaries*/
         IssmDouble det;
+        IssmDouble 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;  
 
         /*Compute invert*/
-        Ainv[0]=   A[3]/det; /* =  d/det */
-        Ainv[1]= - A[1]/det; /* = -b/det */
-        Ainv[2]= - A[2]/det; /* = -c/det */
-        Ainv[3]=   A[0]/det; /* =  a/det */
+        Ainv[0]=   A[3]*det_reciprocal; /* =  d/det */
+        Ainv[1]= - A[1]*det_reciprocal; /* = -b/det */
+        Ainv[2]= - A[2]*det_reciprocal; /* = -c/det */
+        Ainv[3]=   A[0]*det_reciprocal; /* =  a/det */
 
 }/*}}}*/
@@ -441 +459 @@
         /*Intermediaries*/
         IssmDouble det;
+        IssmDouble det_reciprocal;
 
         /*Compute determinant*/
@@ -446 +465 @@
         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;  
+
         /*Compute invert*/
-        Ainv[0]=(A[4]*A[8]-A[5]*A[7])/det; /* = (e*i-f*h)/det */
-        Ainv[1]=(A[2]*A[7]-A[1]*A[8])/det; /* = (c*h-b*i)/det */
-        Ainv[2]=(A[1]*A[5]-A[2]*A[4])/det; /* = (b*f-c*e)/det */
-        Ainv[3]=(A[5]*A[6]-A[3]*A[8])/det; /* = (f*g-d*i)/det */
-        Ainv[4]=(A[0]*A[8]-A[2]*A[6])/det; /* = (a*i-c*g)/det */
-        Ainv[5]=(A[2]*A[3]-A[0]*A[5])/det; /* = (c*d-a*f)/det */
-        Ainv[6]=(A[3]*A[7]-A[4]*A[6])/det; /* = (d*h-e*g)/det */
-        Ainv[7]=(A[1]*A[6]-A[0]*A[7])/det; /* = (b*g-a*h)/det */
-        Ainv[8]=(A[0]*A[4]-A[1]*A[3])/det; /* = (a*e-b*d)/det */
-
+        Ainv[0]=(A[4]*A[8]-A[5]*A[7])*det_reciprocal; /* = (e*i-f*h)/det */
+        Ainv[1]=(A[2]*A[7]-A[1]*A[8])*det_reciprocal; /* = (c*h-b*i)/det */
+        Ainv[2]=(A[1]*A[5]-A[2]*A[4])*det_reciprocal; /* = (b*f-c*e)/det */
+        Ainv[3]=(A[5]*A[6]-A[3]*A[8])*det_reciprocal; /* = (f*g-d*i)/det */
+        Ainv[4]=(A[0]*A[8]-A[2]*A[6])*det_reciprocal; /* = (a*i-c*g)/det */
+        Ainv[5]=(A[2]*A[3]-A[0]*A[5])*det_reciprocal; /* = (c*d-a*f)/det */
+        Ainv[6]=(A[3]*A[7]-A[4]*A[6])*det_reciprocal; /* = (d*h-e*g)/det */
+        Ainv[7]=(A[1]*A[6]-A[0]*A[7])*det_reciprocal; /* = (b*g-a*h)/det */
+        Ainv[8]=(A[0]*A[4]-A[1]*A[3])*det_reciprocal; /* = (a*e-b*d)/det */
 }/*}}}*/
 void Matrix3x3Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B){/*{{{*/
@@ -543 +564 @@
         /*Intermediaries*/
         IssmDouble det;
+        IssmDouble det_reciprocal;
 
         /*Compute determinant*/
@@ -548 +570 @@
         //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;
+
         /*Compute adjoint matrix*/
         Matrix4x4Adjoint(Ainv,A);
 
         /*Scalte adjoint matrix to get inverse*/
-        for(int i=0;i<4*4;i++) Ainv[i] = Ainv[i]/det;
+        for(int i=0;i<4*4;i++) Ainv[i] = Ainv[i]*det_reciprocal;
 }/*}}}*/
 void Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble *B){/*{{{*/