Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13511)
+++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13512)
@@ -180,4 +180,26 @@
 }
 /*}}}*/
+int EDF_fov_reverse_for_solverx(int m, int p, double **dpp_U, int n, double **dpp_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);
+	gsl_permutation *perm_p = gsl_permutation_alloc (m);
+	int permSign;
+	gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
+	for (int weightsRowIndex=0;weightsRowIndex<p;++weightsRowIndex) {
+		// the adjoint of the solution is our right-hand side
+		gsl_vector_view x_bar=gsl_vector_view_array(dpp_U[weightsRowIndex],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(dpp_Z[weightsRowIndex]+m*m,m);
+		gsl_linalg_LU_solve (&aTransposed.matrix, perm_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) dpp_Z[weightsRowIndex][i*m+j]-=dpp_Z[weightsRowIndex][m*m+i]*dp_y[j];
+	}
+	gsl_permutation_free(perm_p);
+	xDelete(transposed);
+	return 0;
+}
+/*}}}*/
 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
 	// pack inputs to conform to the EDF-prescribed interface