Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13513)
+++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13514)
@@ -71,111 +71,111 @@
 #ifdef _HAVE_ADOLC_
 int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/
-    SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
-    return 0;
+	SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
+	return 0;
 } /*}}}*/
 int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/
 #ifdef _HAVE_GSL_
-//  for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;
-//  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;
-  // the matrix will be modified by LU decomposition. Use gsl_A copy
-  double* Acopy = xNew<double>(m*m);
-  xMemCpy(Acopy,inVal,m*m);
-  /*Initialize gsl matrices and vectors: */
-  gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
-  gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
-  gsl_permutation *perm_p = gsl_permutation_alloc (m);
-  int  signPerm;
-  // factorize
-  gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
-  gsl_vector *gsl_x_p = gsl_vector_alloc (m);
-  // solve for the value
-  gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
-  /*Copy result*/
-  xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
-  gsl_vector_free(gsl_x_p);
-//  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;
-  // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
-  // compute the RHS
-  double* r=xNew<double>(m);
-  for (int i=0; i<m; i++) {
-    r[i]=inDeriv[m*m+i]; // this is db[i]
-    for (int j=0;j<m; j++) {
-      r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
-    }
-  }
-  gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
-  gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
-  gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
-  xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
-  gsl_vector_free(gsl_dx_p);
-  xDelete(r);
-  gsl_permutation_free(perm_p);
-  xDelete(Acopy);
-  #endif
-  return 0;
+	//  for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;
+	//  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;
+	// the matrix will be modified by LU decomposition. Use gsl_A copy
+	double* Acopy = xNew<double>(m*m);
+	xMemCpy(Acopy,inVal,m*m);
+	/*Initialize gsl matrices and vectors: */
+	gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
+	gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
+	gsl_permutation *perm_p = gsl_permutation_alloc (m);
+	int  signPerm;
+	// factorize
+	gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
+	gsl_vector *gsl_x_p = gsl_vector_alloc (m);
+	// solve for the value
+	gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
+	/*Copy result*/
+	xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
+	gsl_vector_free(gsl_x_p);
+	//  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;
+	// solve for the derivatives acc. to A * dx = r  with r=db - dA * x
+	// compute the RHS
+	double* r=xNew<double>(m);
+	for (int i=0; i<m; i++) {
+		r[i]=inDeriv[m*m+i]; // this is db[i]
+		for (int j=0;j<m; j++) {
+			r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
+		}
+	}
+	gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
+	gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
+	gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
+	xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
+	gsl_vector_free(gsl_dx_p);
+	xDelete(r);
+	gsl_permutation_free(perm_p);
+	xDelete(Acopy);
+#endif
+	return 0;
 } /*}}}*/
 int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/
 #ifdef _HAVE_GSL_
-  // the matrix will be modified by LU decomposition. Use gsl_A copy
-  double* Acopy = xNew<double>(m*m);
-  xMemCpy(Acopy,inVal,m*m);
-  /*Initialize gsl matrices and vectors: */
-  gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
-  gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
-  gsl_permutation *perm_p = gsl_permutation_alloc (m);
-  int  signPerm;
-  // factorize
-  gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
-  gsl_vector *gsl_x_p = gsl_vector_alloc (m);
-  // solve for the value
-  gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
-  /*Copy result*/
-  xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
-  gsl_vector_free(gsl_x_p);
-  // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
-  double* r=xNew<double>(m);
-  gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
-  for (int dir=0;dir<directionCount;++dir) {
-    // compute the RHS
-    for (int i=0; i<m; i++) {
-      r[i]=inDeriv[m*m+i][dir]; // this is db[i]
-      for (int j=0;j<m; j++) {
-        r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
-      }
-    }
-    gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
-    gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
-    // reuse r
-    xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
-    for (int i=0; i<m; i++) {
-      outDeriv[i][dir]=r[i];
-    }
-  }
-  gsl_vector_free(gsl_dx_p);
-  xDelete(r);
-  gsl_permutation_free(perm_p);
-  xDelete(Acopy);
-  #endif
-  return 0;
+	// the matrix will be modified by LU decomposition. Use gsl_A copy
+	double* Acopy = xNew<double>(m*m);
+	xMemCpy(Acopy,inVal,m*m);
+	/*Initialize gsl matrices and vectors: */
+	gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
+	gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
+	gsl_permutation *perm_p = gsl_permutation_alloc (m);
+	int  signPerm;
+	// factorize
+	gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
+	gsl_vector *gsl_x_p = gsl_vector_alloc (m);
+	// solve for the value
+	gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
+	/*Copy result*/
+	xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
+	gsl_vector_free(gsl_x_p);
+	// solve for the derivatives acc. to A * dx = r  with r=db - dA * x
+	double* r=xNew<double>(m);
+	gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
+	for (int dir=0;dir<directionCount;++dir) {
+		// compute the RHS
+		for (int i=0; i<m; i++) {
+			r[i]=inDeriv[m*m+i][dir]; // this is db[i]
+			for (int j=0;j<m; j++) {
+				r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
+			}
+		}
+		gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
+		gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
+		// reuse r
+		xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
+		for (int i=0; i<m; i++) {
+			outDeriv[i][dir]=r[i];
+		}
+	}
+	gsl_vector_free(gsl_dx_p);
+	xDelete(r);
+	gsl_permutation_free(perm_p);
+	xDelete(Acopy);
+#endif
+	return 0;
 }
 /*}}}*/
 int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z, double* dp_x, double* dp_y) { /*{{{*/
-  // copy to transpose the matrix
-  double* transposed=xNew<double>(m*m);
-  for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
-  gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
-  // the adjoint of the solution is our right-hand side
-  gsl_vector_view x_bar=gsl_vector_view_array(dp_U,m);
-  // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
-  gsl_vector_view b_bar=gsl_vector_view_array(dp_Z+m*m,m);
-  gsl_permutation *p = gsl_permutation_alloc (m);
-  int permSign;
-  gsl_linalg_LU_decomp (&aTransposed.matrix, p, &permSign);
-  gsl_linalg_LU_solve (&aTransposed.matrix, p, &x_bar.vector, &b_bar.vector);
-  // now do the adjoint of the matrix
-  for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) dp_Z[i*m+j]-=dp_Z[m*m+i]*dp_y[j];
-  gsl_permutation_free(p);
-  xDelete(transposed);
-  return 0;
+	// copy to transpose the matrix
+	double* transposed=xNew<double>(m*m);
+	for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
+	gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
+	// the adjoint of the solution is our right-hand side
+	gsl_vector_view x_bar=gsl_vector_view_array(dp_U,m);
+	// the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
+	gsl_vector_view b_bar=gsl_vector_view_array(dp_Z+m*m,m);
+	gsl_permutation *p = gsl_permutation_alloc (m);
+	int permSign;
+	gsl_linalg_LU_decomp (&aTransposed.matrix, p, &permSign);
+	gsl_linalg_LU_solve (&aTransposed.matrix, p, &x_bar.vector, &b_bar.vector);
+	// now do the adjoint of the matrix
+	for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) dp_Z[i*m+j]-=dp_Z[m*m+i]*dp_y[j];
+	gsl_permutation_free(p);
+	xDelete(transposed);
+	return 0;
 }
 /*}}}*/