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);