Index: /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 17898)
+++ /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 17899)
@@ -3,4 +3,5 @@
  */ 
 
+#include <config.h>
 #include "./cores.h"
 #include "../toolkits/toolkits.h"
@@ -15,32 +16,40 @@
 extern "C" void *ctcabe_; // DIS mode : Conversion
 extern "C" void *euclid_; // Scalar product
-typedef void (*SimulFunc) (long*, long *, double [], double *, double [], long [], float [], double []);
-extern "C" void m1qn3_ (void f(long*, long *, double [], double *, double [], long [], float [], double []),
+typedef void (*SimulFunc) (long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs);
+extern "C" void m1qn3_ (void f(long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs),
 			void **, void **, void **,
 			long *, double [], double *, double [], double*, double *,
 			double *, char [], long *, long *, long *, long *, long *, long *, long [], double [], long *,
-			long *, long *, long [], float [], double []
-			);
+			long *, long *, long [], float [],void* );
 
-void simul(long* indic,long* n,double x[2],double* f,double g[2],long izs[1],float rzs[1],double dzs[1]){
-	/*minimization of x^2+cy^4 where c>0 is a parameter*/
-	double c=dzs[0];
-	*f=x[0]*x[0]+c*x[1]*x[1]*x[1]*x[1];
-	g[0]=2.*x[0];
-	g[1]=4.*c*x[1]*x[1]*x[1];
-}
+/*Cost function prototype*/
+void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs);
 
 void controlm1qn3_core(FemModel* femmodel){
 
 	/*Intermediaries*/
-	long   omode; /*m1qn3 output flag*/
-	double f;     /*cost function value*/
+	long         omode;
+	double       f; 
+	int          nsteps,maxiter;
+	int          intn,num_controls,solution_type;
+	IssmDouble  *X  = NULL;
+	IssmDouble  *G  = NULL;
+	IssmDouble  *XL = NULL;
+	IssmDouble  *XU = NULL;
 
+	/*Recover some parameters*/
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+	femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum);
+	femmodel->parameters->SetParam(false,SaveResultsEnum);
+	maxiter=nsteps*10;
+
+	/*Initialize M1QN3 parameters*/
 	if(VerboseControl())_printf0_("   Initialize M1QN3 parameters\n");
 	SimulFunc costfuncion  = &simul;    /*Cost function address*/
 	void**    prosca       = &euclid_;  /*Dot product function (euclid is the default)*/
 	char      normtype[3];              /*Norm type: dfn = scalar product defined by prosca*/ strcpy(normtype, "dfn");
-	double    dzs[1];                   /*Arrays used by m1qn3 subroutines*/
 	long      izs[5];                   /*Arrays used by m1qn3 subroutines*/
+	long      iz[5];                    /*Integer m1qn3 working array of size 5*/
 	float     rzs[1];                   /*Arrays used by m1qn3 subroutines*/
 	long      impres       = 0;         /*verbosity level*/
@@ -48,19 +57,22 @@
 	long      indic        = 4;         /*compute f and g*/
 	long      reverse      = 0;         /*reverse or direct mode*/
-	long      iz[5];                    /*Integer m1qn3 working array of size 5*/
 	long      io           = 6;         /*Channel number for the output*/
 
 	/*Optimization criterions*/
-	double    dxmin        = 1.e-10;    /*Resolution for the solution x*/
-	double    epsrel       = 1.e-5;     /*Gradient stopping criterion in ]0 1[ -> |gk|/|g1| < epsrel*/
-	long      niter        = 200;       /*Maximum number of iterations*/
-	long      nsim         = 200;       /*Maximum number of function calls*/
+	double    dxmin        = 1.e-1;    /*Resolution for the solution x*/
+	double    epsrel       = 1.e-4;     /*Gradient stopping criterion in ]0 1[ -> |gk|/|g1| < epsrel*/
+	long      niter        = long(nsteps); /*Maximum number of iterations*/
+	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");
+	X = Xpetsc->ToMPISerial();
+	Xpetsc->GetSize(&intn);
+	delete Xpetsc;
 
 	/*Get problem dimension and initialize gradient and initial guess*/
-	long    n = 2;
-	double* g = xNew<double>(n);
-	double* x = xNew<double>(n);
-	for(int i=0;i<n;i++) x[i]=5.;
-	dzs[0] = 10; //c = 10 function parameter
+	long n = long(intn);
+	G = xNew<double>(n);
 
 	/*Allocate m1qn3 working arrays (see doc)*/
@@ -70,11 +82,16 @@
 
 	if(VerboseControl())_printf0_("   Computing initial solution\n");
-	simul(&indic,&n,x,&f,g,izs,rzs,&dzs[0]);
+	_printf0_("\n");
+	_printf0_("Cost function f(x)   |  List of contributions\n");
+	_printf0_("_____________________________________________\n");
+	indic = 0; //no adjoint required
+	simul(&indic,&n,X,&f,G,izs,rzs,(void*)femmodel);
 	double f1=f;
 
+	indic = 4; //adjoint and gradient required
 	m1qn3_(costfuncion,prosca,&ctonbe_,&ctcabe_,
-				&n,x,&f,g,&dxmin,&f1,
+				&n,X,&f,G,&dxmin,&f1,
 				&epsrel,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz,
-				&reverse,&indic,izs,rzs,dzs);
+				&reverse,&indic,izs,rzs,(void*)femmodel);
 
 	switch(int(omode)){
@@ -90,15 +107,78 @@
 	}
 
-	_printf0_(" == Final cost function = "<< f <<"\n");
-	//_printf0_(" == Final x = ["<<x[0]<<" "<<x[1]<<"]\n");
+	/*Get solution*/
+	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
+	femmodel->OutputControlsx(&femmodel->results);
+	femmodel->results->AddObject(new GenericExternalResult<double>(JEnum,f,1,0));
+
+	/*Finalize*/
+	if(VerboseControl()) _printf0_("   preparing final solution\n");
+	femmodel->parameters->SetParam(true,SaveResultsEnum);
+	void (*solutioncore)(FemModel*)=NULL;
+	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
+	solutioncore(femmodel);
 
 	/*Clean-up and return*/
-	xDelete<double>(g);
-	xDelete<double>(x);
+	xDelete<double>(G);
+	xDelete<double>(X);
 	xDelete<double>(dz);
 }
+
+/*Cost function definition*/
+void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){
+
+	/*Recover Femmodel*/
+	double      f;
+	int         intn,solution_type;
+	FemModel   *femmodel  = (FemModel*)dzs;
+
+	/*Recover responses*/
+	int         num_responses;
+	int        *responses = NULL;
+	femmodel->parameters->FindParam(&responses,&num_responses,InversionCostFunctionsEnum);
+
+	/*Update control input*/
+	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+
+	/*Recover some parameters*/
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+
+	/*Compute solution and adjoint*/
+	void (*solutioncore)(FemModel*)=NULL;
+	void (*adjointcore)(FemModel*)=NULL;
+	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
+	solutioncore(femmodel);
+
+	/*Compute objective function*/
+	femmodel->CostFunctionx(pf);
+	_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
+
+	/*Retrieve objective functions independently*/
+	for(int i=0;i<num_responses;i++){
+		femmodel->Responsex(&f,EnumToStringx(responses[i]));
+		_printf0_(" "<<setw(12)<<setprecision(7)<<f);
+	}
+	_printf0_("\n");
+
+	if(indic==0) return; /*dry run, no gradient required*/
+
+	/*Compute Adjoint*/
+	AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
+	adjointcore(femmodel);
+
+	/*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];
+
+	/*Clean-up and return*/
+	xDelete<int>(responses);
+	xDelete<IssmDouble>(G2);
+}
+
 #else
 void controlm1qn3_core(FemModel* femmodel){
 	_error_("M1QN3 not installed");
 }
-#endif //_HAVE_TAO_ 
+#endif //_HAVE_M1QN3_
Index: /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp	(revision 17898)
+++ /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp	(revision 17899)
@@ -66,2 +66,19 @@
 	xDelete<int>(control_type);
 }
+void Gradjx(IssmDouble** pgradient,IssmDouble** pnorm_list, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){
+
+	/*output: */
+	IssmDouble* gradient=NULL;
+
+	/*intermediary: */
+	Vector<IssmDouble>* vec_gradient=NULL;
+
+	Gradjx(&vec_gradient,pnorm_list,elements,nodes, vertices,loads,materials,parameters);
+	gradient=vec_gradient->ToMPISerial();
+
+	/*Free ressources:*/
+	delete vec_gradient;
+
+	/*Assign output pointers:*/
+	*pgradient=gradient;
+}
Index: /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h	(revision 17898)
+++ /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h	(revision 17899)
@@ -10,4 +10,5 @@
 /* local prototypes: */
 void Gradjx(Vector<IssmDouble>** pgrad_g,IssmDouble** pgrad_norm,Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters);
+void Gradjx(IssmDouble** pgrad_g,IssmDouble** pgrad_norm,Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters);
 
 #endif  /* _GRADJX_H */