Changeset 26556


Ignore:
Timestamp:
11/09/21 05:34:35 (3 years ago)
Author:
Mathieu Morlighem
Message:

CHG: can now run with CoDiPack 1.9 and 2.0

Location:
issm/trunk-jpl
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/m4/issm_options.m4

    r26392 r26556  
    967967        dnl CoDiPack libraries and header files
    968968        if test "x${HAVE_CODIPACK}" == "xyes"; then
     969
     970                AC_MSG_CHECKING(for CoDiPack version)
     971                CODIPACK_MAJOR=`cat ${CODIPACK_ROOT}/include/codi.hpp | grep "#define CODI_MAJOR_VERSION" | sed 's/#define CODI_MAJOR_VERSION//' | sed 's/ //g'`
     972                CODIPACK_MINOR=`cat ${CODIPACK_ROOT}/include/codi.hpp | grep "#define CODI_MINOR_VERSION" | sed 's/#define CODI_MINOR_VERSION//' | sed 's/ //g'`
     973                if test -z "${CODIPACK_MAJOR}"; then
     974                         CODIPACK_MAJOR=`cat ${CODIPACK_ROOT}/include/codi/configure.h | grep "#define CODI_MAJOR_VERSION" | sed 's/#define CODI_MAJOR_VERSION//' | sed 's/ //g'`
     975                         CODIPACK_MINOR=`cat ${CODIPACK_ROOT}/include/codi/configure.h | grep "#define CODI_MINOR_VERSION" | sed 's/#define CODI_MINOR_VERSION//' | sed 's/ //g'`
     976                fi
     977                if test -z "${CODIPACK_MAJOR}"; then
     978                        AC_MSG_ERROR([Couldn't determine CoDiPack version])
     979                fi
     980                AC_DEFINE_UNQUOTED([_CODIPACK_MAJOR_], ${CODIPACK_MAJOR}, [CoDiPack version major])
     981                AC_DEFINE_UNQUOTED([_CODIPACK_MINOR_], ${CODIPACK_MINOR}, [CoDiPack version minor])
     982                AC_MSG_RESULT([${CODIPACK_MAJOR}.${CODIPACK_MINOR}])
     983
    969984                CODIPACKINCL="-I${CODIPACK_ROOT}/include"
    970985                AC_DEFINE([_HAVE_CODIPACK_], [1], [with CoDiPack in ISSM src])

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r26487 r26556  
    22992299
    23002300                        #if defined(_HAVE_CODIPACK_)
     2301                        #if _CODIPACK_MAJOR_==2
    23012302                        auto& tape_codi = IssmDouble::getTape();
     2303                        #elif _CODIPACK_MAJOR_==1
     2304                        auto& tape_codi = IssmDouble::getGlobalTape();
     2305                        #else
     2306                        #error "_CODIPACK_MAJOR_ not supported"
     2307                        #endif
    23022308                        #endif
    23032309
    … …  
    23112317                                                tape_codi.registerOutput(output_value);
    23122318                                                dependents[i] = output_value.getValue();
     2319                                                #if _CODIPACK_MAJOR_==2
    23132320                                                codi_global.output_indices.push_back(output_value.getIdentifier());
     2321                                                #elif _CODIPACK_MAJOR_==1
     2322                                                codi_global.output_indices.push_back(output_value.getGradientData());
     2323                                                #else
     2324                                                #error "_CODIPACK_MAJOR_ not supported"
     2325                                                #endif
    23142326                                        #else
    23152327                                                output_value>>=dependents[i];

  • issm/trunk-jpl/src/c/classes/IoModel.cpp

    r26468 r26556  
    20632063                                scalar=pscalar;
    20642064                        }
     2065                        #if _CODIPACK_MAJOR_==2
    20652066                        auto& tape_codi = IssmDouble::getTape();
    20662067                        tape_codi.registerInput(scalar);
    20672068                        codi_global.input_indices.push_back(scalar.getIdentifier());
     2069                        #elif _CODIPACK_MAJOR_==1
     2070                        auto& tape_codi = IssmDouble::getGlobalTape();
     2071                        tape_codi.registerInput(scalar);
     2072                        codi_global.input_indices.push_back(scalar.getGradientData());
     2073                        #else
     2074                        #error "_CODIPACK_MAJOR_ not supported"
     2075                        #endif
     2076
    20682077                #else
    20692078                        if(X){
    … …  
    21362145                        #if defined(_HAVE_CODIPACK_)
    21372146                                // FIXME codi here we just assign instead of using "operator <<="
     2147                                #if _CODIPACK_MAJOR_==2
    21382148                                auto& tape_codi = IssmDouble::getTape();
     2149                                #elif _CODIPACK_MAJOR_==1
     2150                                auto& tape_codi = IssmDouble::getGlobalTape();
     2151                                #else
     2152                                #error "_CODIPACK_MAJOR_ not supported"
     2153                                #endif
     2154
    21392155                                if(X){
    21402156                                        for (int i=0;i<M*N;i++) {
    21412157                                                matrix[i]=X[Xcount+i];
    21422158                                                tape_codi.registerInput(matrix[i]);
     2159                                                #if _CODIPACK_MAJOR_==2
    21432160                                                codi_global.input_indices.push_back(matrix[i].getIdentifier());
     2161                                                #elif _CODIPACK_MAJOR_==1
     2162                                                codi_global.input_indices.push_back(matrix[i].getGradientData());
     2163                                                #else
     2164                                                #error "_CODIPACK_MAJOR_ not supported"
     2165                                                #endif
     2166
    21442167                                        }
    21452168                                }
    … …  
    21482171                                                matrix[i]=buffer[i];
    21492172                                                tape_codi.registerInput(matrix[i]);
     2173                                                #if _CODIPACK_MAJOR_==2
    21502174                                                codi_global.input_indices.push_back(matrix[i].getIdentifier());
     2175                                                #elif _CODIPACK_MAJOR_==1
     2176                                                codi_global.input_indices.push_back(matrix[i].getGradientData());
     2177                                                #else
     2178                                                #error "_CODIPACK_MAJOR_ not supported"
     2179                                                #endif
    21512180                                        }
    21522181                                }
    … …  
    29612990                 * -> These parameters are not read for the CoDiPack ISSM version!
    29622991                 */
     2992                #if _CODIPACK_MAJOR_==2
    29632993                auto& tape_codi = IssmDouble::getTape();
     2994                #elif _CODIPACK_MAJOR_==1
     2995                auto& tape_codi = IssmDouble::getGlobalTape();
     2996                #else
     2997                #error "_CODIPACK_MAJOR_ not supported"
     2998                #endif
     2999
    29643000                tape_codi.setActive();
    29653001                #if _AD_TAPE_ALLOC_

  • issm/trunk-jpl/src/c/cores/ad_core.cpp

    r26432 r26556  
    326326
    327327                        /*First, stop tracing: */
     328                        #if _CODIPACK_MAJOR_==2
    328329                        auto& tape_codi = IssmDouble::getTape();
     330                        #elif _CODIPACK_MAJOR_==1
     331                        auto& tape_codi = IssmDouble::getGlobalTape();
     332                        #else
     333                        #error "_CODIPACK_MAJOR_ not supported"
     334                        #endif
     335
    329336                        tape_codi.setPassive();
    330337

  • issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp

    r26523 r26556  
    1414extern CoDi_global codi_global;
    1515#include <sstream> // for output of the CoDiPack tape
     16#include <fenv.h>
    1617void transient_ad(FemModel* femmodel);
    1718#endif
    … …  
    7576                 * -> These parameters are not read for the CoDiPack ISSM version!
    7677                 */
     78                #if _CODIPACK_MAJOR_==2
    7779                auto& tape_codi = IssmDouble::getTape();
     80                #elif _CODIPACK_MAJOR_==1
     81                auto& tape_codi = IssmDouble::getGlobalTape();
     82                #else
     83                #error "_CODIPACK_MAJOR_ not supported"
     84                #endif
     85
    7886                tape_codi.setActive();
    7987                #if _AD_TAPE_ALLOC_
    … …  
    153161                #endif
    154162                std::stringstream out_s;
     163                #if _CODIPACK_MAJOR_==2
    155164                IssmDouble::getTape().printStatistics(out_s);
     165                #elif _CODIPACK_MAJOR_==1
     166                IssmDouble::getGlobalTape().printStatistics(out_s);
     167                #else
     168                #error "_CODIPACK_MAJOR_ not supported"
     169                #endif
    156170                _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str());
    157171                #endif
    … …  
    159173
    160174        #elif defined(_HAVE_CODIPACK_)
     175
     176        /*Get Tape*/
     177        #if _CODIPACK_MAJOR_==2
    161178        auto& tape_codi = IssmDouble::getTape();
     179        #elif _CODIPACK_MAJOR_==1
     180        auto& tape_codi = IssmDouble::getGlobalTape();
     181        #else
     182        #error "_CODIPACK_MAJOR_ not supported"
     183        #endif
     184
    162185        tape_codi.setPassive();
    163186        if(VerboseAutodiff()){
    … …  
    239262        }
    240263        #elif defined(_HAVE_CODIPACK_)
     264
     265        /*Get tape*/
     266        #if _CODIPACK_MAJOR_==2
    241267        auto& tape_codi = IssmDouble::getTape();
     268        #elif _CODIPACK_MAJOR_==1
     269        auto& tape_codi = IssmDouble::getGlobalTape();
     270        #else
     271        #error "_CODIPACK_MAJOR_ not supported"
     272        #endif
     273
    242274        codi_global.input_indices.clear();
    243275        if(my_rank==0){
    … …  
    245277                        aX[i]=X[i];
    246278                        tape_codi.registerInput(aX[i]);
     279                        #if _CODIPACK_MAJOR_==2
    247280                        codi_global.input_indices.push_back(aX[i].getIdentifier());
     281                        #elif _CODIPACK_MAJOR_==1
     282                        codi_global.input_indices.push_back(aX[i].getGradientData());
     283                        #else
     284                        #error "_CODIPACK_MAJOR_ not supported"
     285                        #endif
     286
    248287                }
    249288        }
    … …  
    285324                        tape_codi.registerOutput(output_value);
    286325                        dependents[i] = output_value.getValue();
     326                        #if _CODIPACK_MAJOR_==2
    287327                        codi_global.output_indices.push_back(output_value.getIdentifier());
     328                        #elif _CODIPACK_MAJOR_==1
     329                        codi_global.output_indices.push_back(output_value.getGradientData());
     330                        #else
     331                        #error "_CODIPACK_MAJOR_ not supported"
     332                        #endif
    288333
    289334                        #elif defined(_HAVE_ADOLC_)
    … …  
    377422                        tape_codi.setGradient(codi_global.output_indices[dep_index],1.0);
    378423                }
     424                feclearexcept(FE_ALL_EXCEPT);
     425                feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
    379426                tape_codi.evaluate();
    380427
    … …  
    446493        }
    447494        Gnorm = sqrt(Gnorm);
     495        _assert_(!xIsNan(Gnorm));
    448496
    449497        /*Print info*/

  • issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp

    r26478 r26556  
    5050                 * -> These parameters are not read for the CoDiPack ISSM version!
    5151                 */
     52                #if _CODIPACK_MAJOR_==2
    5253                auto& tape_codi = IssmDouble::getTape();
     54                #elif _CODIPACK_MAJOR_==1
     55                auto& tape_codi = IssmDouble::getGlobalTape();
     56                #else
     57                #error "_CODIPACK_MAJOR_ not supported"
     58                #endif
     59
    5360                tape_codi.setActive();
    5461                #if _AD_TAPE_ALLOC_
    … …  
    127134                #endif
    128135                std::stringstream out_s;
     136                #if _CODIPACK_MAJOR_==2
    129137                IssmDouble::getTape().printStatistics(out_s);
     138                #elif _CODIPACK_MAJOR_==1
     139                IssmDouble::getGlobalTape().printStatistics(out_s);
     140                #else
     141                #error "_CODIPACK_MAJOR_ not supported"
     142                #endif
    130143                _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str());
    131144                #endif
    … …  
    133146
    134147        #elif defined(_HAVE_CODIPACK_)
     148
     149        #if _CODIPACK_MAJOR_==2
    135150        auto& tape_codi = IssmDouble::getTape();
     151        #elif _CODIPACK_MAJOR_==1
     152        auto& tape_codi = IssmDouble::getGlobalTape();
     153        #else
     154        #error "_CODIPACK_MAJOR_ not supported"
     155        #endif
     156
    136157        tape_codi.setPassive();
    137158        if(VerboseAutodiff()){
    … …  
    191212        simul_starttrace2(femmodel);
    192213        IssmDouble* aX=xNew<IssmDouble>(n);
     214
     215        #if _CODIPACK_MAJOR_==2
    193216        auto& tape_codi = IssmDouble::getTape();
     217        #elif _CODIPACK_MAJOR_==1
     218        auto& tape_codi = IssmDouble::getGlobalTape();
     219        #else
     220        #error "_CODIPACK_MAJOR_ not supported"
     221        #endif
     222
    194223        codi_global.input_indices.clear();
    195224        if(my_rank==0){
    … …  
    197226                        aX[i]=X0[i];
    198227                        tape_codi.registerInput(aX[i]);
     228                        #if _CODIPACK_MAJOR_==2
    199229                        codi_global.input_indices.push_back(aX[i].getIdentifier());
     230                        #elif _CODIPACK_MAJOR_==1
     231                        codi_global.input_indices.push_back(aX[i].getGradientData());
     232                        #else
     233                        #error "_CODIPACK_MAJOR_ not supported"
     234                        #endif
     235
    200236                }
    201237        }
    … …  
    231267                        tape_codi.registerOutput(output_value);
    232268                        dependents[i] = output_value.getValue();
     269                        #if _CODIPACK_MAJOR_==2
    233270                        codi_global.output_indices.push_back(output_value.getIdentifier());
     271                        #elif _CODIPACK_MAJOR_==1
     272                        codi_global.output_indices.push_back(output_value.getGradientData());
     273                        #else
     274                        #error "_CODIPACK_MAJOR_ not supported"
     275                        #endif
     276
    234277                        J+=output_value;
    235278                }

  • issm/trunk-jpl/src/c/cores/transient_core.cpp

    r26526 r26556  
    379379
    380380        /*Start tracing*/
     381        #if _CODIPACK_MAJOR_==2
    381382        auto& tape_codi = IssmDouble::getTape();
     383        #elif _CODIPACK_MAJOR_==1
     384        auto& tape_codi = IssmDouble::getGlobalTape();
     385        #else
     386        #error "_CODIPACK_MAJOR_ not supported"
     387        #endif
    382388        tape_codi.setActive();
    383389

  • issm/trunk-jpl/src/c/shared/io/Marshalling/Marshalling.cpp

    r26534 r26556  
    7676        this->identifiers  = identifiers_in;
    7777        this->size_max         = size_max_in;
     78        #if _CODIPACK_MAJOR_==2
    7879        this->tape_codi    = &(IssmDouble::getTape());
     80        #elif _CODIPACK_MAJOR_==1
     81        this->tape_codi    = &(IssmDouble::getGlobalTape());
     82        #else
     83        #error "_CODIPACK_MAJOR_ not supported"
     84        #endif
     85
    7986}/*}}}*/
    8087void RegisterInputFunctor::Echo(void){/*{{{*/
    … …  
    8693        _assert_(!xIsNan<IssmDouble>(value));
    8794        this->tape_codi->registerInput(value);
     95        #if _CODIPACK_MAJOR_==2
    8896        this->identifiers[this->double_count] = value.getIdentifier();
     97        #elif _CODIPACK_MAJOR_==1
     98        this->identifiers[this->double_count] = value.getGradientData();
     99        #else
     100        #error "_CODIPACK_MAJOR_ not supported"
     101        #endif
     102
    89103        this->double_count++;
    90104}/*}}}*/
    … …  
    95109                        _assert_(!xIsNan<IssmDouble>(value[i]));
    96110                        this->tape_codi->registerInput(value[i]);
     111                        #if _CODIPACK_MAJOR_==2
    97112                        this->identifiers[this->double_count] = value[i].getIdentifier();
     113                        #elif _CODIPACK_MAJOR_==1
     114                        this->identifiers[this->double_count] = value[i].getGradientData();
     115                        #else
     116                        #error "_CODIPACK_MAJOR_ not supported"
     117                        #endif
     118
    98119                        this->double_count++;
    99120                }
    … …  
    103124RegisterOutputFunctor::RegisterOutputFunctor(void) : MarshallHandle(AD_REGISTEROUTPUT){/*{{{*/
    104125        this->double_count = 0;
     126        #if _CODIPACK_MAJOR_==2
    105127        this->tape_codi    = &(IssmDouble::getTape());
     128        #elif _CODIPACK_MAJOR_==1
     129        this->tape_codi    = &(IssmDouble::getGlobalTape());
     130        #else
     131        #error "_CODIPACK_MAJOR_ not supported"
     132        #endif
    106133}/*}}}*/
    107134void RegisterOutputFunctor::Echo(void){/*{{{*/
    … …  
    126153SetAdjointFunctor::SetAdjointFunctor(double* adjoint_in,int size_max_in) : MarshallHandle(AD_SETADJOINT){/*{{{*/
    127154        this->double_count = 0;
    128         this->tape_codi    = &(IssmDouble::getTape());
    129155        this->adjoint      = adjoint_in;
    130156        this->size_max     = size_max_in;
     157        #if _CODIPACK_MAJOR_==2
     158        this->tape_codi    = &(IssmDouble::getTape());
     159        #elif _CODIPACK_MAJOR_==1
     160        this->tape_codi    = &(IssmDouble::getGlobalTape());
     161        #else
     162        #error "_CODIPACK_MAJOR_ not supported"
     163        #endif
    131164}/*}}}*/
    132165void SetAdjointFunctor::Echo(void){/*{{{*/

  • issm/trunk-jpl/src/c/shared/io/Marshalling/Marshalling.h

    r26432 r26556  
    108108};/*}}}*/
    109109#if defined(_HAVE_CODIPACK_) && !defined(_WRAPPERS_)
     110#if _CODIPACK_MAJOR_==2
     111using Tape = typename IssmDouble::Tape;
     112#elif _CODIPACK_MAJOR_==1
     113using Tape = typename IssmDouble::TapeType;
     114#else
     115#error "_CODIPACK_MAJOR_ not supported"
     116#endif
    110117class CountDoublesFunctor:    public MarshallHandle{ /*{{{*/
    111118
    … …  
    128135                int *identifiers;
    129136                int  size_max;
    130                 IssmDouble::Tape *tape_codi;
     137                Tape *tape_codi;
    131138
    132139        public:
    … …  
    142149        private:
    143150                int   double_count;
    144                 IssmDouble::Tape* tape_codi;
     151                Tape *tape_codi;
    145152
    146153        public:
    … …  
    155162
    156163        private:
    157                 int                   double_count;
    158                 int                   size_max;
    159                 IssmDouble::Tape* tape_codi;
    160                 double*               adjoint;
     164                int     double_count;
     165                int     size_max;
     166                Tape   *tape_codi;
     167                double *adjoint;
    161168
    162169        public:

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

    r26433 r26556  
    3535inline void getVectorGradData(const Real* vec, Data* dataVec, int n) {
    3636  for(int i = 0; i < n; ++i) {
     37         #if _CODIPACK_MAJOR_==2
    3738    dataVec[i]=vec[i].getIdentifier();
     39         #elif _CODIPACK_MAJOR_==1
     40         dataVec[i]=vec[i].getGradientData();
     41         #else
     42         #error "_CODIPACK_MAJOR_ not supported"
     43         #endif
    3844  }
    3945}
    … …  
    4349  for(int i = 0; i < n; ++i) {
    4450    pasVec[i]=vec[i].getValue();
     51         #if _CODIPACK_MAJOR_==2
    4552    dataVec[i]=vec[i].getIdentifier();
     53         #elif _CODIPACK_MAJOR_==1
     54         dataVec[i]=vec[i].getGradientData();
     55         #else
     56         #error "_CODIPACK_MAJOR_ not supported"
     57         #endif
    4658  }
    4759}
    … …  
    5062inline void getPrimalAndGradData(const Real& value, Passive& pas, Data& data) {
    5163  pas=value.getValue();
     64  #if _CODIPACK_MAJOR_==2
    5265  data=value.getIdentifier();
     66  #elif _CODIPACK_MAJOR_==1
     67  data=value.getGradientData();
     68  #else
     69  #error "_CODIPACK_MAJOR_ not supported"
     70  #endif
    5371}
    5472
    5573template<typename Real, typename Data>
    5674inline void registerVector(Real* vec, Data* dataVec, int n) {
    57   typename Real::Tape& tape = Real::getTape();
     75        #if _CODIPACK_MAJOR_==2
     76        typename Real::Tape& tape = Real::getTape();
     77        #elif _CODIPACK_MAJOR_==1
     78        typename Real::TapeType& tape = Real::getGlobalTape();
     79        #else
     80        #error "_CODIPACK_MAJOR_ not supported"
     81        #endif
    5882
    5983  for(int i = 0; i < n; ++i) {
    6084    tape.registerInput(vec[i]);
     85         #if _CODIPACK_MAJOR_==2
    6186    dataVec[i]=vec[i].getIdentifier();
     87         #elif _CODIPACK_MAJOR_==1
     88         dataVec[i]=vec[i].getGradientData();
     89         #else
     90         #error "_CODIPACK_MAJOR_ not supported"
     91         #endif
    6292  }
    6393}

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

    r26468 r26556  
    254254
    255255#ifdef _HAVE_CODIPACK_
     256#if _CODIPACK_MAJOR_==2
    256257using Tape = typename IssmDouble::Tape;
    257258using AccessInterface = codi::VectorAccessInterface<typename Tape::Real, typename Tape::Identifier>;
    … …  
    385386}
    386387/*}}}*/
     388#elif _CODIPACK_MAJOR_==1
     389void SolverxSeq_codi_b(void* tape_in,void* data_in,void* ra) {/*{{{*/
     390
     391        /*recast data_in and tape*/
     392        codi::DataStore* data = (codi::DataStore*)data_in;
     393        //IssmDouble::TapeType& tape = (IssmDouble::TapeType&)tape_in;
     394        IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
     395
     396        IssmDouble::Real* valueATrans;
     397        IssmDouble::GradientData* indexATrans;
     398        IssmDouble::GradientData* indexB;
     399        IssmDouble::Real* valueX;
     400        IssmDouble::GradientData* indexX;
     401        int n;
     402
     403        data->getData(valueATrans);
     404        data->getData(indexATrans);
     405        data->getData(indexB);
     406        data->getData(valueX);
     407        data->getData(indexX);
     408        data->getData(n);
     409
     410
     411        // create the adjoint vector for x and reset the adjoint values on the tape
     412        IssmDouble::GradientValue* adjX = xNew<IssmDouble::GradientValue>(n);
     413        getVectorAdjoint(tape, indexX, adjX, n);
     414
     415        IssmDouble::GradientValue* sol  = xNew<IssmDouble::GradientValue>(n);
     416        SolverxSeq(sol, valueATrans, adjX, n);
     417
     418        updateVectorAdjoint(tape, indexB, sol, n);
     419        for(int i=0; i<n; ++i) {
     420                for (int j=0; j<n; ++j) {
     421                        // we access the transposed matrix here because we stored the indices in a transposed way
     422                        updateAdjoint(tape, indexATrans[i*n+j], -sol[j]*valueX[i]);
     423                }
     424        }
     425
     426        xDelete(sol);
     427        xDelete(adjX);
     428}
     429/*}}}*/
     430void SolverxSeq_codi_delete(void* tape_in,void* data_in) {/*{{{*/
     431
     432        /*recast data_in*/
     433        codi::DataStore* data = (codi::DataStore*)data_in;
     434
     435        IssmDouble::Real* valueATrans;
     436        IssmDouble::GradientData* indexATrans;
     437        IssmDouble::GradientData* indexB;
     438        IssmDouble::Real* valueX;
     439        IssmDouble::GradientData* indexX;
     440        int n;
     441
     442        data->getData(valueATrans);
     443        data->getData(indexATrans);
     444        data->getData(indexB);
     445        data->getData(valueX);
     446        data->getData(indexX);
     447        data->getData(n);
     448
     449        xDelete(valueATrans);
     450        xDelete(indexATrans);
     451        xDelete(indexB);
     452        xDelete(valueX);
     453        xDelete(indexX);
     454
     455        delete data;
     456}
     457/*}}}*/
     458void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
     459        IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
     460        codi::DataStore* dataHandler = NULL;
     461
     462        if(tape.isActive()) {
     463                dataHandler = new codi::DataStore();
     464
     465                // create the index vector and the double data for A and B
     466                IssmDouble::Real* valueATrans = xNew<IssmDouble::Real>(n*n);
     467                IssmDouble::GradientData* indexATrans = xNew<IssmDouble::GradientData>(n*n);
     468
     469                // read the data for matrix in a transposed fashion
     470                for (int i=0; i<n; ++i) {
     471                        for (int j=0; j<n; ++j) {
     472                                getPrimalAndGradData(A[i*n+j], valueATrans[j*n+i], indexATrans[j*n+i]);
     473                        }
     474                }
     475
     476                // read the data from B (primal values are not required vor B
     477                IssmDouble::GradientData* indexB = xNew<IssmDouble::GradientData>(n);
     478                getVectorGradData(B, indexB, n);
     479
     480                dataHandler->addData(valueATrans);
     481                dataHandler->addData(indexATrans);
     482                dataHandler->addData(indexB);
     483        }
     484
     485        // unpack the primal values from the matrix and the vector
     486        IssmDouble::Real* valueA = xNew<IssmDouble::Real>(n*n);
     487        IssmDouble::Real* valueB = xNew<IssmDouble::Real>(n);
     488        // read the data from A and B
     489        getVectorPrimal(A, valueA, n*n);
     490        getVectorPrimal(B, valueB, n);
     491
     492        // create the placeholder for X and solve the system
     493        IssmDouble::Real* valueX = xNew<IssmDouble::Real>(n);
     494        SolverxSeq(valueX, valueA, valueB, n);
     495
     496        // pack the values into x
     497        setVectorPrimal(X, valueX, n);
     498
     499        if(tape.isActive()) {
     500                // create the index vector X and register x as active variables
     501                IssmDouble::GradientData* indexX = xNew<IssmDouble::GradientData>(n);
     502                registerVector(X, indexX, n);
     503
     504                dataHandler->addData(valueX);
     505                dataHandler->addData(indexX);
     506
     507                // store other arguments
     508                dataHandler->addData(n);
     509
     510                tape.pushExternalFunctionHandle(&SolverxSeq_codi_b, dataHandler, &SolverxSeq_codi_delete);
     511        }
     512        else{
     513                // if the tape is active valueX is stored in the dataHandler and deleted in the reverse sweep
     514                xDelete(valueX);
     515        }
     516
     517        xDelete(valueB);
     518        xDelete(valueA);
     519}
     520/*}}}*/
     521#else
     522#error "_CODIPACK_MAJOR_ not supported"
     523#endif
    387524void 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){ /*{{{*/
    388525
    … …  
    400537}
    401538/*}}}*/
    402 #endif
     539
     540#endif

  • issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.h

    r26432 r26556  
    3434                        //using MpiTypes = CoDiMpiTypes<IssmDouble>;
    3535                        /*New implementation*/
     36                        #if _CODIPACK_MAJOR_==2
    3637                        #include <codi/tools/mpi/codiMpiTypes.hpp>
    37          using MpiTypes = codi::CoDiMpiTypes<IssmDouble>;
     38                        using MpiTypes = codi::CoDiMpiTypes<IssmDouble>;
     39
     40                        #elif _CODIPACK_MAJOR_==1
     41                        #include <codi/externals/codiMpiTypes.hpp>
     42                        using MpiTypes = CoDiMpiTypes<IssmDouble>;
     43
     44                        #else
     45                        #error "_CODIPACK_MAJOR_ not supported"
     46                        #endif
     47
    3848                        extern MpiTypes* mpiTypes;
    3949                        #define AMPI_ADOUBLE mpiTypes->MPI_TYPE

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

    r26433 r26556  
    485485
    486486#ifdef _HAVE_CODIPACK_
     487#if _CODIPACK_MAJOR_==2
    487488using Tape = typename IssmDouble::Tape;
    488489using AccessInterface = codi::VectorAccessInterface<typename Tape::Real, typename Tape::Identifier>;
    489 
    490490void MumpsSolve_codi_b(Tape* tape, void* data_in, AccessInterface* ra) {/*{{{*/
    491491
    … …  
    648648}
    649649/*}}}*/
    650 #endif
     650#elif _CODIPACK_MAJOR_==1
     651void MumpsSolve_codi_b(void* tape_in,void* data_in,void* ra) {/*{{{*/
     652
     653        /*recast data_in and tape*/
     654        codi::DataStore* data = (codi::DataStore*)data_in;
     655        //IssmDouble::TapeType& tape = (IssmDouble::TapeType&)tape_in;
     656        IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
     657
     658
     659        IssmDouble::Real* valueATrans;
     660        IssmDouble::GradientData* indexATrans;
     661        int* irnATrans;
     662        int* jcnATrans;
     663        IssmDouble::GradientData* indexB;
     664        IssmDouble::Real* valueX;
     665        IssmDouble::GradientData* indexX;
     666        int n;
     667        int nnz;
     668        int local_nnz;
     669        Parameters* parameters;
     670
     671        data->getData(valueATrans);
     672        data->getData(indexATrans);
     673        data->getData(irnATrans);
     674        data->getData(jcnATrans);
     675        data->getData(indexB);
     676        data->getData(valueX);
     677        data->getData(indexX);
     678        data->getData(n);
     679        data->getData(nnz);
     680        data->getData(local_nnz);
     681        data->getData(parameters);
     682
     683        // create the adjoint vector for x and reset the adjoint values on the tape
     684        IssmDouble::GradientValue* adjX = xNew<IssmDouble::GradientValue>(n);
     685        getVectorAdjoint(tape, indexX, adjX, n);
     686
     687        MumpsSolve(n, nnz, local_nnz, irnATrans, jcnATrans, valueATrans, adjX, parameters);
     688        // adjX contains now the solution
     689
     690        updateVectorAdjoint(tape, indexB, adjX, n);
     691
     692        // bcast dp_y (the solution of the forward system)
     693        ISSM_MPI_Bcast(valueX,n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
     694        // bcast the adjoint of the right-hand-side, i.e. this solution
     695        ISSM_MPI_Bcast(adjX,n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
     696
     697        for(int i=0; i<local_nnz; ++i) {
     698                // we access the transposed matrix here because we stored the indices in a transposed way
     699                // -1 is substracted because jcn and irn are stored with fortran indexing
     700                updateAdjoint(tape, indexATrans[i], -adjX[jcnATrans[i]-1]*valueX[irnATrans[i]-1]);
     701        }
     702
     703        xDelete(adjX);
     704}
     705/*}}}*/
     706void MumpsSolve_codi_delete(void* tape_in,void* data_in) {/*{{{*/
     707
     708        /*recast data_in*/
     709        codi::DataStore* data = (codi::DataStore*)data_in;
     710
     711        IssmDouble::Real* valueATrans;
     712        IssmDouble::GradientData* indexATrans;
     713        int* irnATrans;
     714        int* jcnATrans;
     715        IssmDouble::GradientData* indexB;
     716        IssmDouble::Real* valueX;
     717        IssmDouble::GradientData* indexX;
     718        int n;
     719        int nnz;
     720        int local_nnz;
     721        Parameters* parameters;
     722
     723        data->getData(valueATrans);
     724        data->getData(indexATrans);
     725        data->getData(irnATrans);
     726        data->getData(jcnATrans);
     727        data->getData(indexB);
     728        data->getData(valueX);
     729        data->getData(indexX);
     730        data->getData(n);
     731        data->getData(nnz);
     732        data->getData(local_nnz);
     733        data->getData(parameters);
     734
     735        xDelete(valueATrans);
     736        xDelete(indexATrans);
     737        xDelete(irnATrans);
     738        xDelete(jcnATrans);
     739        xDelete(indexB);
     740        xDelete(valueX);
     741        xDelete(indexX);
     742
     743        delete data;
     744}
     745/*}}}*/
     746void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){/*{{{*/
     747        IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
     748        codi::DataStore* dataHandler = NULL;
     749
     750        if(tape.isActive()) {
     751                dataHandler = new codi::DataStore();
     752
     753                // create the index and double vector for the matrix
     754                IssmDouble::Real* valueATrans = xNew<IssmDouble::Real>(local_nnz);
     755                IssmDouble::GradientData* indexATrans = xNew<IssmDouble::GradientData>(local_nnz);
     756                int* irnATrans = xNew<int>(local_nnz);
     757                int* jcnATrans = xNew<int>(local_nnz);
     758
     759                // read the data for the matrix A in a transposed fashion
     760                for (int i=0; i<local_nnz; ++i) {
     761                        getPrimalAndGradData(a_loc[i], valueATrans[i], indexATrans[i]);
     762                        irnATrans[i]=jcn_loc[i];  // transposed store
     763                        jcnATrans[i]=irn_loc[i];  // transposed store
     764                }
     765
     766                // create the index vector for a (primal values are not needed for a)
     767                IssmDouble::GradientData* indexB = xNew<IssmDouble::GradientData>(n);
     768                getVectorGradData(rhs, indexB, n);
     769
     770                dataHandler->addData(valueATrans);
     771                dataHandler->addData(indexATrans);
     772                dataHandler->addData(irnATrans);
     773                dataHandler->addData(jcnATrans);
     774                dataHandler->addData(indexB);
     775        }
     776
     777        // unpack the primal values from the matrix and the vector
     778        IssmDouble::Real* valueA = xNew<IssmDouble::Real>(local_nnz);
     779        IssmDouble::Real* valueB = xNew<IssmDouble::Real>(n);
     780        // read the data from A and B
     781        getVectorPrimal(a_loc, valueA, local_nnz);
     782        getVectorPrimal(rhs, valueB, n);
     783
     784        MumpsSolve(n, nnz, local_nnz, irn_loc, jcn_loc, valueA, valueB, parameters);
     785        // valueB contains now the solution
     786
     787        // pack the values into rhs
     788        setVectorPrimal(rhs, valueB, n);
     789
     790        if(tape.isActive()) {
     791                // create the index vector X and register x as active variables
     792                IssmDouble::GradientData* indexX = xNew<IssmDouble::GradientData>(n);
     793                registerVector(rhs, indexX, n);
     794
     795                dataHandler->addData(valueB); // contains the values from x
     796                dataHandler->addData(indexX);
     797
     798                // store other arguments
     799                dataHandler->addData(n);
     800                dataHandler->addData(nnz);
     801                dataHandler->addData(local_nnz);
     802                dataHandler->addData(parameters); // we assume here that parameters is still intact when the reverse run is called
     803
     804                //tape.pushExternalFunction(&MumpsSolve_codi_b, dataHandler, &MumpsSolve_codi_delete);
     805                tape.pushExternalFunctionHandle(&MumpsSolve_codi_b,(void*)dataHandler, &MumpsSolve_codi_delete);
     806        }
     807        else{
     808                // if the tape is active valueB is stored in the dataHandler and deleted in the reverse sweep
     809                xDelete(valueB);
     810        }
     811
     812        xDelete(valueA);
     813}
     814/*}}}*/
     815#else
     816#error "_CODIPACK_MAJOR_ not supported"
     817#endif
     818#endif
Note: See TracChangeset for help on using the changeset viewer.