Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13497)
+++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13498)
@@ -161,4 +161,20 @@
 /*}}}*/
 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;
 }
@@ -189,6 +205,5 @@
 	int              s;
 	gsl_matrix_view  a;
-	gsl_vector_view  b;
-	gsl_vector      *x = NULL;
+	gsl_vector_view  b,x;
 	gsl_permutation *p = NULL;
 //	for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl;
@@ -201,23 +216,15 @@
 	a = gsl_matrix_view_array (Acopy,n,n);
 	b = gsl_vector_view_array (B,n);
-	x = gsl_vector_alloc (n);
+	x = gsl_vector_view_array (X,n);
 
 	/*Run LU and solve: */
 	p = gsl_permutation_alloc (n);
 	gsl_linalg_LU_decomp (&a.matrix, p, &s);
-	gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
-
-	//printf ("x = \n");
-	//gsl_vector_fprintf (stdout, x, "%g");
-
-	/*Copy result*/
-	xMemCpy(X,gsl_vector_ptr(x,0),n);
-//      for (int i=0; i<n; ++i) std::cout << "SolverxSeq x["<< i << "]=" << X[i] << std::endl;
+	gsl_linalg_LU_solve (&a.matrix, p, &b.vector, &x.vector);
 
 	/*Clean up and assign output pointer*/
 	xDelete(Acopy);
 	gsl_permutation_free(p);
-	gsl_vector_free(x);
-#endif
-}
-/*}}}*/
+#endif
+}
+/*}}}*/