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

../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp

@@ -179 +179 @@
   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: