Index: ../trunk-jpl/src/c/Makefile.am =================================================================== --- ../trunk-jpl/src/c/Makefile.am (revision 18895) +++ ../trunk-jpl/src/c/Makefile.am (revision 18896) @@ -345,6 +345,7 @@ ./cores/WrapperCorePointerFromSolutionEnum.cpp\ ./cores/CorePointerFromSolutionEnum.cpp\ ./cores/ad_core.cpp\ + ./cores/adgradient_core.cpp\ ./main/EnvironmentInit.cpp\ ./main/EnvironmentFinalize.cpp\ ./analyses/EnumToAnalysis.h\ Index: ../trunk-jpl/src/c/cores/adgradient_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/adgradient_core.cpp (revision 0) +++ ../trunk-jpl/src/c/cores/adgradient_core.cpp (revision 18896) @@ -0,0 +1,145 @@ +/*!\file adgradient_core + * \brief: compute gradient for all scalar depenendents, then sum them up as output. This relies mainly on the fos_reverse + * driver, hence the incapacity to merge this with ad_core.cpp. + */ + +/*Includes: {{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" +/*}}}*/ + +void adgradient_core(FemModel* femmodel){ + + /*diverse: */ + int i; + int dummy; + int num_dependents=0; + int num_dependents_old=0; + int num_independents=0; + bool isautodiff = false; + int aDepIndex=0; + int my_rank=IssmComm::GetRank(); + + /*state variables: */ + IssmDouble *axp = NULL; + IssmPDouble *xp = NULL; + + /*intermediary: */ + IssmPDouble *aWeightVector=NULL; + IssmPDouble *weightVectorTimesJac=NULL; + + /*output: */ + IssmPDouble *totalgradient=NULL; + + /*AD mode on?: */ + femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); + + if(isautodiff){ + + #ifdef _HAVE_ADOLC_ + if(VerboseAutodiff())_printf0_(" start ad core\n"); + + /*First, stop tracing: */ + trace_off(); + + /*retrieve parameters: */ + femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); + femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum); + + /*if no dependents, no point in running a driver: */ + if(!(num_dependents*num_independents)) return; + + /*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/ + num_dependents_old=num_dependents; + if (my_rank!=0){ + num_dependents=0; num_independents=0; + } + + /*retrieve state variable: */ + femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum); + + /* driver argument */ + xp=xNew(num_independents); + for(i=0;i(axp[i]); + } + + /*get the EDF pointer:*/ + ext_diff_fct *anEDF_for_solverx_p=xDynamicCast * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p; + + /* these are always needed regardless of the interpreter */ + anEDF_for_solverx_p->dp_x=xNew(anEDF_for_solverx_p->max_n); + anEDF_for_solverx_p->dp_y=xNew(anEDF_for_solverx_p->max_m); + + /* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so + * as to generate num_dependents gradients: */ + + /*Initialize outputs: */ + totalgradient=xNewZeroInit(num_independents); + + for(aDepIndex=0;aDepIndex(num_dependents); + if (my_rank==0) aWeightVector[aDepIndex]=1.0; + + /*initialize output gradient: */ + weightVectorTimesJac=xNew(num_independents); + + /*set the forward method function pointer: */ + #ifdef _HAVE_GSL_ + anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx; + #endif + #ifdef _HAVE_MUMPS_ + anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF; + #endif + + anEDF_for_solverx_p->dp_U=xNew(anEDF_for_solverx_p->max_m); + anEDF_for_solverx_p->dp_Z=xNew(anEDF_for_solverx_p->max_n); + + /*call driver: */ + fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac ); + + /*Add to totalgradient: */ + if(my_rank==0)for(i=0;iresults->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,totalgradient,num_independents,1,1,0.0)); + + if(VerboseAutodiff())_printf0_(" end ad core\n"); + + /* delete the allocated space for the parameters and free ressources:{{{*/ + xDelete(anEDF_for_solverx_p->dp_x); + xDelete(anEDF_for_solverx_p->dp_X); + xDelete(anEDF_for_solverx_p->dpp_X); + xDelete(anEDF_for_solverx_p->dp_y); + xDelete(anEDF_for_solverx_p->dp_Y); + xDelete(anEDF_for_solverx_p->dpp_Y); + xDelete(anEDF_for_solverx_p->dp_U); + xDelete(anEDF_for_solverx_p->dpp_U); + xDelete(anEDF_for_solverx_p->dp_Z); + xDelete(anEDF_for_solverx_p->dpp_Z); + xDelete(xp); + xDelete(totalgradient); + xDelete(axp); /*}}}*/ + #else + _error_("Should not be requesting AD drivers when an AD library is not available!"); + #endif + } +} Index: ../trunk-jpl/src/c/cores/cores.h =================================================================== --- ../trunk-jpl/src/c/cores/cores.h (revision 18895) +++ ../trunk-jpl/src/c/cores/cores.h (revision 18896) @@ -42,6 +42,7 @@ void transient_core(FemModel* femmodel); void dakota_core(FemModel* femmodel); void ad_core(FemModel* femmodel); +void adgradient_core(FemModel* femmodel); void dummy_core(FemModel* femmodel); void gia_core(FemModel* femmodel); void damage_core(FemModel* femmodel); Index: ../trunk-jpl/src/c/cores/controlad_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/controlad_core.cpp (revision 18895) +++ ../trunk-jpl/src/c/cores/controlad_core.cpp (revision 18896) @@ -215,7 +215,7 @@ _printf0_("f(x) = "<