Index: /issm/trunk-jpl/m4/issm_options.m4 =================================================================== --- /issm/trunk-jpl/m4/issm_options.m4 (revision 26555) +++ /issm/trunk-jpl/m4/issm_options.m4 (revision 26556) @@ -967,4 +967,19 @@ dnl CoDiPack libraries and header files if test "x${HAVE_CODIPACK}" == "xyes"; then + + AC_MSG_CHECKING(for CoDiPack version) + CODIPACK_MAJOR=`cat ${CODIPACK_ROOT}/include/codi.hpp | grep "#define CODI_MAJOR_VERSION" | sed 's/#define CODI_MAJOR_VERSION//' | sed 's/ //g'` + CODIPACK_MINOR=`cat ${CODIPACK_ROOT}/include/codi.hpp | grep "#define CODI_MINOR_VERSION" | sed 's/#define CODI_MINOR_VERSION//' | sed 's/ //g'` + if test -z "${CODIPACK_MAJOR}"; then + CODIPACK_MAJOR=`cat ${CODIPACK_ROOT}/include/codi/configure.h | grep "#define CODI_MAJOR_VERSION" | sed 's/#define CODI_MAJOR_VERSION//' | sed 's/ //g'` + CODIPACK_MINOR=`cat ${CODIPACK_ROOT}/include/codi/configure.h | grep "#define CODI_MINOR_VERSION" | sed 's/#define CODI_MINOR_VERSION//' | sed 's/ //g'` + fi + if test -z "${CODIPACK_MAJOR}"; then + AC_MSG_ERROR([Couldn't determine CoDiPack version]) + fi + AC_DEFINE_UNQUOTED([_CODIPACK_MAJOR_], ${CODIPACK_MAJOR}, [CoDiPack version major]) + AC_DEFINE_UNQUOTED([_CODIPACK_MINOR_], ${CODIPACK_MINOR}, [CoDiPack version minor]) + AC_MSG_RESULT([${CODIPACK_MAJOR}.${CODIPACK_MINOR}]) + CODIPACKINCL="-I${CODIPACK_ROOT}/include" AC_DEFINE([_HAVE_CODIPACK_], [1], [with CoDiPack in ISSM src]) Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26555) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26556) @@ -2299,5 +2299,11 @@ #if defined(_HAVE_CODIPACK_) + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif #endif @@ -2311,5 +2317,11 @@ tape_codi.registerOutput(output_value); dependents[i] = output_value.getValue(); + #if _CODIPACK_MAJOR_==2 codi_global.output_indices.push_back(output_value.getIdentifier()); + #elif _CODIPACK_MAJOR_==1 + codi_global.output_indices.push_back(output_value.getGradientData()); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif #else output_value>>=dependents[i]; Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 26555) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 26556) @@ -2063,7 +2063,16 @@ scalar=pscalar; } + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); tape_codi.registerInput(scalar); codi_global.input_indices.push_back(scalar.getIdentifier()); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + tape_codi.registerInput(scalar); + codi_global.input_indices.push_back(scalar.getGradientData()); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + #else if(X){ @@ -2136,10 +2145,24 @@ #if defined(_HAVE_CODIPACK_) // FIXME codi here we just assign instead of using "operator <<=" + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + if(X){ for (int i=0;i These parameters are not read for the CoDiPack ISSM version! */ + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + tape_codi.setActive(); #if _AD_TAPE_ALLOC_ Index: /issm/trunk-jpl/src/c/cores/ad_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/ad_core.cpp (revision 26555) +++ /issm/trunk-jpl/src/c/cores/ad_core.cpp (revision 26556) @@ -326,5 +326,12 @@ /*First, stop tracing: */ + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + tape_codi.setPassive(); Index: /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 26555) +++ /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 26556) @@ -14,4 +14,5 @@ extern CoDi_global codi_global; #include // for output of the CoDiPack tape +#include void transient_ad(FemModel* femmodel); #endif @@ -75,5 +76,12 @@ * -> These parameters are not read for the CoDiPack ISSM version! */ + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + tape_codi.setActive(); #if _AD_TAPE_ALLOC_ @@ -153,5 +161,11 @@ #endif std::stringstream out_s; + #if _CODIPACK_MAJOR_==2 IssmDouble::getTape().printStatistics(out_s); + #elif _CODIPACK_MAJOR_==1 + IssmDouble::getGlobalTape().printStatistics(out_s); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str()); #endif @@ -159,5 +173,14 @@ #elif defined(_HAVE_CODIPACK_) + + /*Get Tape*/ + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + tape_codi.setPassive(); if(VerboseAutodiff()){ @@ -239,5 +262,14 @@ } #elif defined(_HAVE_CODIPACK_) + + /*Get tape*/ + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + codi_global.input_indices.clear(); if(my_rank==0){ @@ -245,5 +277,12 @@ aX[i]=X[i]; tape_codi.registerInput(aX[i]); + #if _CODIPACK_MAJOR_==2 codi_global.input_indices.push_back(aX[i].getIdentifier()); + #elif _CODIPACK_MAJOR_==1 + codi_global.input_indices.push_back(aX[i].getGradientData()); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + } } @@ -285,5 +324,11 @@ tape_codi.registerOutput(output_value); dependents[i] = output_value.getValue(); + #if _CODIPACK_MAJOR_==2 codi_global.output_indices.push_back(output_value.getIdentifier()); + #elif _CODIPACK_MAJOR_==1 + codi_global.output_indices.push_back(output_value.getGradientData()); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif #elif defined(_HAVE_ADOLC_) @@ -377,4 +422,6 @@ tape_codi.setGradient(codi_global.output_indices[dep_index],1.0); } + feclearexcept(FE_ALL_EXCEPT); + feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); tape_codi.evaluate(); @@ -446,4 +493,5 @@ } Gnorm = sqrt(Gnorm); + _assert_(!xIsNan(Gnorm)); /*Print info*/ Index: /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 26555) +++ /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 26556) @@ -50,5 +50,12 @@ * -> These parameters are not read for the CoDiPack ISSM version! */ + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + tape_codi.setActive(); #if _AD_TAPE_ALLOC_ @@ -127,5 +134,11 @@ #endif std::stringstream out_s; + #if _CODIPACK_MAJOR_==2 IssmDouble::getTape().printStatistics(out_s); + #elif _CODIPACK_MAJOR_==1 + IssmDouble::getGlobalTape().printStatistics(out_s); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str()); #endif @@ -133,5 +146,13 @@ #elif defined(_HAVE_CODIPACK_) + + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + tape_codi.setPassive(); if(VerboseAutodiff()){ @@ -191,5 +212,13 @@ simul_starttrace2(femmodel); IssmDouble* aX=xNew(n); + + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + codi_global.input_indices.clear(); if(my_rank==0){ @@ -197,5 +226,12 @@ aX[i]=X0[i]; tape_codi.registerInput(aX[i]); + #if _CODIPACK_MAJOR_==2 codi_global.input_indices.push_back(aX[i].getIdentifier()); + #elif _CODIPACK_MAJOR_==1 + codi_global.input_indices.push_back(aX[i].getGradientData()); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + } } @@ -231,5 +267,12 @@ tape_codi.registerOutput(output_value); dependents[i] = output_value.getValue(); + #if _CODIPACK_MAJOR_==2 codi_global.output_indices.push_back(output_value.getIdentifier()); + #elif _CODIPACK_MAJOR_==1 + codi_global.output_indices.push_back(output_value.getGradientData()); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + J+=output_value; } Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 26555) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 26556) @@ -379,5 +379,11 @@ /*Start tracing*/ + #if _CODIPACK_MAJOR_==2 auto& tape_codi = IssmDouble::getTape(); + #elif _CODIPACK_MAJOR_==1 + auto& tape_codi = IssmDouble::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif tape_codi.setActive(); Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/Marshalling.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/Marshalling.cpp (revision 26555) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/Marshalling.cpp (revision 26556) @@ -76,5 +76,12 @@ this->identifiers = identifiers_in; this->size_max = size_max_in; + #if _CODIPACK_MAJOR_==2 this->tape_codi = &(IssmDouble::getTape()); + #elif _CODIPACK_MAJOR_==1 + this->tape_codi = &(IssmDouble::getGlobalTape()); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + }/*}}}*/ void RegisterInputFunctor::Echo(void){/*{{{*/ @@ -86,5 +93,12 @@ _assert_(!xIsNan(value)); this->tape_codi->registerInput(value); + #if _CODIPACK_MAJOR_==2 this->identifiers[this->double_count] = value.getIdentifier(); + #elif _CODIPACK_MAJOR_==1 + this->identifiers[this->double_count] = value.getGradientData(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + this->double_count++; }/*}}}*/ @@ -95,5 +109,12 @@ _assert_(!xIsNan(value[i])); this->tape_codi->registerInput(value[i]); + #if _CODIPACK_MAJOR_==2 this->identifiers[this->double_count] = value[i].getIdentifier(); + #elif _CODIPACK_MAJOR_==1 + this->identifiers[this->double_count] = value[i].getGradientData(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + this->double_count++; } @@ -103,5 +124,11 @@ RegisterOutputFunctor::RegisterOutputFunctor(void) : MarshallHandle(AD_REGISTEROUTPUT){/*{{{*/ this->double_count = 0; + #if _CODIPACK_MAJOR_==2 this->tape_codi = &(IssmDouble::getTape()); + #elif _CODIPACK_MAJOR_==1 + this->tape_codi = &(IssmDouble::getGlobalTape()); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif }/*}}}*/ void RegisterOutputFunctor::Echo(void){/*{{{*/ @@ -126,7 +153,13 @@ SetAdjointFunctor::SetAdjointFunctor(double* adjoint_in,int size_max_in) : MarshallHandle(AD_SETADJOINT){/*{{{*/ this->double_count = 0; - this->tape_codi = &(IssmDouble::getTape()); this->adjoint = adjoint_in; this->size_max = size_max_in; + #if _CODIPACK_MAJOR_==2 + this->tape_codi = &(IssmDouble::getTape()); + #elif _CODIPACK_MAJOR_==1 + this->tape_codi = &(IssmDouble::getGlobalTape()); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif }/*}}}*/ void SetAdjointFunctor::Echo(void){/*{{{*/ Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/Marshalling.h =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/Marshalling.h (revision 26555) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/Marshalling.h (revision 26556) @@ -108,4 +108,11 @@ };/*}}}*/ #if defined(_HAVE_CODIPACK_) && !defined(_WRAPPERS_) +#if _CODIPACK_MAJOR_==2 +using Tape = typename IssmDouble::Tape; +#elif _CODIPACK_MAJOR_==1 +using Tape = typename IssmDouble::TapeType; +#else +#error "_CODIPACK_MAJOR_ not supported" +#endif class CountDoublesFunctor: public MarshallHandle{ /*{{{*/ @@ -128,5 +135,5 @@ int *identifiers; int size_max; - IssmDouble::Tape *tape_codi; + Tape *tape_codi; public: @@ -142,5 +149,5 @@ private: int double_count; - IssmDouble::Tape* tape_codi; + Tape *tape_codi; public: @@ -155,8 +162,8 @@ private: - int double_count; - int size_max; - IssmDouble::Tape* tape_codi; - double* adjoint; + int double_count; + int size_max; + Tape *tape_codi; + double *adjoint; public: Index: /issm/trunk-jpl/src/c/toolkits/codipack/CoDiPackCommon.hpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/codipack/CoDiPackCommon.hpp (revision 26555) +++ /issm/trunk-jpl/src/c/toolkits/codipack/CoDiPackCommon.hpp (revision 26556) @@ -35,5 +35,11 @@ inline void getVectorGradData(const Real* vec, Data* dataVec, int n) { for(int i = 0; i < n; ++i) { + #if _CODIPACK_MAJOR_==2 dataVec[i]=vec[i].getIdentifier(); + #elif _CODIPACK_MAJOR_==1 + dataVec[i]=vec[i].getGradientData(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif } } @@ -43,5 +49,11 @@ for(int i = 0; i < n; ++i) { pasVec[i]=vec[i].getValue(); + #if _CODIPACK_MAJOR_==2 dataVec[i]=vec[i].getIdentifier(); + #elif _CODIPACK_MAJOR_==1 + dataVec[i]=vec[i].getGradientData(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif } } @@ -50,14 +62,32 @@ inline void getPrimalAndGradData(const Real& value, Passive& pas, Data& data) { pas=value.getValue(); + #if _CODIPACK_MAJOR_==2 data=value.getIdentifier(); + #elif _CODIPACK_MAJOR_==1 + data=value.getGradientData(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif } template inline void registerVector(Real* vec, Data* dataVec, int n) { - typename Real::Tape& tape = Real::getTape(); + #if _CODIPACK_MAJOR_==2 + typename Real::Tape& tape = Real::getTape(); + #elif _CODIPACK_MAJOR_==1 + typename Real::TapeType& tape = Real::getGlobalTape(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif for(int i = 0; i < n; ++i) { tape.registerInput(vec[i]); + #if _CODIPACK_MAJOR_==2 dataVec[i]=vec[i].getIdentifier(); + #elif _CODIPACK_MAJOR_==1 + dataVec[i]=vec[i].getGradientData(); + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif } } Index: /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp (revision 26555) +++ /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp (revision 26556) @@ -254,4 +254,5 @@ #ifdef _HAVE_CODIPACK_ +#if _CODIPACK_MAJOR_==2 using Tape = typename IssmDouble::Tape; using AccessInterface = codi::VectorAccessInterface; @@ -385,4 +386,140 @@ } /*}}}*/ +#elif _CODIPACK_MAJOR_==1 +void SolverxSeq_codi_b(void* tape_in,void* data_in,void* ra) {/*{{{*/ + + /*recast data_in and tape*/ + codi::DataStore* data = (codi::DataStore*)data_in; + //IssmDouble::TapeType& tape = (IssmDouble::TapeType&)tape_in; + IssmDouble::TapeType& tape = IssmDouble::getGlobalTape(); + + IssmDouble::Real* valueATrans; + IssmDouble::GradientData* indexATrans; + IssmDouble::GradientData* indexB; + IssmDouble::Real* valueX; + IssmDouble::GradientData* indexX; + int n; + + data->getData(valueATrans); + data->getData(indexATrans); + data->getData(indexB); + data->getData(valueX); + data->getData(indexX); + data->getData(n); + + + // create the adjoint vector for x and reset the adjoint values on the tape + IssmDouble::GradientValue* adjX = xNew(n); + getVectorAdjoint(tape, indexX, adjX, n); + + IssmDouble::GradientValue* sol = xNew(n); + SolverxSeq(sol, valueATrans, adjX, n); + + updateVectorAdjoint(tape, indexB, sol, n); + for(int i=0; igetData(valueATrans); + data->getData(indexATrans); + data->getData(indexB); + data->getData(valueX); + data->getData(indexX); + data->getData(n); + + xDelete(valueATrans); + xDelete(indexATrans); + xDelete(indexB); + xDelete(valueX); + xDelete(indexX); + + delete data; +} +/*}}}*/ +void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/ + IssmDouble::TapeType& tape = IssmDouble::getGlobalTape(); + codi::DataStore* dataHandler = NULL; + + if(tape.isActive()) { + dataHandler = new codi::DataStore(); + + // create the index vector and the double data for A and B + IssmDouble::Real* valueATrans = xNew(n*n); + IssmDouble::GradientData* indexATrans = xNew(n*n); + + // read the data for matrix in a transposed fashion + for (int i=0; i(n); + getVectorGradData(B, indexB, n); + + dataHandler->addData(valueATrans); + dataHandler->addData(indexATrans); + dataHandler->addData(indexB); + } + + // unpack the primal values from the matrix and the vector + IssmDouble::Real* valueA = xNew(n*n); + IssmDouble::Real* valueB = xNew(n); + // read the data from A and B + getVectorPrimal(A, valueA, n*n); + getVectorPrimal(B, valueB, n); + + // create the placeholder for X and solve the system + IssmDouble::Real* valueX = xNew(n); + SolverxSeq(valueX, valueA, valueB, n); + + // pack the values into x + setVectorPrimal(X, valueX, n); + + if(tape.isActive()) { + // create the index vector X and register x as active variables + IssmDouble::GradientData* indexX = xNew(n); + registerVector(X, indexX, n); + + dataHandler->addData(valueX); + dataHandler->addData(indexX); + + // store other arguments + dataHandler->addData(n); + + tape.pushExternalFunctionHandle(&SolverxSeq_codi_b, dataHandler, &SolverxSeq_codi_delete); + } + else{ + // if the tape is active valueX is stored in the dataHandler and deleted in the reverse sweep + xDelete(valueX); + } + + xDelete(valueB); + xDelete(valueA); +} +/*}}}*/ +#else +#error "_CODIPACK_MAJOR_ not supported" +#endif void DenseGslSolve(/*output*/ IssmDouble** px,/*stiffness matrix:*/ IssmDouble* Kff, int Kff_M, int Kff_N, /*right hand side load vector: */ IssmDouble* pf, int pf_M, Parameters* parameters){ /*{{{*/ @@ -400,3 +537,4 @@ } /*}}}*/ -#endif + +#endif Index: /issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.h (revision 26555) +++ /issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.h (revision 26556) @@ -34,6 +34,16 @@ //using MpiTypes = CoDiMpiTypes; /*New implementation*/ + #if _CODIPACK_MAJOR_==2 #include - using MpiTypes = codi::CoDiMpiTypes; + using MpiTypes = codi::CoDiMpiTypes; + + #elif _CODIPACK_MAJOR_==1 + #include + using MpiTypes = CoDiMpiTypes; + + #else + #error "_CODIPACK_MAJOR_ not supported" + #endif + extern MpiTypes* mpiTypes; #define AMPI_ADOUBLE mpiTypes->MPI_TYPE Index: /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp (revision 26555) +++ /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp (revision 26556) @@ -485,7 +485,7 @@ #ifdef _HAVE_CODIPACK_ +#if _CODIPACK_MAJOR_==2 using Tape = typename IssmDouble::Tape; using AccessInterface = codi::VectorAccessInterface; - void MumpsSolve_codi_b(Tape* tape, void* data_in, AccessInterface* ra) {/*{{{*/ @@ -648,3 +648,171 @@ } /*}}}*/ -#endif +#elif _CODIPACK_MAJOR_==1 +void MumpsSolve_codi_b(void* tape_in,void* data_in,void* ra) {/*{{{*/ + + /*recast data_in and tape*/ + codi::DataStore* data = (codi::DataStore*)data_in; + //IssmDouble::TapeType& tape = (IssmDouble::TapeType&)tape_in; + IssmDouble::TapeType& tape = IssmDouble::getGlobalTape(); + + + IssmDouble::Real* valueATrans; + IssmDouble::GradientData* indexATrans; + int* irnATrans; + int* jcnATrans; + IssmDouble::GradientData* indexB; + IssmDouble::Real* valueX; + IssmDouble::GradientData* indexX; + int n; + int nnz; + int local_nnz; + Parameters* parameters; + + data->getData(valueATrans); + data->getData(indexATrans); + data->getData(irnATrans); + data->getData(jcnATrans); + data->getData(indexB); + data->getData(valueX); + data->getData(indexX); + data->getData(n); + data->getData(nnz); + data->getData(local_nnz); + data->getData(parameters); + + // create the adjoint vector for x and reset the adjoint values on the tape + IssmDouble::GradientValue* adjX = xNew(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); + + // bcast dp_y (the solution of the forward system) + ISSM_MPI_Bcast(valueX,n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); + // bcast the adjoint of the right-hand-side, i.e. this solution + ISSM_MPI_Bcast(adjX,n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); + + for(int i=0; igetData(valueATrans); + data->getData(indexATrans); + data->getData(irnATrans); + data->getData(jcnATrans); + data->getData(indexB); + data->getData(valueX); + data->getData(indexX); + data->getData(n); + data->getData(nnz); + data->getData(local_nnz); + data->getData(parameters); + + xDelete(valueATrans); + xDelete(indexATrans); + xDelete(irnATrans); + xDelete(jcnATrans); + xDelete(indexB); + xDelete(valueX); + xDelete(indexX); + + delete data; +} +/*}}}*/ +void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){/*{{{*/ + IssmDouble::TapeType& tape = IssmDouble::getGlobalTape(); + codi::DataStore* dataHandler = NULL; + + if(tape.isActive()) { + dataHandler = new codi::DataStore(); + + // create the index and double vector for the matrix + IssmDouble::Real* valueATrans = xNew(local_nnz); + IssmDouble::GradientData* indexATrans = xNew(local_nnz); + int* irnATrans = xNew(local_nnz); + int* jcnATrans = xNew(local_nnz); + + // read the data for the matrix A in a transposed fashion + for (int i=0; i(n); + getVectorGradData(rhs, indexB, n); + + dataHandler->addData(valueATrans); + dataHandler->addData(indexATrans); + dataHandler->addData(irnATrans); + dataHandler->addData(jcnATrans); + dataHandler->addData(indexB); + } + + // unpack the primal values from the matrix and the vector + IssmDouble::Real* valueA = xNew(local_nnz); + IssmDouble::Real* valueB = xNew(n); + // read the data from A and B + getVectorPrimal(a_loc, valueA, local_nnz); + getVectorPrimal(rhs, valueB, n); + + MumpsSolve(n, nnz, local_nnz, irn_loc, jcn_loc, valueA, valueB, parameters); + // valueB contains now the solution + + // pack the values into rhs + setVectorPrimal(rhs, valueB, n); + + if(tape.isActive()) { + // create the index vector X and register x as active variables + IssmDouble::GradientData* indexX = xNew(n); + registerVector(rhs, indexX, n); + + dataHandler->addData(valueB); // contains the values from x + dataHandler->addData(indexX); + + // store other arguments + dataHandler->addData(n); + dataHandler->addData(nnz); + dataHandler->addData(local_nnz); + 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); + } + else{ + // if the tape is active valueB is stored in the dataHandler and deleted in the reverse sweep + xDelete(valueB); + } + + xDelete(valueA); +} +/*}}}*/ +#else +#error "_CODIPACK_MAJOR_ not supported" +#endif +#endif