Changeset 13498

Timestamp:

10/02/12 14:13:50 (12 years ago)

Author:

utke

Message:

CHG implement fos_reverse for the solver and use a view to the solution vector rather than allocating a new one and copying it

File:

• issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp (modified) (3 diffs)

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

@@ -161 +161 @@
 /*}}}*/
 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 +205 @@
         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 +216 @@
         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
+}
+/*}}}*/