Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13296)
+++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13297)
@@ -112,4 +112,50 @@
 /*}}}*/
 
+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*n+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;
+}
+/*}}}*/
+
 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
 	// pack inputs to conform to the EDF-prescribed interface