Index: /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 21448)
+++ /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 21449)
@@ -26,4 +26,13 @@
 void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs);
 
+/*Use struct to provide arguments*/
+typedef struct{
+	FemModel  * femmodel;
+	IssmDouble* Jlist;
+	int         M;
+	int         N;
+	int*        i;
+} m1qn3_struct;
+
 void controlm1qn3_core(FemModel* femmodel){
 
@@ -32,5 +41,5 @@
 	double       f,dxmin,gttol; 
 	int          maxsteps,maxiter;
-	int          intn,numberofvertices,num_controls,solution_type;
+	int          intn,numberofvertices,num_controls,num_cost_functions,solution_type;
 	IssmDouble  *scaling_factors = NULL;
 	IssmDouble  *X  = NULL;
@@ -40,4 +49,5 @@
 	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+	femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum);
 	femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum);
 	femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
@@ -86,5 +96,5 @@
 	}
 
-	/*Allocate m1qn3 working arrays (see doc)*/
+	/*Allocate m1qn3 working arrays (see documentation)*/
 	long      m   = 100;
 	long      ndz = 4*n+m*(2*n+1);
@@ -96,7 +106,15 @@
 	_printf0_("____________________________________________________________________\n");
 
+	/*Prepare structure for m1qn3*/
+	m1qn3_struct mystruct;
+	mystruct.femmodel = femmodel;
+	mystruct.M        = maxiter;
+	mystruct.N        = num_cost_functions+1;
+	mystruct.Jlist    = xNewZeroInit<IssmDouble>(mystruct.M*mystruct.N);
+	mystruct.i        = xNewZeroInit<int>(1);
+
 	/*Initialize Gradient and cost function of M1QN3*/
-	indic = 4; //adjoint and gradient required
-	simul(&indic,&n,X,&f,G,izs,rzs,(void*)femmodel);
+	indic = 4; /*gradient required*/
+	simul(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
 
 	/*Estimation of the expected decrease in f during the first iteration*/
@@ -107,5 +125,5 @@
 				&n,X,&f,G,&dxmin,&df1,
 				&gttol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz,
-				&reverse,&indic,izs,rzs,(void*)femmodel);
+				&reverse,&indic,izs,rzs,(void*)&mystruct);
 
 	switch(int(omode)){
@@ -137,5 +155,5 @@
 	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
 	femmodel->OutputControlsx(&femmodel->results);
-	femmodel->results->AddObject(new GenericExternalResult<double>(femmodel->results->Size()+1,JEnum,f));
+	femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,mystruct.M,mystruct.N,0,0));
 
 	/*Finalize*/
@@ -153,4 +171,6 @@
 	xDelete<double>(XL);
 	xDelete<double>(scaling_factors);
+	xDelete<double>(mystruct.Jlist);
+	xDelete<int>(mystruct.i);
 }
 
@@ -158,17 +178,22 @@
 void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){
 
-	/*Recover Femmodel*/
-	int         solution_type;
-	FemModel   *femmodel  = (FemModel*)dzs;
-
-	/*Recover number of cost functions responses*/
-	int num_responses,num_controls,numberofvertices;
+	/*Recover Arguments*/
+	m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
+	FemModel     *femmodel     = input_struct->femmodel;
+	IssmDouble   *Jlist        = input_struct->Jlist;
+	int           JlistM       = input_struct->M;
+	int           JlistN       = input_struct->N;
+	int          *Jlisti       = input_struct->i;
+
+	/*Recover some parameters*/
+	int num_responses,num_controls,numberofvertices,solution_type;
 	IssmDouble* scaling_factors = NULL;
 	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
 	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
 	femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 	numberofvertices=femmodel->vertices->NumberOfVertices();
 
-	/*Constrain input vector*/
+	/*Constrain input vector and update controls*/
 	IssmDouble  *XL = NULL;
 	IssmDouble  *XU = NULL;
@@ -183,10 +208,5 @@
 		}
 	}
-
-	/*Update control input*/
 	SetControlInputsFromVectorx(femmodel,X);
-
-	/*Recover some parameters*/
-	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 
 	/*Compute solution and adjoint*/
@@ -196,8 +216,17 @@
 	solutioncore(femmodel);
 
+	/*Check size of Jlist to avoid crashes*/
+	_assert_((*Jlisti)<JlistM);
+	_assert_(JlistN==num_responses+1);
+
 	/*Compute objective function*/
-	IssmDouble* Jlist = NULL;
-	femmodel->CostFunctionx(pf,&Jlist,NULL);
+	IssmDouble* Jtemp = NULL;
+	femmodel->CostFunctionx(pf,&Jtemp,NULL);
 	_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];
+	Jlist[(*Jlisti)*JlistN+num_responses] = *pf;
+	xDelete<IssmDouble>(Jtemp);
 
 	if(*indic==0){
@@ -206,8 +235,8 @@
 		/*Retrieve objective functions independently*/
 		_printf0_("            N/A |\n");
-		for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
+		for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
 		_printf0_("\n");
 
-		xDelete<IssmDouble>(Jlist);
+		*Jlisti = (*Jlisti) +1;
 		xDelete<IssmDouble>(XU);
 		xDelete<IssmDouble>(XL);
@@ -241,9 +270,9 @@
 	/*Print info*/
 	_printf0_("       "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
-	for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
+	for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
 	_printf0_("\n");
 
 	/*Clean-up and return*/
-	xDelete<IssmDouble>(Jlist);
+	*Jlisti = (*Jlisti) +1;
 	xDelete<IssmDouble>(XU);
 	xDelete<IssmDouble>(XL);