Changeset 26433

Timestamp:

09/06/21 17:31:41 (4 years ago)

Author:

Mathieu Morlighem

Message:

CHG: transitioning to CoDiPack v2.0

Location:

issm/trunk-jpl/src/c/toolkits

Files:

• codipack/CoDiPackCommon.hpp (modified) (1 diff)
• gsl/DenseGslSolve.cpp (modified) (7 diffs)
• mumps/MumpsSolve.cpp (modified) (8 diffs)

issm/trunk-jpl/src/c/toolkits/codipack/CoDiPackCommon.hpp

@@ -55 +55 @@
 template<typename Real, typename Data>
 inline void registerVector(Real* vec, Data* dataVec, int n) {
-  typename Real::TapeType& tape = Real::getTape();
+  typename Real::Tape& tape = Real::getTape();
 
   for(int i = 0; i < n; ++i) {

issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp

@@ -254 +254 @@
 
 #ifdef _HAVE_CODIPACK_
-void SolverxSeq_codi_b(void* tape_in,void* data_in,void* ra) {/*{{{*/
+using Tape = typename IssmDouble::Tape;
+using AccessInterface = codi::VectorAccessInterface<typename Tape::Real, typename Tape::Identifier>;
+void SolverxSeq_codi_b(Tape* tape,void* data_in, AccessInterface* ra) {/*{{{*/
 
         /*recast data_in and tape*/
-        codi::DataStore* data = (codi::DataStore*)data_in;
-        //IssmDouble::Tape& tape = (IssmDouble::Tape&)tape_in;
-        IssmDouble::Tape& tape = IssmDouble::getTape();
+        codi::ExternalFunctionUserData* data = (codi::ExternalFunctionUserData*)data_in;
 
   IssmDouble::Real* valueATrans;
-  IssmDouble::GradientData* indexATrans;
-  IssmDouble::GradientData* indexB;
+  IssmDouble::Identifier* indexATrans;
+  IssmDouble::Identifier* indexB;
   IssmDouble::Real* valueX;
-  IssmDouble::GradientData* indexX;
+  IssmDouble::Identifier* indexX;
   int n;
 
@@ -275 +275 @@
   data->getData(n);
 
-
   // create the adjoint vector for x and reset the adjoint values on the tape
   IssmDouble::Gradient* adjX = xNew<IssmDouble::Gradient>(n);
-  getVectorAdjoint(tape, indexX, adjX, n);
+  getVectorAdjoint(*tape, indexX, adjX, n);
 
   IssmDouble::Gradient* sol  = xNew<IssmDouble::Gradient>(n);
   SolverxSeq(sol, valueATrans, adjX, n);
 
-  updateVectorAdjoint(tape, indexB, sol, n);
+  updateVectorAdjoint(*tape, indexB, sol, n);
   for(int i=0; i<n; ++i) {
     for (int j=0; j<n; ++j) {
       // we access the transposed matrix here because we stored the indices in a transposed way
-      updateAdjoint(tape, indexATrans[i*n+j], -sol[j]*valueX[i]);
+      updateAdjoint(*tape, indexATrans[i*n+j], -sol[j]*valueX[i]);
     }
   }
@@ -295 +294 @@
 }
 /*}}}*/
-void SolverxSeq_codi_delete(void* tape_in,void* data_in) {/*{{{*/
+void SolverxSeq_codi_delete(Tape* tape,void* data_in) {/*{{{*/
 
         /*recast data_in*/
-        codi::DataStore* data = (codi::DataStore*)data_in;
+        codi::ExternalFunctionUserData* data = (codi::ExternalFunctionUserData*)data_in;
 
   IssmDouble::Real* valueATrans;
-  IssmDouble::GradientData* indexATrans;
-  IssmDouble::GradientData* indexB;
+  IssmDouble::Identifier* indexATrans;
+  IssmDouble::Identifier* indexB;
   IssmDouble::Real* valueX;
-  IssmDouble::GradientData* indexX;
+  IssmDouble::Identifier* indexX;
   int n;
 
@@ -325 +324 @@
 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
   IssmDouble::Tape& tape = IssmDouble::getTape();
-  codi::DataStore* dataHandler = NULL;
+  codi::ExternalFunctionUserData* dataHandler = NULL;
 
   if(tape.isActive()) {
-    dataHandler = new codi::DataStore();
+    dataHandler = new codi::ExternalFunctionUserData();
 
     // create the index vector and the double data for A and B
     IssmDouble::Real* valueATrans = xNew<IssmDouble::Real>(n*n);
-    IssmDouble::GradientData* indexATrans = xNew<IssmDouble::GradientData>(n*n);
+    IssmDouble::Identifier* indexATrans = xNew<IssmDouble::Identifier>(n*n);
 
     // read the data for matrix in a transposed fashion
@@ -342 +341 @@
 
     // read the data from B (primal values are not required vor B
-    IssmDouble::GradientData* indexB = xNew<IssmDouble::GradientData>(n);
+    IssmDouble::Identifier* indexB = xNew<IssmDouble::Identifier>(n);
     getVectorGradData(B, indexB, n);
 
@@ -366 +365 @@
   if(tape.isActive()) {
     // create the index vector X and register x as active variables
-    IssmDouble::GradientData* indexX = xNew<IssmDouble::GradientData>(n);
+    IssmDouble::Identifier* indexX = xNew<IssmDouble::Identifier>(n);
     registerVector(X, indexX, n);
 
@@ -375 +374 @@
     dataHandler->addData(n);
 
-    tape.pushExternalFunctionHandle(&SolverxSeq_codi_b, dataHandler, &SolverxSeq_codi_delete);
+         tape.pushExternalFunction(codi::ExternalFunction<Tape>::create(&SolverxSeq_codi_b,(void*)dataHandler, &SolverxSeq_codi_delete));
   }
   else{

issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp

@@ -485 +485 @@
 
 #ifdef _HAVE_CODIPACK_
-void MumpsSolve_codi_b(void* tape_in,void* data_in,void* ra) {/*{{{*/
+using Tape = typename IssmDouble::Tape;
+using AccessInterface = codi::VectorAccessInterface<typename Tape::Real, typename Tape::Identifier>;
+
+void MumpsSolve_codi_b(Tape* tape, void* data_in, AccessInterface* ra) {/*{{{*/
 
         /*recast data_in and tape*/
-  codi::DataStore* data = (codi::DataStore*)data_in;
-  //IssmDouble::Tape& tape = (IssmDouble::Tape&)tape_in;
-  IssmDouble::Tape& tape = IssmDouble::getTape();
-
+  codi::ExternalFunctionUserData* data = (codi::ExternalFunctionUserData*)data_in;
 
   IssmDouble::Real* valueATrans;
-  IssmDouble::GradientData* indexATrans;
+  IssmDouble::Identifier* indexATrans;
   int* irnATrans;
   int* jcnATrans;
-  IssmDouble::GradientData* indexB;
+  IssmDouble::Identifier* indexB;
   IssmDouble::Real* valueX;
-  IssmDouble::GradientData* indexX;
+  IssmDouble::Identifier* indexX;
   int n;
   int nnz;
@@ -519 +519 @@
   // create the adjoint vector for x and reset the adjoint values on the tape
   IssmDouble::Gradient* adjX = xNew<IssmDouble::Gradient>(n);
-  getVectorAdjoint(tape, indexX, adjX, n);
+  getVectorAdjoint(*tape, indexX, adjX, n);
 
   MumpsSolve(n, nnz, local_nnz, irnATrans, jcnATrans, valueATrans, adjX, parameters);
   // adjX contains now the solution
 
-  updateVectorAdjoint(tape, indexB, adjX, n);
+  updateVectorAdjoint(*tape, indexB, adjX, n);
 
   // bcast dp_y (the solution of the forward system)
@@ -534 +534 @@
     // we access the transposed matrix here because we stored the indices in a transposed way
     // -1 is substracted because jcn and irn are stored with fortran indexing
-    updateAdjoint(tape, indexATrans[i], -adjX[jcnATrans[i]-1]*valueX[irnATrans[i]-1]);
+    updateAdjoint(*tape, indexATrans[i], -adjX[jcnATrans[i]-1]*valueX[irnATrans[i]-1]);
   }
 
@@ -540 +540 @@
 }
 /*}}}*/
-void MumpsSolve_codi_delete(void* tape_in,void* data_in) {/*{{{*/
+void MumpsSolve_codi_delete(Tape* tape,void* data_in) {/*{{{*/
 
         /*recast data_in*/
-        codi::DataStore* data = (codi::DataStore*)data_in;
+        codi::ExternalFunctionUserData* data = (codi::ExternalFunctionUserData*)data_in;
 
   IssmDouble::Real* valueATrans;
-  IssmDouble::GradientData* indexATrans;
+  IssmDouble::Identifier* indexATrans;
   int* irnATrans;
   int* jcnATrans;
-  IssmDouble::GradientData* indexB;
+  IssmDouble::Identifier* indexB;
   IssmDouble::Real* valueX;
-  IssmDouble::GradientData* indexX;
+  IssmDouble::Identifier* indexX;
   int n;
   int nnz;
@@ -582 +582 @@
 void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){/*{{{*/
   IssmDouble::Tape& tape = IssmDouble::getTape();
-  codi::DataStore* dataHandler = NULL;
+  codi::ExternalFunctionUserData* dataHandler = NULL;
 
   if(tape.isActive()) {
-    dataHandler = new codi::DataStore();
+    dataHandler = new codi::ExternalFunctionUserData();
 
     // create the index and double vector for the matrix
     IssmDouble::Real* valueATrans = xNew<IssmDouble::Real>(local_nnz);
-    IssmDouble::GradientData* indexATrans = xNew<IssmDouble::GradientData>(local_nnz);
+    IssmDouble::Identifier* indexATrans = xNew<IssmDouble::Identifier>(local_nnz);
     int* irnATrans = xNew<int>(local_nnz);
     int* jcnATrans = xNew<int>(local_nnz);
@@ -601 +601 @@
 
     // create the index vector for a (primal values are not needed for a)
-    IssmDouble::GradientData* indexB = xNew<IssmDouble::GradientData>(n);
+    IssmDouble::Identifier* indexB = xNew<IssmDouble::Identifier>(n);
     getVectorGradData(rhs, indexB, n);
 
@@ -626 +626 @@
   if(tape.isActive()) {
     // create the index vector X and register x as active variables
-    IssmDouble::GradientData* indexX = xNew<IssmDouble::GradientData>(n);
+    IssmDouble::Identifier* indexX = xNew<IssmDouble::Identifier>(n);
     registerVector(rhs, indexX, n);
 
@@ -638 +638 @@
     dataHandler->addData(parameters); // we assume here that parameters is still intact when the reverse run is called
 
-         //tape.pushExternalFunction(&MumpsSolve_codi_b, dataHandler, &MumpsSolve_codi_delete);
-    tape.pushExternalFunctionHandle(&MumpsSolve_codi_b,(void*)dataHandler, &MumpsSolve_codi_delete);
+         tape.pushExternalFunction(codi::ExternalFunction<Tape>::create(&MumpsSolve_codi_b,(void*)dataHandler, &MumpsSolve_codi_delete));
   }
   else{