source: issm/oecreview/Archive/23185-23389/ISSM-23266-23267.diff@ 23390

Last change on this file since 23390 was 23390, checked in by Mathieu Morlighem, 6 years ago

CHG: added Archive/23185-23389
File size: 14.8 KB

  • ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp

     
    102102                #if defined(_HAVE_AD_)
    103103
    104104        if(isautodiff){
    105 #if _HAVE_ADOLC_
    106                 /*Copy some parameters from IoModel to parameters dataset: {{{*/
     105                #if defined(_HAVE_ADOLC_)
     106                /*Copy some parameters from IoModel to parameters dataset*/
    107107                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.obufsize",AutodiffObufsizeEnum));
    108108                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.cbufsize",AutodiffCbufsizeEnum));
    109109                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.lbufsize",AutodiffLbufsizeEnum));
    … …  
    110110                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.tbufsize",AutodiffTbufsizeEnum));
    111111                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerRatio",AutodiffGcTriggerRatioEnum));
    112112                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerMaxSize",AutodiffGcTriggerMaxSizeEnum));
    113                 /*}}}*/
    114 #endif
     113
     114                #elif defined(_HAVE_CODIPACK_)
     115                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.tapeAlloc",AutodiffTapeAllocEnum));
     116
     117                #else
     118                _error_("not supported yet");
     119                #endif
     120
    115121                /*retrieve driver: {{{*/
    116122                iomodel->FindConstant(&autodiff_driver,"md.autodiff.driver");
    117123                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.driver",AutodiffDriverEnum));

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

     
    1010#include "../modules/modules.h"
    1111#include "../solutionsequences/solutionsequences.h"
    1212
    13 #if defined (_HAVE_M1QN3_) && defined(_HAVE_ADOLC_)
     13#ifdef _HAVE_CODIPACK_
     14extern CoDi_global codi_global;
     15#include <sstream> // for output of the CoDiPack tape
     16#endif
     17
     18#if defined (_HAVE_M1QN3_) && defined(_HAVE_AD_)
    1419/*m1qn3 prototypes*/
    1520extern "C" void *ctonbe_; // DIS mode : Conversion
    1621extern "C" void *ctcabe_; // DIS mode : Conversion
    … …  
    3439/*m1qm3 functions*/
    3540void simul_starttrace(FemModel* femmodel){/*{{{*/
    3641
     42        #if defined(_HAVE_ADOLC_)
    3743        /*Retrive ADOLC parameters*/
    3844        IssmDouble gcTriggerRatio;
    3945        IssmDouble gcTriggerMaxSize;
    … …  
    5662        int keepTaylors=1;
    5763        int my_rank=IssmComm::GetRank();
    5864        trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
     65
     66        #elif defined(_HAVE_CODIPACK_)
     67
     68                //fprintf(stderr, "*** Codipack IoModel::StartTrace\n");
     69                /*
     70                 * FIXME codi
     71                 * - ADOL-C variant uses fine grained tracing with various arguments
     72                 * - ADOL-C variant sets a garbage collection parameter for its tape
     73                 * -> These parameters are not read for the CoDiPack ISSM version!
     74                 */
     75                auto& tape_codi = IssmDouble::getGlobalTape();
     76                tape_codi.setActive();
     77                #if _AD_TAPE_ALLOC_
     78                //alloc_profiler.Tag(StartInit, true);
     79                IssmDouble x_t(1.0), y_t(1.0);
     80                tape_codi.registerInput(y_t);
     81                int codi_allocn = 0;
     82                femmodel->parameters->FindParam(&codi_allocn,AutodiffTapeAllocEnum);
     83                for(int i = 0;i < codi_allocn;++i) {
     84                        x_t = y_t * y_t;
     85                }
     86                /*
     87                std::stringstream out_s;
     88                IssmDouble::getGlobalTape().printStatistics(out_s);
     89                _printf0_("StartTrace::Tape Statistics     : TapeAlloc count=[" << codi_allocn << "]\n" << out_s.str());
     90                */
     91                tape_codi.reset();
     92                //alloc_profiler.Tag(FinishInit, true);
     93                #endif
     94
     95        #else
     96        _error_("not implemented");
     97        #endif
    5998}/*}}}*/
     99void simul_stoptrace(){/*{{{*/
     100
     101        #if defined(_HAVE_ADOLC_)
     102        trace_off();
     103        if(VerboseAutodiff()){ /*{{{*/
     104
     105                #ifdef _HAVE_ADOLC_
     106                int my_rank=IssmComm::GetRank();
     107                size_t  tape_stats[15];
     108                tapestats(my_rank,tape_stats); //reading of tape statistics
     109                int commSize=IssmComm::GetSize();
     110                int *sstats=new int[7];
     111                sstats[0]=tape_stats[NUM_OPERATIONS];
     112                sstats[1]=tape_stats[OP_FILE_ACCESS];
     113                sstats[2]=tape_stats[NUM_LOCATIONS];
     114                sstats[3]=tape_stats[LOC_FILE_ACCESS];
     115                sstats[4]=tape_stats[NUM_VALUES];
     116                sstats[5]=tape_stats[VAL_FILE_ACCESS];
     117                sstats[6]=tape_stats[TAY_STACK_SIZE];
     118                int *rstats=NULL;
     119                if (my_rank==0) rstats=new int[commSize*7];
     120                ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
     121                if (my_rank==0) {
     122                        int offset=50;
     123                        int rOffset=(commSize/10)+1;
     124                        _printf_("   ADOLC statistics: \n");
     125                        _printf_("     "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
     126                        _printf_("     "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
     127                        _printf_("     "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
     128                        _printf_("     operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
     129                        _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
     130                        for (int r=0;r<commSize;++r)
     131                         _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+0] << (rstats[r*7+1]?" ->file":"") << "\n");
     132                        _printf_("     locations: entry size " << sizeof(locint) << " Bytes\n");
     133                        _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
     134                        for (int r=0;r<commSize;++r)
     135                         _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+2] << (rstats[r*7+3]?" ->file":"") << "\n");
     136                        _printf_("     constant values: entry size " << sizeof(double) << " Bytes\n");
     137                        _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
     138                        for (int r=0;r<commSize;++r)
     139                         _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+4] << (rstats[r*7+5]?" ->file":"") << "\n");
     140                        _printf_("     Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
     141                        _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
     142                        for (int r=0;r<commSize;++r)
     143                         _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+6] << (rstats[r*7+6]>tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n");
     144                        delete []rstats;
     145                }
     146                delete [] sstats;
     147                #endif
     148
     149                #ifdef _HAVE_CODIPACK_
     150                #ifdef _AD_TAPE_ALLOC_
     151                //_printf_("Allocation time  P(" << my_rank << "): " << alloc_profiler.DeltaTime(StartInit, FinishInit) << "\n");
     152                #endif
     153                std::stringstream out_s;
     154                IssmDouble::getGlobalTape().printStatistics(out_s);
     155                _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str());
     156                #endif
     157        } /*}}}*/
     158
     159        #elif defined(_HAVE_CODIPACK_)
     160        auto& tape_codi = IssmDouble::getGlobalTape();
     161        tape_codi.setPassive();
     162        if(VerboseAutodiff()){
     163                int my_rank=IssmComm::GetRank();
     164                if(my_rank == 0) {
     165                        // FIXME codi "just because" for now
     166                        tape_codi.printStatistics(std::cout);
     167                        codi_global.print(std::cout);
     168                }
     169        }
     170        #else
     171        _error_("not implemented");
     172        #endif
     173}/*}}}*/
    60174void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
    61175
    62176        /*Get rank*/
    … …  
    72186        int N_add = 0;
    73187        int* control_enum = NULL;
    74188
    75         if (solution_type == TransientSolutionEnum){
    76                 femmodel = input_struct->femmodel->copy();
    77         }
     189        /*In transient, we need to make sure we do not modify femmodel at each iteration, make a copy*/
     190        if(solution_type == TransientSolutionEnum) femmodel = input_struct->femmodel->copy();
    78191
    79192        IssmPDouble*  Jlist        = input_struct->Jlist;
    80193        int           JlistM       = input_struct->M;
    … …  
    116229#else
    117230        IssmDouble* aX=xNew<IssmDouble>(intn);
    118231#endif
     232
     233        #if defined(_HAVE_ADOLC_)
    119234        if(my_rank==0){
    120235                for(int i=0;i<intn;i++){
    121236                        aX[i]<<=X[i];
    122237                }
    123238        }
     239        #elif defined(_HAVE_CODIPACK_)
     240        auto& tape_codi = IssmDouble::getGlobalTape();
     241        if(my_rank==0){
     242                for (int i=0;i<intn;i++) {
     243                        aX[i]=X[i];
     244                        tape_codi.registerInput(aX[i]);
     245                        codi_global.input_indices.push_back(aX[i].getGradientData());
     246                }
     247        }
     248        #else
     249        _error_("not suppoted");
     250        #endif
    124251
    125252        ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    126253        SetControlInputsFromVectorx(femmodel,aX);
    … …  
    154281                DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
    155282                if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
    156283                if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
    157                 if (my_rank==0) {
     284                if(my_rank==0) {
     285
     286                        #if defined(_HAVE_CODIPACK_)
     287                        tape_codi.registerOutput(output_value);
     288                        dependents[i] = output_value.getValue();
     289                        codi_global.output_indices.push_back(output_value.getGradientData());
     290
     291                        #elif defined(_HAVE_ADOLC_)
    158292                        output_value>>=dependents[i];
     293
     294                        #else
     295                        _error_("not suppoted");
     296                        #endif
    159297                        J+=output_value;
    160298                }
    161299        }
    162300
    163301        /*Turning off trace tape*/
    164         trace_off();
     302        simul_stoptrace();
    165303        //time_t now = time(NULL);
    166304        //if(my_rank==0) _printf_("\nTIME: "<<now<<"\n");
    167305
    168         /*Print tape statistics so that user can kill this run if something is off already:*/
    169         if(VerboseAutodiff()){ /*{{{*/
    170 
    171                 #ifdef _HAVE_ADOLC_
    172                 size_t  tape_stats[15];
    173                 tapestats(my_rank,tape_stats); //reading of tape statistics
    174                 int commSize=IssmComm::GetSize();
    175                 int *sstats=new int[7];
    176                 sstats[0]=tape_stats[NUM_OPERATIONS];
    177                 sstats[1]=tape_stats[OP_FILE_ACCESS];
    178                 sstats[2]=tape_stats[NUM_LOCATIONS];
    179                 sstats[3]=tape_stats[LOC_FILE_ACCESS];
    180                 sstats[4]=tape_stats[NUM_VALUES];
    181                 sstats[5]=tape_stats[VAL_FILE_ACCESS];
    182                 sstats[6]=tape_stats[TAY_STACK_SIZE];
    183                 int *rstats=NULL;
    184                 if (my_rank==0) rstats=new int[commSize*7];
    185                 ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
    186                 if (my_rank==0) {
    187                         int offset=50;
    188                         int rOffset=(commSize/10)+1;
    189                         _printf_("   ADOLC statistics: \n");
    190                         _printf_("     "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
    191                         _printf_("     "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
    192                         _printf_("     "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
    193                         _printf_("     operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
    194                         _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
    195                         for (int r=0;r<commSize;++r)
    196                          _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+0] << (rstats[r*7+1]?" ->file":"") << "\n");
    197                         _printf_("     locations: entry size " << sizeof(locint) << " Bytes\n");
    198                         _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
    199                         for (int r=0;r<commSize;++r)
    200                          _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+2] << (rstats[r*7+3]?" ->file":"") << "\n");
    201                         _printf_("     constant values: entry size " << sizeof(double) << " Bytes\n");
    202                         _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
    203                         for (int r=0;r<commSize;++r)
    204                          _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+4] << (rstats[r*7+5]?" ->file":"") << "\n");
    205                         _printf_("     Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
    206                         _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
    207                         for (int r=0;r<commSize;++r)
    208                          _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+6] << (rstats[r*7+6]>tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n");
    209                         delete []rstats;
    210                 }
    211                 delete [] sstats;
    212                 #endif
    213 
    214                 #ifdef _HAVE_CODIPACK_
    215                 #ifdef _AD_TAPE_ALLOC_
    216                 //_printf_("Allocation time  P(" << my_rank << "): " << alloc_profiler.DeltaTime(StartInit, FinishInit) << "\n");
    217                 #endif
    218                 std::stringstream out_s;
    219                 IssmDouble::getGlobalTape().printStatistics(out_s);
    220                 _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str());
    221                 #endif
    222         } /*}}}*/
    223 
    224306        /*diverse: */
    225307        int  dummy;
    226308        int  num_independents=0;
    … …  
    246328                num_independents = 0;
    247329        }
    248330
     331        #if defined(_HAVE_ADOLC_)
     332        /*Get gradient for ADOLC {{{*/
     333
    249334        /*get the EDF pointer:*/
    250335        ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
    251336
    … …  
    289374                xDelete(weightVectorTimesJac);
    290375                xDelete(aWeightVector);
    291376        }
     377        /*}}}*/
     378        #elif defined(_HAVE_CODIPACK_)
     379        /*Get gradient for CoDiPack{{{*/
     380        if(VerboseAutodiff())_printf0_("   CoDiPack fos_reverse\n");
     381        int     aDepIndex=0; /*FIXME: do we really need this?*/
    292382
     383        /*retrieve direction index: */
     384        femmodel->parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum);
     385        if (my_rank==0) {
     386                if (aDepIndex<0 || aDepIndex>=num_dependents || codi_global.output_indices.size() <= aDepIndex){
     387                        _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]");
     388                }
     389                tape_codi.setGradient(codi_global.output_indices[aDepIndex], 1.0);
     390        }
     391        tape_codi.evaluate();
     392
     393        weightVectorTimesJac=xNew<double>(num_independents);
     394        /*call driver: */
     395        auto in_size = codi_global.input_indices.size();
     396        for(size_t i = 0; i < in_size; ++i) {
     397                weightVectorTimesJac[i] = tape_codi.getGradient(codi_global.input_indices[i]);
     398        }
     399
     400        /*Add to totalgradient: */
     401        totalgradient=xNewZeroInit<IssmPDouble>(num_independents_old);
     402        if(my_rank==0) for(int i=0;i<num_independents;i++) {
     403                totalgradient[i]+=weightVectorTimesJac[i];
     404        }
     405
     406        /*free resources :*/
     407        xDelete(weightVectorTimesJac);
     408        /*}}}*/
     409        #else
     410        _error_("not suppoted");
     411        #endif
     412
    293413        /*Broadcast gradient to other ranks*/
    294414        ISSM_MPI_Bcast(totalgradient,num_independents_old,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
    295415        /*Check size of Jlist to avoid crashes*/
    … …  
    550670}/*}}}*/
    551671
    552672#else
    553 void controladm1qn3_core(FemModel* femmodel){_error_("M1QN3 or ADOLC not installed");}
     673void controladm1qn3_core(FemModel* femmodel){_error_("M1QN3 or ADOLC/CoDiPack not installed");}
    554674#endif //_HAVE_M1QN3_
Note: See TracBrowser for help on using the repository browser.