source: issm/oecreview/Archive/23185-23389/ISSM-23305-23306.diff@ 23390

../trunk-jpl/src/c/cores/controlm1qn3_core.cpp

@@ -10 +10 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
 
-#if defined (_HAVE_M1QN3_) && !defined(_HAVE_AD_)
+#if defined (_HAVE_M1QN3_)
 /*m1qn3 prototypes {{{*/
 extern "C" void *ctonbe_; // DIS mode : Conversion
 extern "C" void *ctcabe_; // DIS mode : Conversion
@@ -27 +27 @@
 
 /*Use struct to provide arguments*/
 typedef struct{
-        FemModel  * femmodel;
-        IssmDouble* Jlist;
+        FemModel   * femmodel;
+        IssmPDouble* Jlist;
         int         M;
         int         N;
         int*        i;
@@ -43 +43 @@
         int          maxsteps,maxiter;
         int          intn,numberofvertices,num_controls,num_cost_functions,solution_type;
         IssmDouble  *scaling_factors = NULL;
-        IssmDouble  *X  = NULL;
-        IssmDouble  *G  = NULL;
+        double      *X  = NULL;
+        double      *G  = NULL;
 
         /*Recover some parameters*/
         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
@@ -52 +52 @@
         femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum);
         femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum);
         femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
-        femmodel->parameters->FindParam(&dxmin,InversionDxminEnum);
-        femmodel->parameters->FindParam(&gttol,InversionGttolEnum);
+        femmodel->parameters->FindParamAndMakePassive(&dxmin,InversionDxminEnum);
+        femmodel->parameters->FindParamAndMakePassive(&gttol,InversionGttolEnum);
         femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
         femmodel->parameters->SetParam(false,SaveResultsEnum);
         numberofvertices=femmodel->vertices->NumberOfVertices();
@@ -77 +77 @@
         long nsim  = long(maxiter);/*Maximum number of function calls*/
 
         /*Get initial guess*/
-        Vector<IssmDouble> *Xpetsc = NULL;
-        GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
+        Vector<double> *Xpetsc = NULL;
+        GetPassiveVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
         X = Xpetsc->ToMPISerial();
         Xpetsc->GetSize(&intn);
         delete Xpetsc;
@@ -92 +92 @@
         for(int c=0;c<num_controls;c++){
                 for(int i=0;i<numberofvertices;i++){
                         int index = numberofvertices*c+i;
-                        X[index] = X[index]/scaling_factors[c];
+                        X[index] = X[index]/reCast<double>(scaling_factors[c]);
                 }
         }
 
@@ -111 +111 @@
         mystruct.femmodel = femmodel;
         mystruct.M        = maxiter;
         mystruct.N        = num_cost_functions+1;
-        mystruct.Jlist    = xNewZeroInit<IssmDouble>(mystruct.M*mystruct.N);
+        mystruct.Jlist    = xNewZeroInit<IssmPDouble>(mystruct.M*mystruct.N);
         mystruct.i        = xNewZeroInit<int>(1);
 
         /*Initialize Gradient and cost function of M1QN3*/
@@ -140 +140 @@
         }
 
         /*Constrain solution vector*/
-        IssmDouble  *XL = NULL;
-        IssmDouble  *XU = NULL;
-        GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
-        GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+        double  *XL = NULL;
+        double  *XU = NULL;
+        GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
+        GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
         for(int c=0;c<num_controls;c++){
                 for(int i=0;i<numberofvertices;i++){
                         int index = numberofvertices*c+i;
-                        X[index] = X[index]*scaling_factors[c];
+                        X[index] = X[index]*reCast<double>(scaling_factors[c]);
                         if(X[index]>XU[index]) X[index]=XU[index];
                         if(X[index]<XL[index]) X[index]=XL[index];
                 }
         }
+
+        /*Set X as our new control (need to recast)*/
+        #ifdef _HAVE_AD_
+        IssmDouble* aX=xNew<IssmDouble>(intn);
+        IssmDouble* aG=xNew<IssmDouble>(intn);
+        for(int i=0;i<intn;i++) {
+                aX[i] = reCast<IssmDouble>(X[i]); 
+                aG[i] = reCast<IssmDouble>(G[i]);
+        }
+        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
+        SetControlInputsFromVectorx(femmodel,aX);
+        xDelete(aX);
+        xDelete(aG);
+        #else
         SetControlInputsFromVectorx(femmodel,X);
         ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
+        #endif
+
         femmodel->OutputControlsx(&femmodel->results);
         femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0));
 
@@ -170 +186 @@
         xDelete<double>(dz);
         xDelete<double>(XU);
         xDelete<double>(XL);
-        xDelete<double>(scaling_factors);
+        xDelete<IssmDouble>(scaling_factors);
         xDelete<double>(mystruct.Jlist);
         xDelete<int>(mystruct.i);
 }/*}}}*/
@@ -179 +195 @@
         /*Recover Arguments*/
         m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
         FemModel     *femmodel     = input_struct->femmodel;
-        IssmDouble   *Jlist        = input_struct->Jlist;
+        IssmPDouble  *Jlist        = input_struct->Jlist;
         int           JlistM       = input_struct->M;
         int           JlistN       = input_struct->N;
         int          *Jlisti       = input_struct->i;
@@ -194 +210 @@
         numberofvertices=femmodel->vertices->NumberOfVertices();
 
         /*Constrain input vector and update controls*/
-        IssmDouble  *XL = NULL;
-        IssmDouble  *XU = NULL;
+        double *XL = NULL;
+        double *XU = NULL;
+        #ifdef _HAVE_AD_
+        GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
+        GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+        #else
         GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
         GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+        #endif
         for(int c=0;c<num_controls;c++){
                 for(int i=0;i<numberofvertices;i++){
                         int index = numberofvertices*c+i;
-                        X[index] = X[index]*scaling_factors[c];
+                        X[index] = X[index]*reCast<double>(scaling_factors[c]);
                         if(X[index]>XU[index]) X[index]=XU[index];
                         if(X[index]<XL[index]) X[index]=XL[index];
                 }
         }
+        #ifdef _HAVE_AD_
+        IssmDouble* aX=xNew<IssmDouble>(*n);
+        for(int i=0;i<*n;i++) aX[i] = reCast<IssmDouble>(X[i]); 
+        SetControlInputsFromVectorx(femmodel,aX);
+        xDelete(aX);
+        #else
         SetControlInputsFromVectorx(femmodel,X);
+        #endif
 
         /*Compute solution and adjoint*/
         void (*solutioncore)(FemModel*)=NULL;
@@ -220 +248 @@
 
         /*Compute objective function*/
         IssmDouble* Jtemp = NULL;
-        femmodel->CostFunctionx(pf,&Jtemp,NULL);
+        IssmDouble  J;
+        femmodel->CostFunctionx(&J,&Jtemp,NULL);
+        *pf = reCast<double>(J);
         _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
 
         /*Record cost function values and delete Jtemp*/
-        for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = Jtemp[i];
+        for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = reCast<double>(Jtemp[i]);
         Jlist[(*Jlisti)*JlistN+num_responses] = *pf;
         xDelete<IssmDouble>(Jtemp);
 
@@ -237 +267 @@
                 _printf0_("\n");
 
                 *Jlisti = (*Jlisti) +1;
-                xDelete<IssmDouble>(XU);
-                xDelete<IssmDouble>(XL);
+                xDelete<double>(XU);
+                xDelete<double>(XL);
                 return;
         }
 
@@ -249 +279 @@
         /*Compute gradient*/
         IssmDouble* G2 = NULL;
         Gradjx(&G2,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-        for(long i=0;i<*n;i++) G[i] = -G2[i];
+        for(long i=0;i<*n;i++) G[i] = -reCast<double>(G2[i]);
         xDelete<IssmDouble>(G2);
 
         /*Constrain Gradient*/
@@ -259 +289 @@
                         int index = numberofvertices*c+i;
                         if(X[index]>=XU[index]) G[index]=0.;
                         if(X[index]<=XL[index]) G[index]=0.;
-                        G[index] = G[index]*scaling_factors[c];
-                        X[index] = X[index]/scaling_factors[c];
+                        G[index] = G[index]*reCast<double>(scaling_factors[c]);
+                        X[index] = X[index]/reCast<double>(scaling_factors[c]);
                         Gnorm += G[index]*G[index];
                 }
         }
@@ -273 +303 @@
 
         /*Clean-up and return*/
         *Jlisti = (*Jlisti) +1;
-        xDelete<IssmDouble>(XU);
-        xDelete<IssmDouble>(XL);
+        xDelete<double>(XU);
+        xDelete<double>(XL);
         xDelete<IssmDouble>(scaling_factors);
 }/*}}}*/
 
 #else
-void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed (or you might have turned AD on)");}
+void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");}
 #endif //_HAVE_M1QN3_