Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
===================================================================
--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 13141)
+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 13142)
@@ -16,57 +16,48 @@
 /*}}}*/
 
 /*FUNCTION TripleMultiply {{{*/
-int TripleMultiply( IssmDouble* a, int nrowa, int ncola, int itrna, IssmDouble* b, int nrowb, int ncolb, int itrnb, IssmDouble* c, int nrowc, int ncolc, int itrnc, IssmDouble* d, int iaddd){
+int TripleMultiply(IssmDouble* a, int nrowa, int ncola, int itrna, IssmDouble* b, int nrowb, int ncolb, int itrnb, IssmDouble* c, int nrowc, int ncolc, int itrnc, IssmDouble* d, int iaddd){
 	/*TripleMultiply    Perform triple matrix product a*b*c+d.*/
-	
-	int idima,idimb,idimc,idimd;
-	IssmDouble* dtemp;
 
-/*  set up dimensions for triple product  */
+	int         idima,idimb,idimc,idimd;
+	IssmDouble *dtemp;
 
-	if (!itrna) {
+	/*  set up dimensions for triple product  */
+
+	if (!itrna){
 		idima=nrowa;
 		idimb=ncola;
 	}
-	else {
+	else{
 		idima=ncola;
 		idimb=nrowa;
 	}
 
-	if (!itrnb) {
-		if (nrowb != idimb) {
-			_error_("Matrix A and B inner vectors not equal size.");
-		}
+	if (!itrnb){
+		if (nrowb != idimb) _error_("Matrix A and B inner vectors not equal size.");
 		idimc=ncolb;
 	}
-	else {
-		if (ncolb != idimb) {
-			_error_("Matrix A and B inner vectors not equal size.");
-		}
+	else{
+		if (ncolb != idimb) _error_("Matrix A and B inner vectors not equal size.");
 		idimc=nrowb;
 	}
 
 	if (!itrnc) {
-		if (nrowc != idimc) {
-			_error_("Matrix B and C inner vectors not equal size.");
-		}
+		if (nrowc != idimc) _error_("Matrix B and C inner vectors not equal size.");
 		idimd=ncolc;
 	}
-	else {
-		if (ncolc != idimc) {
-			_error_("Matrix B and C inner vectors not equal size.");
-		}
+	else{
+		if (ncolc != idimc) _error_("Matrix B and C inner vectors not equal size.");
 		idimd=nrowc;
 	}
 
-/*  perform the matrix triple product in the order that minimizes the
-	number of multiplies and the temporary space used, noting that
-	(a*b)*c requires ac(b+d) multiplies and ac IssmDoubles, and a*(b*c)
-	requires bd(a+c) multiplies and bd IssmDoubles (both are the same for
-	a symmetric triple product)  */
+	/*  perform the matrix triple product in the order that minimizes the
+		 number of multiplies and the temporary space used, noting that
+		 (a*b)*c requires ac(b+d) multiplies and ac IssmDoubles, and a*(b*c)
+		 requires bd(a+c) multiplies and bd IssmDoubles (both are the same for
+		 a symmetric triple product)  */
 
-/*  multiply (a*b)*c+d  */
-
+	/*  multiply (a*b)*c+d  */
 	if (idima*idimc*(idimb+idimd) <= idimb*idimd*(idima+idimc)) {
 		dtemp=xNew<IssmDouble>(idima*idimc);
 
@@ -75,9 +66,8 @@
 		xDelete<IssmDouble>(dtemp);
 	}
 
-/*  multiply a*(b*c)+d  */
-
-	else {
+	/*  multiply a*(b*c)+d  */
+	else{
 		dtemp=xNew<IssmDouble>(idimb*idimd);
 
 		MatrixMultiply(b,nrowb,ncolb,itrnb,c,nrowc,ncolc,itrnc,dtemp,0);
@@ -88,13 +78,13 @@
 	return 1;
 }/*}}}*/
 /*FUNCTION MatrixMuliply {{{*/
-int MatrixMultiply( IssmDouble* a, int nrowa, int ncola, int itrna, IssmDouble* b, int nrowb, int ncolb, int itrnb, IssmDouble* c, int iaddc ){
+int MatrixMultiply(IssmDouble* a, int nrowa, int ncola, int itrna, IssmDouble* b, int nrowb, int ncolb, int itrnb, IssmDouble* c, int iaddc ){
 	/*MatrixMultiply    Perform matrix multiplication a*b+c.*/
 	int noerr=1;
 	int i,j,k,ipta,iptb,iptc;
 	int nrowc,ncolc,iinca,jinca,iincb,jincb,ntrma,ntrmb,nterm;
 
-/*  set up dimensions and increments for matrix a  */
+	/*  set up dimensions and increments for matrix a  */
 	if (!itrna) {
 		nrowc=nrowa;
 		ntrma=ncola;
@@ -108,7 +98,7 @@
 		jinca=ncola;
 	}
 
-/*  set up dimensions and increments for matrix b  */
+	/*  set up dimensions and increments for matrix b  */
 	if (!itrnb) {
 		ncolc=ncolb;
 		ntrmb=nrowb;
@@ -122,34 +112,25 @@
 		jincb=ncolb;
 	}
 
-	if (ntrma != ntrmb) {
-		_error_("Matrix A and B inner vectors not equal size");
-	    noerr=0;	
-		return noerr;
-	}
-	else
-		nterm=ntrma;
+	if (ntrma != ntrmb) _error_("Matrix A and B inner vectors not equal size");
 
-/*  zero matrix c, if not being added to product  */
+	nterm=ntrma;
 
-	if (!iaddc)
-		for (i=0; i<nrowc*ncolc; i++)
-			*(c+i)=0.;
+	/*  zero matrix c, if not being added to product  */
+	if (!iaddc) for (i=0;i<nrowc*ncolc;i++) c[i]=0.;
 
-/*  perform the matrix multiplication  */
-
+	/*  perform the matrix multiplication  */
 	iptc=0;
-	for (i=0; i<nrowc; i++) {
-		for (j=0; j<ncolc; j++) {
+	for (i=0; i<nrowc; i++){
+		for (j=0; j<ncolc; j++){
 			ipta=i*iinca;
 			iptb=j*jincb;
 
-			for (k=0; k<nterm; k++) {
-				*(c+iptc)+=*(a+ipta)**(b+iptb);
+			for (k=0; k<nterm; k++){
+				c[iptc]+=a[ipta]*b[iptb];
 				ipta+=jinca;
 				iptb+=iincb;
 			}
-
 			iptc++;
 		}
 	}
@@ -158,16 +139,16 @@
 }/*}}}*/
 /*FUNCTION MatrixInverse {{{*/
 int MatrixInverse( IssmDouble* a, int ndim, int nrow, IssmDouble* b, int nvec, IssmDouble* pdet ){
-/* MatrixInverse    Perform matrix inversion and linear equation solution.
+	/* MatrixInverse    Perform matrix inversion and linear equation solution.
 
-	This function uses Gaussian elimination on the original matrix
-	augmented by an identity matrix of the same size to calculate
-	the inverse (see for example, "Modern Methods of Engineering
-	Computation", Sec. 6.4).  By noting how the matrices are
-	unpopulated and repopulated, the calculation may be done in place.
+		This function uses Gaussian elimination on the original matrix
+		augmented by an identity matrix of the same size to calculate
+		the inverse (see for example, "Modern Methods of Engineering
+		Computation", Sec. 6.4).  By noting how the matrices are
+		unpopulated and repopulated, the calculation may be done in place.
 
-	Gaussian elimination is inherently inefficient, and so this is
-	intended for small matrices.  */
+		Gaussian elimination is inherently inefficient, and so this is
+		intended for small matrices.  */
 	int noerr=1;
 	int i,j,k,ipt,jpt,irow,icol,ipiv,ncol;
 	int *pivrc1,*pivrc2,*pindx;
@@ -183,7 +164,7 @@
 	_assert_(!(ndim==2 && nrow==2));
 	_assert_(!(ndim==3 && nrow==3));
 
-/*  initialize local variables and arrays  */
+	/*  initialize local variables and arrays  */
 
 	ncol=nrow;
 	det=1.;
@@ -191,23 +172,23 @@
 	pivrc2 =xNew<int>(nrow);
 	pindx =xNew<int>(nrow);
 
-/*  loop over the rows/columns of the matrix  */
+	/*  loop over the rows/columns of the matrix  */
 
 	for (i=0; i<nrow; i++) {
 
-/*  search for pivot, finding the term with the greatest magnitude
-	in the rows/columns not yet used  */
+		/*  search for pivot, finding the term with the greatest magnitude
+			 in the rows/columns not yet used  */
 
 		pivot=0.;
 		for (j=0; j<nrow; j++)
-			if (!pindx[j])
-				for (k=0; k<ncol; k++)
-					if (!pindx[k])
-						if (fabs(a[j*ndim+k]) > fabs(pivot)) {
-							irow=j;
-							icol=k;
-							pivot=a[j*ndim+k];
-						}
+		 if (!pindx[j])
+		  for (k=0; k<ncol; k++)
+			if (!pindx[k])
+			 if (fabs(a[j*ndim+k]) > fabs(pivot)) {
+				 irow=j;
+				 icol=k;
+				 pivot=a[j*ndim+k];
+			 }
 
 		if (fabs(pivot) < DBL_EPSILON) {
 			xDelete<int>(pivrc1);
@@ -224,15 +205,15 @@
 		ipiv=icol;
 		pindx[ipiv]++;
 
-//		_printf_(true,"pivot for i=%d: irow=%d, icol=%d, pindx[%d]=%d\n",
-//				 i,irow,icol,ipiv,pindx[ipiv]);
+		//		_printf_(true,"pivot for i=%d: irow=%d, icol=%d, pindx[%d]=%d\n",
+		//				 i,irow,icol,ipiv,pindx[ipiv]);
 
-/*  switch rows to put pivot element on diagonal, noting that the
-	column stays the same and the determinant changes sign  */
+		/*  switch rows to put pivot element on diagonal, noting that the
+			 column stays the same and the determinant changes sign  */
 
 		if (irow != icol) {
-//			_printf_(true,"row switch for i=%d: irow=%d, icol=%d\n",
-//					 i,irow,icol);
+			//			_printf_(true,"row switch for i=%d: irow=%d, icol=%d\n",
+			//					 i,irow,icol);
 
 			ipt=irow*ndim;
 			jpt=icol*ndim;
@@ -253,28 +234,28 @@
 			det=-det;
 		}
 
-/*  divide pivot row by pivot element, noting that the original
-	matrix will have 1 on the diagonal, which will be discarded,
-	and the augmented matrix will start with 1 from the identity
-	matrix and then have 1/pivot, which is part of the inverse.  */
+		/*  divide pivot row by pivot element, noting that the original
+			 matrix will have 1 on the diagonal, which will be discarded,
+			 and the augmented matrix will start with 1 from the identity
+			 matrix and then have 1/pivot, which is part of the inverse.  */
 
 		a[ipiv*ndim+ipiv]=1.;
 
 		ipt=ipiv*ndim;
 		for (k=0; k<ncol; k++)
-			a[ipt+k]/=pivot;
+		 a[ipt+k]/=pivot;
 
 		ipt=ipiv*nvec;
 		for (k=0; k<nvec; k++)
-			b[ipt+k]/=pivot;
+		 b[ipt+k]/=pivot;
 
-/*  reduce non-pivot rows such that they will have 0 in the pivot
-	column, which will be discarded, and the augmented matrix will
-	start with 0 from the identity matrix and then have non-zero
-	in the corresponding column, which is part of the inverse.
-	only one column of the augmented matrix is populated at a time,
-	which corresponds to the only column of the original matrix
-	being zeroed, so that the inverse may be done in place.  */
+		/*  reduce non-pivot rows such that they will have 0 in the pivot
+			 column, which will be discarded, and the augmented matrix will
+			 start with 0 from the identity matrix and then have non-zero
+			 in the corresponding column, which is part of the inverse.
+			 only one column of the augmented matrix is populated at a time,
+			 which corresponds to the only column of the original matrix
+			 being zeroed, so that the inverse may be done in place.  */
 
 		for (j=0; j<nrow; j++) {
 			if (j == ipiv) continue;
@@ -286,25 +267,25 @@
 				ipt=j   *ndim;
 				jpt=ipiv*ndim;
 				for (k=0; k<ncol; k++)
-					a[ipt+k]-=dtemp*a[jpt+k];
+				 a[ipt+k]-=dtemp*a[jpt+k];
 
 				ipt=j   *nvec;
 				jpt=ipiv*nvec;
 				for (k=0; k<nvec; k++)
-					b[ipt+k]-=dtemp*b[jpt+k];
+				 b[ipt+k]-=dtemp*b[jpt+k];
 			}
 		}
 
-/*  for a diagonal matrix, the determinant is the product of the
-	diagonal terms, and so it may be accumulated from the pivots,
-	noting that switching rows changes the sign as above  */
+		/*  for a diagonal matrix, the determinant is the product of the
+			 diagonal terms, and so it may be accumulated from the pivots,
+			 noting that switching rows changes the sign as above  */
 
 		det*=pivot;
 	}
 
-/*  switch columns back in reverse order, noting that a row switch
-	in the original matrix corresponds to a column switch in the
-	inverse matrix  */
+	/*  switch columns back in reverse order, noting that a row switch
+		 in the original matrix corresponds to a column switch in the
+		 inverse matrix  */
 
 	for (i=0; i<nrow; i++) {
 		j=(nrow-1)-i;
@@ -313,8 +294,8 @@
 			irow=pivrc1[j];
 			icol=pivrc2[j];
 
-//			_printf_(true,"column switch back for j=%d: irow=%d, icol=%d\n",
-//					 j,irow,icol);
+			//			_printf_(true,"column switch back for j=%d: irow=%d, icol=%d\n",
+			//					 j,irow,icol);
 
 			ipt=0;
 			for (k=0; k<nrow; k++) {
@@ -362,7 +343,7 @@
 	/*Compute determinant of a 3x3 matrix*/
 
 	/*det = a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g)*/
-   *Adet= A[0]*A[4]*A[8]-A[0]*A[5]*A[7]-A[3]*A[1]*A[8]+A[3]*A[2]*A[7]+A[6]*A[1]*A[5]-A[6]*A[2]*A[4];
+	*Adet= A[0]*A[4]*A[8]-A[0]*A[5]*A[7]-A[3]*A[1]*A[8]+A[3]*A[2]*A[7]+A[6]*A[1]*A[5]-A[6]*A[2]*A[4];
 }
 /*}}}*/
 /*FUNCTION Matrix3x3Invert(IssmDouble* Ainv,IssmDouble* A) {{{*/