Index: /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 20727)
+++ /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 20728)
@@ -22,5 +22,5 @@
 
 	int         idima,idimb,idimc,idimd;
-	IssmDouble *dtemp;
+	IssmDouble dtemp_static[600];
 	_assert_(a && b && c && d);
 
@@ -52,5 +52,5 @@
 		idimd=nrowc;
 	}
-
+	
 	/*  perform the matrix triple product in the order that minimizes the
 		 number of multiplies and the temporary space used, noting that
@@ -61,18 +61,32 @@
 	/*  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,14 +338,18 @@
 	/*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,4 +459,5 @@
 	/*Intermediaries*/
 	IssmDouble det;
+	IssmDouble det_reciprocal;
 
 	/*Compute determinant*/
@@ -446,15 +465,17 @@
 	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,4 +564,5 @@
 	/*Intermediaries*/
 	IssmDouble det;
+	IssmDouble det_reciprocal;
 
 	/*Compute determinant*/
@@ -548,9 +570,12 @@
 	//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){/*{{{*/