source: issm/oecreview/Archive/13393-13976/ISSM-13511-13512.diff@ 14312

Index: ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp   (revision 13511)
+++ ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp   (revision 13512)
@@ -179,6 +179,28 @@
   return 0;
 }
 /*}}}*/
+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
         // ensure a contiguous block of locations: