Index: /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 20728)
+++ /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 20729)
@@ -22,5 +22,7 @@
 
 	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);
 
@@ -51,4 +53,13 @@
 		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];
 	}
 	
@@ -61,34 +72,16 @@
 	/*  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;
 }/*}}}*/
@@ -337,6 +330,5 @@
 
 	/*Intermediaries*/
-	IssmDouble det;
-	IssmDouble det_reciprocal;
+	IssmDouble det,det_reciprocal;
 
 	/*Compute determinant*/
@@ -345,5 +337,5 @@
 	
 	/*Multiplication is faster than divsion, so we multiply by the reciprocal*/
-	det_reciprocal = 1/det;  
+	det_reciprocal = 1./det;  
 
 	/*Compute invert*/
@@ -458,6 +450,5 @@
 
 	/*Intermediaries*/
-	IssmDouble det;
-	IssmDouble det_reciprocal;
+	IssmDouble det,det_reciprocal;
 
 	/*Compute determinant*/
@@ -466,5 +457,5 @@
 
 	/*Multiplication is faster than divsion, so we multiply by the reciprocal*/
-	det_reciprocal = 1/det;  
+	det_reciprocal = 1./det;  
 
 	/*Compute invert*/
@@ -563,13 +554,12 @@
 
 	/*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*/