source: issm/oecreview/Archive/13393-13976/ISSM-13513-13514.diff

Index: ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp   (revision 13513)
+++ ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp   (revision 13514)
@@ -70,113 +70,113 @@
 
 #ifdef _HAVE_ADOLC_
 int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/
-    SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
-    return 0;
+       SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
+       return 0;
 } /*}}}*/
 int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/
 #ifdef _HAVE_GSL_
-//  for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;
-//  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;
-  // the matrix will be modified by LU decomposition. Use gsl_A copy
-  double* Acopy = xNew<double>(m*m);
-  xMemCpy(Acopy,inVal,m*m);
-  /*Initialize gsl matrices and vectors: */
-  gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
-  gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
-  gsl_permutation *perm_p = gsl_permutation_alloc (m);
-  int  signPerm;
-  // factorize
-  gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
-  gsl_vector *gsl_x_p = gsl_vector_alloc (m);
-  // solve for the value
-  gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
-  /*Copy result*/
-  xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
-  gsl_vector_free(gsl_x_p);
-//  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;
-  // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
-  // compute the RHS
-  double* r=xNew<double>(m);
-  for (int i=0; i<m; i++) {
-    r[i]=inDeriv[m*m+i]; // this is db[i]
-    for (int j=0;j<m; j++) {
-      r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
-    }
-  }
-  gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
-  gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
-  gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
-  xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
-  gsl_vector_free(gsl_dx_p);
-  xDelete(r);
-  gsl_permutation_free(perm_p);
-  xDelete(Acopy);
-  #endif
-  return 0;
+       //  for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;
+       //  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;
+       // the matrix will be modified by LU decomposition. Use gsl_A copy
+       double* Acopy = xNew<double>(m*m);
+       xMemCpy(Acopy,inVal,m*m);
+       /*Initialize gsl matrices and vectors: */
+       gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
+       gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
+       gsl_permutation *perm_p = gsl_permutation_alloc (m);
+       int  signPerm;
+       // factorize
+       gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
+       gsl_vector *gsl_x_p = gsl_vector_alloc (m);
+       // solve for the value
+       gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
+       /*Copy result*/
+       xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
+       gsl_vector_free(gsl_x_p);
+       //  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;
+       // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
+       // compute the RHS
+       double* r=xNew<double>(m);
+       for (int i=0; i<m; i++) {
+               r[i]=inDeriv[m*m+i]; // this is db[i]
+               for (int j=0;j<m; j++) {
+                       r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
+               }
+       }
+       gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
+       gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
+       gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
+       xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
+       gsl_vector_free(gsl_dx_p);
+       xDelete(r);
+       gsl_permutation_free(perm_p);
+       xDelete(Acopy);
+#endif
+       return 0;
 } /*}}}*/
 int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/
 #ifdef _HAVE_GSL_
-  // the matrix will be modified by LU decomposition. Use gsl_A copy
-  double* Acopy = xNew<double>(m*m);
-  xMemCpy(Acopy,inVal,m*m);
-  /*Initialize gsl matrices and vectors: */
-  gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
-  gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
-  gsl_permutation *perm_p = gsl_permutation_alloc (m);
-  int  signPerm;
-  // factorize
-  gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
-  gsl_vector *gsl_x_p = gsl_vector_alloc (m);
-  // solve for the value
-  gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
-  /*Copy result*/
-  xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
-  gsl_vector_free(gsl_x_p);
-  // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
-  double* r=xNew<double>(m);
-  gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
-  for (int dir=0;dir<directionCount;++dir) {
-    // compute the RHS
-    for (int i=0; i<m; i++) {
-      r[i]=inDeriv[m*m+i][dir]; // this is db[i]
-      for (int j=0;j<m; j++) {
-        r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
-      }
-    }
-    gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
-    gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
-    // reuse r
-    xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
-    for (int i=0; i<m; i++) {
-      outDeriv[i][dir]=r[i];
-    }
-  }
-  gsl_vector_free(gsl_dx_p);
-  xDelete(r);
-  gsl_permutation_free(perm_p);
-  xDelete(Acopy);
-  #endif
-  return 0;
+       // the matrix will be modified by LU decomposition. Use gsl_A copy
+       double* Acopy = xNew<double>(m*m);
+       xMemCpy(Acopy,inVal,m*m);
+       /*Initialize gsl matrices and vectors: */
+       gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
+       gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
+       gsl_permutation *perm_p = gsl_permutation_alloc (m);
+       int  signPerm;
+       // factorize
+       gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
+       gsl_vector *gsl_x_p = gsl_vector_alloc (m);
+       // solve for the value
+       gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
+       /*Copy result*/
+       xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
+       gsl_vector_free(gsl_x_p);
+       // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
+       double* r=xNew<double>(m);
+       gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
+       for (int dir=0;dir<directionCount;++dir) {
+               // compute the RHS
+               for (int i=0; i<m; i++) {
+                       r[i]=inDeriv[m*m+i][dir]; // this is db[i]
+                       for (int j=0;j<m; j++) {
+                               r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
+                       }
+               }
+               gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
+               gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
+               // reuse r
+               xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
+               for (int i=0; i<m; i++) {
+                       outDeriv[i][dir]=r[i];
+               }
+       }
+       gsl_vector_free(gsl_dx_p);
+       xDelete(r);
+       gsl_permutation_free(perm_p);
+       xDelete(Acopy);
+#endif
+       return 0;
 }
 /*}}}*/
 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;
+       // 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;
 }
 /*}}}*/
 int EDF_fov_reverse_for_solverx(int m, int p, double **dpp_U, int n, double **dpp_Z, double* dp_x, double* dp_y) { /*{{{*/