Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 23266) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 23267) @@ -102,8 +102,8 @@ #if defined(_HAVE_AD_) if(isautodiff){ -#if _HAVE_ADOLC_ - /*Copy some parameters from IoModel to parameters dataset: {{{*/ + #if defined(_HAVE_ADOLC_) + /*Copy some parameters from IoModel to parameters dataset*/ parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.obufsize",AutodiffObufsizeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.cbufsize",AutodiffCbufsizeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.lbufsize",AutodiffLbufsizeEnum)); @@ -110,8 +110,14 @@ parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.tbufsize",AutodiffTbufsizeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerRatio",AutodiffGcTriggerRatioEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerMaxSize",AutodiffGcTriggerMaxSizeEnum)); - /*}}}*/ -#endif + + #elif defined(_HAVE_CODIPACK_) + parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.tapeAlloc",AutodiffTapeAllocEnum)); + + #else + _error_("not supported yet"); + #endif + /*retrieve driver: {{{*/ iomodel->FindConstant(&autodiff_driver,"md.autodiff.driver"); parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.driver",AutodiffDriverEnum)); Index: ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 23266) +++ ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 23267) @@ -10,7 +10,12 @@ #include "../modules/modules.h" #include "../solutionsequences/solutionsequences.h" -#if defined (_HAVE_M1QN3_) && defined(_HAVE_ADOLC_) +#ifdef _HAVE_CODIPACK_ +extern CoDi_global codi_global; +#include // for output of the CoDiPack tape +#endif + +#if defined (_HAVE_M1QN3_) && defined(_HAVE_AD_) /*m1qn3 prototypes*/ extern "C" void *ctonbe_; // DIS mode : Conversion extern "C" void *ctcabe_; // DIS mode : Conversion @@ -34,6 +39,7 @@ /*m1qm3 functions*/ void simul_starttrace(FemModel* femmodel){/*{{{*/ + #if defined(_HAVE_ADOLC_) /*Retrive ADOLC parameters*/ IssmDouble gcTriggerRatio; IssmDouble gcTriggerMaxSize; @@ -56,7 +62,115 @@ int keepTaylors=1; int my_rank=IssmComm::GetRank(); trace_on(my_rank,keepTaylors,reCast(obufsize),reCast(lbufsize),reCast(cbufsize),reCast(tbufsize),skipFileDeletion); + + #elif defined(_HAVE_CODIPACK_) + + //fprintf(stderr, "*** Codipack IoModel::StartTrace\n"); + /* + * FIXME codi + * - ADOL-C variant uses fine grained tracing with various arguments + * - ADOL-C variant sets a garbage collection parameter for its tape + * -> These parameters are not read for the CoDiPack ISSM version! + */ + auto& tape_codi = IssmDouble::getGlobalTape(); + tape_codi.setActive(); + #if _AD_TAPE_ALLOC_ + //alloc_profiler.Tag(StartInit, true); + IssmDouble x_t(1.0), y_t(1.0); + tape_codi.registerInput(y_t); + int codi_allocn = 0; + femmodel->parameters->FindParam(&codi_allocn,AutodiffTapeAllocEnum); + for(int i = 0;i < codi_allocn;++i) { + x_t = y_t * y_t; + } + /* + std::stringstream out_s; + IssmDouble::getGlobalTape().printStatistics(out_s); + _printf0_("StartTrace::Tape Statistics : TapeAlloc count=[" << codi_allocn << "]\n" << out_s.str()); + */ + tape_codi.reset(); + //alloc_profiler.Tag(FinishInit, true); + #endif + + #else + _error_("not implemented"); + #endif }/*}}}*/ +void simul_stoptrace(){/*{{{*/ + + #if defined(_HAVE_ADOLC_) + trace_off(); + if(VerboseAutodiff()){ /*{{{*/ + + #ifdef _HAVE_ADOLC_ + int my_rank=IssmComm::GetRank(); + size_t tape_stats[15]; + tapestats(my_rank,tape_stats); //reading of tape statistics + int commSize=IssmComm::GetSize(); + int *sstats=new int[7]; + sstats[0]=tape_stats[NUM_OPERATIONS]; + sstats[1]=tape_stats[OP_FILE_ACCESS]; + sstats[2]=tape_stats[NUM_LOCATIONS]; + sstats[3]=tape_stats[LOC_FILE_ACCESS]; + sstats[4]=tape_stats[NUM_VALUES]; + sstats[5]=tape_stats[VAL_FILE_ACCESS]; + sstats[6]=tape_stats[TAY_STACK_SIZE]; + int *rstats=NULL; + if (my_rank==0) rstats=new int[commSize*7]; + ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm()); + if (my_rank==0) { + int offset=50; + int rOffset=(commSize/10)+1; + _printf_(" ADOLC statistics: \n"); + _printf_(" "<file":"") << "\n"); + _printf_(" locations: entry size " << sizeof(locint) << " Bytes\n"); + _printf_(" "<file":"") << "\n"); + _printf_(" constant values: entry size " << sizeof(double) << " Bytes\n"); + _printf_(" "<file":"") << "\n"); + _printf_(" Taylor stack: entry size " << sizeof(revreal) << " Bytes\n"); + _printf_(" "<tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n"); + delete []rstats; + } + delete [] sstats; + #endif + + #ifdef _HAVE_CODIPACK_ + #ifdef _AD_TAPE_ALLOC_ + //_printf_("Allocation time P(" << my_rank << "): " << alloc_profiler.DeltaTime(StartInit, FinishInit) << "\n"); + #endif + std::stringstream out_s; + IssmDouble::getGlobalTape().printStatistics(out_s); + _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str()); + #endif + } /*}}}*/ + + #elif defined(_HAVE_CODIPACK_) + auto& tape_codi = IssmDouble::getGlobalTape(); + tape_codi.setPassive(); + if(VerboseAutodiff()){ + int my_rank=IssmComm::GetRank(); + if(my_rank == 0) { + // FIXME codi "just because" for now + tape_codi.printStatistics(std::cout); + codi_global.print(std::cout); + } + } + #else + _error_("not implemented"); + #endif +}/*}}}*/ void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/ /*Get rank*/ @@ -72,9 +186,8 @@ int N_add = 0; int* control_enum = NULL; - if (solution_type == TransientSolutionEnum){ - femmodel = input_struct->femmodel->copy(); - } + /*In transient, we need to make sure we do not modify femmodel at each iteration, make a copy*/ + if(solution_type == TransientSolutionEnum) femmodel = input_struct->femmodel->copy(); IssmPDouble* Jlist = input_struct->Jlist; int JlistM = input_struct->M; @@ -116,11 +229,25 @@ #else IssmDouble* aX=xNew(intn); #endif + + #if defined(_HAVE_ADOLC_) if(my_rank==0){ for(int i=0;iGetObjectByOffset(i); if(solution_type==TransientSolutionEnum) output_value = dep->GetValue(); if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel); - if (my_rank==0) { + if(my_rank==0) { + + #if defined(_HAVE_CODIPACK_) + tape_codi.registerOutput(output_value); + dependents[i] = output_value.getValue(); + codi_global.output_indices.push_back(output_value.getGradientData()); + + #elif defined(_HAVE_ADOLC_) output_value>>=dependents[i]; + + #else + _error_("not suppoted"); + #endif J+=output_value; } } /*Turning off trace tape*/ - trace_off(); + simul_stoptrace(); //time_t now = time(NULL); //if(my_rank==0) _printf_("\nTIME: "<file":"") << "\n"); - _printf_(" locations: entry size " << sizeof(locint) << " Bytes\n"); - _printf_(" "<file":"") << "\n"); - _printf_(" constant values: entry size " << sizeof(double) << " Bytes\n"); - _printf_(" "<file":"") << "\n"); - _printf_(" Taylor stack: entry size " << sizeof(revreal) << " Bytes\n"); - _printf_(" "<tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n"); - delete []rstats; - } - delete [] sstats; - #endif - - #ifdef _HAVE_CODIPACK_ - #ifdef _AD_TAPE_ALLOC_ - //_printf_("Allocation time P(" << my_rank << "): " << alloc_profiler.DeltaTime(StartInit, FinishInit) << "\n"); - #endif - std::stringstream out_s; - IssmDouble::getGlobalTape().printStatistics(out_s); - _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str()); - #endif - } /*}}}*/ - /*diverse: */ int dummy; int num_independents=0; @@ -246,6 +328,9 @@ num_independents = 0; } + #if defined(_HAVE_ADOLC_) + /*Get gradient for ADOLC {{{*/ + /*get the EDF pointer:*/ ext_diff_fct *anEDF_for_solverx_p=xDynamicCast * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p; @@ -289,7 +374,42 @@ xDelete(weightVectorTimesJac); xDelete(aWeightVector); } + /*}}}*/ + #elif defined(_HAVE_CODIPACK_) + /*Get gradient for CoDiPack{{{*/ + if(VerboseAutodiff())_printf0_(" CoDiPack fos_reverse\n"); + int aDepIndex=0; /*FIXME: do we really need this?*/ + /*retrieve direction index: */ + femmodel->parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum); + if (my_rank==0) { + if (aDepIndex<0 || aDepIndex>=num_dependents || codi_global.output_indices.size() <= aDepIndex){ + _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]"); + } + tape_codi.setGradient(codi_global.output_indices[aDepIndex], 1.0); + } + tape_codi.evaluate(); + + weightVectorTimesJac=xNew(num_independents); + /*call driver: */ + auto in_size = codi_global.input_indices.size(); + for(size_t i = 0; i < in_size; ++i) { + weightVectorTimesJac[i] = tape_codi.getGradient(codi_global.input_indices[i]); + } + + /*Add to totalgradient: */ + totalgradient=xNewZeroInit(num_independents_old); + if(my_rank==0) for(int i=0;i