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
RevLine 
[23390]1Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 23266)
4+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 23267)
5@@ -102,8 +102,8 @@
6 #if defined(_HAVE_AD_)
7
8 if(isautodiff){
9-#if _HAVE_ADOLC_
10- /*Copy some parameters from IoModel to parameters dataset: {{{*/
11+ #if defined(_HAVE_ADOLC_)
12+ /*Copy some parameters from IoModel to parameters dataset*/
13 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.obufsize",AutodiffObufsizeEnum));
14 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.cbufsize",AutodiffCbufsizeEnum));
15 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.lbufsize",AutodiffLbufsizeEnum));
16@@ -110,8 +110,14 @@
17 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.tbufsize",AutodiffTbufsizeEnum));
18 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerRatio",AutodiffGcTriggerRatioEnum));
19 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerMaxSize",AutodiffGcTriggerMaxSizeEnum));
20- /*}}}*/
21-#endif
22+
23+ #elif defined(_HAVE_CODIPACK_)
24+ parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.tapeAlloc",AutodiffTapeAllocEnum));
25+
26+ #else
27+ _error_("not supported yet");
28+ #endif
29+
30 /*retrieve driver: {{{*/
31 iomodel->FindConstant(&autodiff_driver,"md.autodiff.driver");
32 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.driver",AutodiffDriverEnum));
33Index: ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp
34===================================================================
35--- ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 23266)
36+++ ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 23267)
37@@ -10,7 +10,12 @@
38 #include "../modules/modules.h"
39 #include "../solutionsequences/solutionsequences.h"
40
41-#if defined (_HAVE_M1QN3_) && defined(_HAVE_ADOLC_)
42+#ifdef _HAVE_CODIPACK_
43+extern CoDi_global codi_global;
44+#include <sstream> // for output of the CoDiPack tape
45+#endif
46+
47+#if defined (_HAVE_M1QN3_) && defined(_HAVE_AD_)
48 /*m1qn3 prototypes*/
49 extern "C" void *ctonbe_; // DIS mode : Conversion
50 extern "C" void *ctcabe_; // DIS mode : Conversion
51@@ -34,6 +39,7 @@
52 /*m1qm3 functions*/
53 void simul_starttrace(FemModel* femmodel){/*{{{*/
54
55+ #if defined(_HAVE_ADOLC_)
56 /*Retrive ADOLC parameters*/
57 IssmDouble gcTriggerRatio;
58 IssmDouble gcTriggerMaxSize;
59@@ -56,7 +62,115 @@
60 int keepTaylors=1;
61 int my_rank=IssmComm::GetRank();
62 trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
63+
64+ #elif defined(_HAVE_CODIPACK_)
65+
66+ //fprintf(stderr, "*** Codipack IoModel::StartTrace\n");
67+ /*
68+ * FIXME codi
69+ * - ADOL-C variant uses fine grained tracing with various arguments
70+ * - ADOL-C variant sets a garbage collection parameter for its tape
71+ * -> These parameters are not read for the CoDiPack ISSM version!
72+ */
73+ auto& tape_codi = IssmDouble::getGlobalTape();
74+ tape_codi.setActive();
75+ #if _AD_TAPE_ALLOC_
76+ //alloc_profiler.Tag(StartInit, true);
77+ IssmDouble x_t(1.0), y_t(1.0);
78+ tape_codi.registerInput(y_t);
79+ int codi_allocn = 0;
80+ femmodel->parameters->FindParam(&codi_allocn,AutodiffTapeAllocEnum);
81+ for(int i = 0;i < codi_allocn;++i) {
82+ x_t = y_t * y_t;
83+ }
84+ /*
85+ std::stringstream out_s;
86+ IssmDouble::getGlobalTape().printStatistics(out_s);
87+ _printf0_("StartTrace::Tape Statistics : TapeAlloc count=[" << codi_allocn << "]\n" << out_s.str());
88+ */
89+ tape_codi.reset();
90+ //alloc_profiler.Tag(FinishInit, true);
91+ #endif
92+
93+ #else
94+ _error_("not implemented");
95+ #endif
96 }/*}}}*/
97+void simul_stoptrace(){/*{{{*/
98+
99+ #if defined(_HAVE_ADOLC_)
100+ trace_off();
101+ if(VerboseAutodiff()){ /*{{{*/
102+
103+ #ifdef _HAVE_ADOLC_
104+ int my_rank=IssmComm::GetRank();
105+ size_t tape_stats[15];
106+ tapestats(my_rank,tape_stats); //reading of tape statistics
107+ int commSize=IssmComm::GetSize();
108+ int *sstats=new int[7];
109+ sstats[0]=tape_stats[NUM_OPERATIONS];
110+ sstats[1]=tape_stats[OP_FILE_ACCESS];
111+ sstats[2]=tape_stats[NUM_LOCATIONS];
112+ sstats[3]=tape_stats[LOC_FILE_ACCESS];
113+ sstats[4]=tape_stats[NUM_VALUES];
114+ sstats[5]=tape_stats[VAL_FILE_ACCESS];
115+ sstats[6]=tape_stats[TAY_STACK_SIZE];
116+ int *rstats=NULL;
117+ if (my_rank==0) rstats=new int[commSize*7];
118+ ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
119+ if (my_rank==0) {
120+ int offset=50;
121+ int rOffset=(commSize/10)+1;
122+ _printf_(" ADOLC statistics: \n");
123+ _printf_(" "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
124+ _printf_(" "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
125+ _printf_(" "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
126+ _printf_(" operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
127+ _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
128+ for (int r=0;r<commSize;++r)
129+ _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");
130+ _printf_(" locations: entry size " << sizeof(locint) << " Bytes\n");
131+ _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
132+ for (int r=0;r<commSize;++r)
133+ _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");
134+ _printf_(" constant values: entry size " << sizeof(double) << " Bytes\n");
135+ _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
136+ for (int r=0;r<commSize;++r)
137+ _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");
138+ _printf_(" Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
139+ _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
140+ for (int r=0;r<commSize;++r)
141+ _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");
142+ delete []rstats;
143+ }
144+ delete [] sstats;
145+ #endif
146+
147+ #ifdef _HAVE_CODIPACK_
148+ #ifdef _AD_TAPE_ALLOC_
149+ //_printf_("Allocation time P(" << my_rank << "): " << alloc_profiler.DeltaTime(StartInit, FinishInit) << "\n");
150+ #endif
151+ std::stringstream out_s;
152+ IssmDouble::getGlobalTape().printStatistics(out_s);
153+ _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str());
154+ #endif
155+ } /*}}}*/
156+
157+ #elif defined(_HAVE_CODIPACK_)
158+ auto& tape_codi = IssmDouble::getGlobalTape();
159+ tape_codi.setPassive();
160+ if(VerboseAutodiff()){
161+ int my_rank=IssmComm::GetRank();
162+ if(my_rank == 0) {
163+ // FIXME codi "just because" for now
164+ tape_codi.printStatistics(std::cout);
165+ codi_global.print(std::cout);
166+ }
167+ }
168+ #else
169+ _error_("not implemented");
170+ #endif
171+}/*}}}*/
172 void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
173
174 /*Get rank*/
175@@ -72,9 +186,8 @@
176 int N_add = 0;
177 int* control_enum = NULL;
178
179- if (solution_type == TransientSolutionEnum){
180- femmodel = input_struct->femmodel->copy();
181- }
182+ /*In transient, we need to make sure we do not modify femmodel at each iteration, make a copy*/
183+ if(solution_type == TransientSolutionEnum) femmodel = input_struct->femmodel->copy();
184
185 IssmPDouble* Jlist = input_struct->Jlist;
186 int JlistM = input_struct->M;
187@@ -116,11 +229,25 @@
188 #else
189 IssmDouble* aX=xNew<IssmDouble>(intn);
190 #endif
191+
192+ #if defined(_HAVE_ADOLC_)
193 if(my_rank==0){
194 for(int i=0;i<intn;i++){
195 aX[i]<<=X[i];
196 }
197 }
198+ #elif defined(_HAVE_CODIPACK_)
199+ auto& tape_codi = IssmDouble::getGlobalTape();
200+ if(my_rank==0){
201+ for (int i=0;i<intn;i++) {
202+ aX[i]=X[i];
203+ tape_codi.registerInput(aX[i]);
204+ codi_global.input_indices.push_back(aX[i].getGradientData());
205+ }
206+ }
207+ #else
208+ _error_("not suppoted");
209+ #endif
210
211 ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
212 SetControlInputsFromVectorx(femmodel,aX);
213@@ -154,73 +281,28 @@
214 DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
215 if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
216 if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
217- if (my_rank==0) {
218+ if(my_rank==0) {
219+
220+ #if defined(_HAVE_CODIPACK_)
221+ tape_codi.registerOutput(output_value);
222+ dependents[i] = output_value.getValue();
223+ codi_global.output_indices.push_back(output_value.getGradientData());
224+
225+ #elif defined(_HAVE_ADOLC_)
226 output_value>>=dependents[i];
227+
228+ #else
229+ _error_("not suppoted");
230+ #endif
231 J+=output_value;
232 }
233 }
234
235 /*Turning off trace tape*/
236- trace_off();
237+ simul_stoptrace();
238 //time_t now = time(NULL);
239 //if(my_rank==0) _printf_("\nTIME: "<<now<<"\n");
240
241- /*Print tape statistics so that user can kill this run if something is off already:*/
242- if(VerboseAutodiff()){ /*{{{*/
243-
244- #ifdef _HAVE_ADOLC_
245- size_t tape_stats[15];
246- tapestats(my_rank,tape_stats); //reading of tape statistics
247- int commSize=IssmComm::GetSize();
248- int *sstats=new int[7];
249- sstats[0]=tape_stats[NUM_OPERATIONS];
250- sstats[1]=tape_stats[OP_FILE_ACCESS];
251- sstats[2]=tape_stats[NUM_LOCATIONS];
252- sstats[3]=tape_stats[LOC_FILE_ACCESS];
253- sstats[4]=tape_stats[NUM_VALUES];
254- sstats[5]=tape_stats[VAL_FILE_ACCESS];
255- sstats[6]=tape_stats[TAY_STACK_SIZE];
256- int *rstats=NULL;
257- if (my_rank==0) rstats=new int[commSize*7];
258- ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
259- if (my_rank==0) {
260- int offset=50;
261- int rOffset=(commSize/10)+1;
262- _printf_(" ADOLC statistics: \n");
263- _printf_(" "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
264- _printf_(" "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
265- _printf_(" "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
266- _printf_(" operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
267- _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
268- for (int r=0;r<commSize;++r)
269- _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");
270- _printf_(" locations: entry size " << sizeof(locint) << " Bytes\n");
271- _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
272- for (int r=0;r<commSize;++r)
273- _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");
274- _printf_(" constant values: entry size " << sizeof(double) << " Bytes\n");
275- _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
276- for (int r=0;r<commSize;++r)
277- _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");
278- _printf_(" Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
279- _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
280- for (int r=0;r<commSize;++r)
281- _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");
282- delete []rstats;
283- }
284- delete [] sstats;
285- #endif
286-
287- #ifdef _HAVE_CODIPACK_
288- #ifdef _AD_TAPE_ALLOC_
289- //_printf_("Allocation time P(" << my_rank << "): " << alloc_profiler.DeltaTime(StartInit, FinishInit) << "\n");
290- #endif
291- std::stringstream out_s;
292- IssmDouble::getGlobalTape().printStatistics(out_s);
293- _printf0_("CoDiPack Profiling::Tape Statistics :\n" << out_s.str());
294- #endif
295- } /*}}}*/
296-
297 /*diverse: */
298 int dummy;
299 int num_independents=0;
300@@ -246,6 +328,9 @@
301 num_independents = 0;
302 }
303
304+ #if defined(_HAVE_ADOLC_)
305+ /*Get gradient for ADOLC {{{*/
306+
307 /*get the EDF pointer:*/
308 ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
309
310@@ -289,7 +374,42 @@
311 xDelete(weightVectorTimesJac);
312 xDelete(aWeightVector);
313 }
314+ /*}}}*/
315+ #elif defined(_HAVE_CODIPACK_)
316+ /*Get gradient for CoDiPack{{{*/
317+ if(VerboseAutodiff())_printf0_(" CoDiPack fos_reverse\n");
318+ int aDepIndex=0; /*FIXME: do we really need this?*/
319
320+ /*retrieve direction index: */
321+ femmodel->parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum);
322+ if (my_rank==0) {
323+ if (aDepIndex<0 || aDepIndex>=num_dependents || codi_global.output_indices.size() <= aDepIndex){
324+ _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]");
325+ }
326+ tape_codi.setGradient(codi_global.output_indices[aDepIndex], 1.0);
327+ }
328+ tape_codi.evaluate();
329+
330+ weightVectorTimesJac=xNew<double>(num_independents);
331+ /*call driver: */
332+ auto in_size = codi_global.input_indices.size();
333+ for(size_t i = 0; i < in_size; ++i) {
334+ weightVectorTimesJac[i] = tape_codi.getGradient(codi_global.input_indices[i]);
335+ }
336+
337+ /*Add to totalgradient: */
338+ totalgradient=xNewZeroInit<IssmPDouble>(num_independents_old);
339+ if(my_rank==0) for(int i=0;i<num_independents;i++) {
340+ totalgradient[i]+=weightVectorTimesJac[i];
341+ }
342+
343+ /*free resources :*/
344+ xDelete(weightVectorTimesJac);
345+ /*}}}*/
346+ #else
347+ _error_("not suppoted");
348+ #endif
349+
350 /*Broadcast gradient to other ranks*/
351 ISSM_MPI_Bcast(totalgradient,num_independents_old,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
352 /*Check size of Jlist to avoid crashes*/
353@@ -550,5 +670,5 @@
354 }/*}}}*/
355
356 #else
357-void controladm1qn3_core(FemModel* femmodel){_error_("M1QN3 or ADOLC not installed");}
358+void controladm1qn3_core(FemModel* femmodel){_error_("M1QN3 or ADOLC/CoDiPack not installed");}
359 #endif //_HAVE_M1QN3_
Note: See TracBrowser for help on using the repository browser.