Index: /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp	(revision 22828)
+++ /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp	(revision 22829)
@@ -30,4 +30,7 @@
 	int          N;
 	int*         i;
+	IssmPDouble* X_best;
+	IssmPDouble* G_best;
+	IssmPDouble* J_best;
 } m1qn3_struct;
 
@@ -77,8 +80,11 @@
 	}
 
-	IssmPDouble  *Jlist        = input_struct->Jlist;
+	IssmPDouble*  Jlist        = input_struct->Jlist;
 	int           JlistM       = input_struct->M;
 	int           JlistN       = input_struct->N;
-	int          *Jlisti       = input_struct->i;
+	int*          Jlisti       = input_struct->i;
+	IssmPDouble*  J_best       = input_struct->J_best;
+	IssmPDouble*  X_best			= input_struct->X_best;
+	IssmPDouble*  G_best			= input_struct->G_best;
 	int           intn         = (int)*n;
 
@@ -292,4 +298,19 @@
 	Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J);
 
+	
+	if(*J_best<0 || J<*J_best){
+		*J_best = reCast<IssmPDouble>(J);
+		for(int i=0;i<intn;i++){
+			X_best[i] = reCast<IssmPDouble>(X[i]);
+			G_best[i] = reCast<IssmPDouble>(G[i]);
+		}
+	}
+
+/*
+	IssmDouble* test = xNew<IssmDouble>(intn);
+	femmodel->parameters->FindParam(&test,NULL,InversionXBestEnum);
+	for(int i=0;i<10;i++)_printf_("X "<<test[i]<<"\n");
+	xDelete<IssmDouble>(intn);
+*/
 	if(*indic==0){
 		/*dry run, no gradient required*/
@@ -413,5 +434,5 @@
 	long      ndz = 4*n+m*(2*n+1);
 	double*   dz  = xNew<double>(ndz);
-
+	IssmDouble J_best = -10.;
 	if(VerboseControl())_printf0_("   Computing initial solution\n");
 	_printf0_("\n");
@@ -426,4 +447,7 @@
 	mystruct.Jlist    = xNewZeroInit<IssmPDouble>(mystruct.M*mystruct.N);
 	mystruct.i        = xNewZeroInit<int>(1);
+	mystruct.J_best   = xNewZeroInit<IssmPDouble>(1);
+	mystruct.X_best	= xNewZeroInit<IssmPDouble>(intn);
+	mystruct.G_best   = xNewZeroInit<IssmPDouble>(intn);
 	/*Initialize Gradient and cost function of M1QN3*/
 	indic = 4; /*gradient required*/
@@ -465,19 +489,29 @@
 		N_add +=N[c];
 	}
-		
+	
+	IssmDouble* aX_best = NULL;
+	IssmDouble* aG_best = NULL;
+
+	femmodel->parameters->FindParam(&aX_best,NULL,InversionXBestEnum);
+	femmodel->parameters->FindParam(&aG_best,NULL,InversionGBestEnum);
+	
 	/*Set X as our new control*/
 	IssmDouble* aX=xNew<IssmDouble>(intn);
+	IssmDouble* aG=xNew<IssmDouble>(intn);
+	double* X_best = xNew<double>(intn);
+	double* G_best = xNew<double>(intn);
+
 	for(int i=0;i<intn;i++) {
 		aX[i] = reCast<IssmDouble>(X[i]); 
-		}
+		aG[i] = reCast<IssmDouble>(G[i]);
+		X_best[i] = reCast<double>(aX_best[i]);	
+		G_best[i] = reCast<double>(aG_best[i]);	
+		}
+
+	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
 	SetControlInputsFromVectorx(femmodel,aX);
+	
 	xDelete(aX);
 
-	/*Set final gradient in inputs*/
-	IssmDouble* aG=xNew<IssmDouble>(intn);
-	for(int i=0;i<intn;i++) {
-		aG[i] = reCast<IssmDouble>(G[i]);
-	}
-	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
 
 	if (solution_type == TransientSolutionEnum){
@@ -493,12 +527,19 @@
 			GenericExternalResult<IssmPDouble*>* X_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,control_enum[i],&X[offset],N[i],numberofvertices,1,0.);
 
+			GenericExternalResult<IssmPDouble*>* Gbest_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,Outputdefinition90Enum+i,&G_best[offset],N[i],numberofvertices,1,0.);
+			GenericExternalResult<IssmPDouble*>* Xbest_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,Outputdefinition80Enum+i,&X_best[offset],N[i],numberofvertices,1,0.);
+
 			/*transpose for consistency with MATLAB's formating*/
 			G_output->Transpose();
 			X_output->Transpose();
+			Gbest_output->Transpose();
+			Xbest_output->Transpose();
 
 			/*Add to results*/
 			femmodel->results->AddObject(G_output);
 			femmodel->results->AddObject(X_output);
-
+			femmodel->results->AddObject(Gbest_output);
+			femmodel->results->AddObject(Xbest_output);
+			
 			offset += N[i]*numberofvertices;
 		}
@@ -524,4 +565,6 @@
 	xDelete<double>(G);
 	xDelete<double>(X);
+	xDelete<double>(X_best);
+	xDelete<double>(G_best);
 	xDelete<double>(dz);
 	xDelete<double>(XU);