Index: ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 23305)
+++ ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 23306)
@@ -10,7 +10,7 @@
 #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,8 +27,8 @@
 
 /*Use struct to provide arguments*/
 typedef struct{
-	FemModel  * femmodel;
-	IssmDouble* Jlist;
+	FemModel   * femmodel;
+	IssmPDouble* Jlist;
 	int         M;
 	int         N;
 	int*        i;
@@ -43,8 +43,8 @@
 	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,8 +52,8 @@
 	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,8 +77,8 @@
 	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,7 +92,7 @@
 	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,7 +111,7 @@
 	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,20 +140,36 @@
 	}
 
 	/*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,7 +186,7 @@
 	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,7 +195,7 @@
 	/*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,19 +210,31 @@
 	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,11 +248,13 @@
 
 	/*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,8 +267,8 @@
 		_printf0_("\n");
 
 		*Jlisti = (*Jlisti) +1;
-		xDelete<IssmDouble>(XU);
-		xDelete<IssmDouble>(XL);
+		xDelete<double>(XU);
+		xDelete<double>(XL);
 		return;
 	}
 
@@ -249,7 +279,7 @@
 	/*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,8 +289,8 @@
 			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,11 +303,11 @@
 
 	/*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_