Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13294)
+++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13295)
@@ -66,9 +66,50 @@
 
 #ifdef _HAVE_ADOLC_
-int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){/*{{{*/
-    if(m*(m+1)!=n)_error_("Stiffness matrix should be square!");
-    SolverxSeq(y,x, x+m*m, m);
+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;
-}/*}}}*/
+}
+/*}}}*/
+
+int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, 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
+  // 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*n+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;
+}
+/*}}}*/
+
 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
 	// pack inputs to conform to the EDF-prescribed interface