Index: issm/trunk-jpl/src/c/classes/Misfit.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Misfit.cpp	(revision 22437)
+++ issm/trunk-jpl/src/c/classes/Misfit.cpp	(revision 22438)
@@ -12,5 +12,4 @@
 
 #include "./classes.h"
-#include "./FemModel.h"
 #include "./ExternalResults/ExternalResult.h"
 #include "./ExternalResults/Results.h"
@@ -18,4 +17,5 @@
 #include "./Elements/Element.h"
 #include "./Elements/Elements.h"
+#include "./FemModel.h"
 #include "../modules/SurfaceAreax/SurfaceAreax.h"
 #include "../classes/Params/Parameters.h"
Index: issm/trunk-jpl/src/c/classes/Misfit.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Misfit.h	(revision 22437)
+++ issm/trunk-jpl/src/c/classes/Misfit.h	(revision 22438)
@@ -8,5 +8,12 @@
 /*Headers:*/
 #include "./Definition.h"
+#include "../datastructures/datastructures.h"
+#include "./Elements/Element.h"
+#include "./Elements/Elements.h"
 #include "./FemModel.h"
+#include "../modules/SurfaceAreax/SurfaceAreax.h"
+#include "../classes/Params/Parameters.h"
+#include "../classes/Inputs/Input.h"
+#include "../classes/gauss/Gauss.h"
 
 IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,int output_enum);
Index: issm/trunk-jpl/src/c/classes/Nodalvalue.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Nodalvalue.h	(revision 22437)
+++ issm/trunk-jpl/src/c/classes/Nodalvalue.h	(revision 22438)
@@ -9,5 +9,12 @@
 /*{{{*/
 #include "./Definition.h"
+#include "../datastructures/datastructures.h"
+#include "./Elements/Element.h"
+#include "./Elements/Elements.h"
 #include "./FemModel.h"
+#include "../modules/SurfaceAreax/SurfaceAreax.h"
+#include "../classes/Params/Parameters.h"
+#include "../classes/Inputs/Input.h"
+#include "../classes/gauss/Gauss.h"
 /*}}}*/
 
Index: issm/trunk-jpl/src/c/classes/Numberedcostfunction.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Numberedcostfunction.h	(revision 22437)
+++ issm/trunk-jpl/src/c/classes/Numberedcostfunction.h	(revision 22438)
@@ -8,5 +8,19 @@
 /*Headers:*/
 #include "./Definition.h"
+#include "../datastructures/datastructures.h"
+#include "./Elements/Element.h"
+#include "./Elements/Elements.h"
 #include "./FemModel.h"
+#include "./ExternalResults/ExternalResult.h"
+#include "./ExternalResults/Results.h"
+#include "../modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h"
+#include "../modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h"
+#include "../modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h"
+#include "../modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h"
+#include "../modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h"
+#include "../modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h"
+#include "../modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h"
+#include "../modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h"
+#include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
 
 
@@ -15,28 +29,28 @@
 class Numberedcostfunction: public Object, public Definition{
 
-public: 
+	public: 
 
-int   definitionenum;
-char* name;
-int   number_cost_functions;
-int*  cost_functions_list;
+		int   definitionenum;
+		char* name;
+		int   number_cost_functions;
+		int*  cost_functions_list;
+		
+		/*Numberedcostfunction constructors, destructors :*/
+		Numberedcostfunction();
+		Numberedcostfunction(char* in_name, int in_definitionenum,int number_cost_functions_in,int* cost_functions_list_in);
+		~Numberedcostfunction();
 
-/*Numberedcostfunction constructors, destructors :*/
-Numberedcostfunction();
-Numberedcostfunction(char* in_name, int in_definitionenum,int number_cost_functions_in,int* cost_functions_list_in);
-~Numberedcostfunction();
+		/*Object virtual function resolutoin: */
+		Object*	copy();
+		void		DeepEcho(void);
+		void		Echo(void);
+		int		Id(void);
+		void		Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
+		int		ObjectEnum(void);
 
-/*Object virtual function resolutoin: */
-Object*	copy();
-void		DeepEcho(void);
-void		Echo(void);
-int		Id(void);
-void		Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
-int		ObjectEnum(void);
-
-/*Definition virtual function resolutoin: */
-int		DefinitionEnum();
-char*		Name();
-IssmDouble Response(FemModel* femmodel);
+		/*Definition virtual function resolutoin: */
+		int		DefinitionEnum();
+		char*		Name();
+		IssmDouble Response(FemModel* femmodel);
 };
 
Index: issm/trunk-jpl/src/c/classes/Regionaloutput.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Regionaloutput.cpp	(revision 22437)
+++ issm/trunk-jpl/src/c/classes/Regionaloutput.cpp	(revision 22438)
@@ -15,4 +15,5 @@
 #include "./Elements/Element.h"
 #include "./Elements/Elements.h"
+#include "./FemModel.h"
 #include "../classes/Params/Parameters.h"
 
Index: issm/trunk-jpl/src/c/classes/Regionaloutput.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Regionaloutput.h	(revision 22437)
+++ issm/trunk-jpl/src/c/classes/Regionaloutput.h	(revision 22438)
@@ -9,5 +9,9 @@
 /*{{{*/
 #include "./Definition.h"
+#include "../datastructures/datastructures.h"
+#include "./Elements/Element.h"
+#include "./Elements/Elements.h"
 #include "./FemModel.h"
+#include "../classes/Params/Parameters.h"
 
 /*}}}*/
Index: issm/trunk-jpl/src/c/cores/ad_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/ad_core.cpp	(revision 22437)
+++ issm/trunk-jpl/src/c/cores/ad_core.cpp	(revision 22438)
@@ -26,5 +26,5 @@
 	int     num_dependents=0;
 	int     num_independents=0;
-	bool    isautodiff,iscontrol;
+	bool    isautodiff       = false;
 	char   *driver           = NULL;
 	size_t  tape_stats[15];
@@ -37,7 +37,6 @@
 	/*AD mode on?: */
 	femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
-	femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum);
-
-	if(isautodiff && !iscontrol){
+
+	if(isautodiff){
 
 		#ifdef _HAVE_ADOLC_
Index: issm/trunk-jpl/src/c/cores/controlADm1qn3_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/controlADm1qn3_core.cpp	(revision 22437)
+++ 	(revision )
@@ -1,470 +1,0 @@
-/*!\file: controlm1qn3_core.cpp
- * \brief: core of the control solution 
- */ 
-
-#include <config.h>
-#include "./cores.h"
-#include "../toolkits/toolkits.h"
-#include "../classes/classes.h"
-#include "../shared/shared.h"
-#include "../modules/modules.h"
-#include "../solutionsequences/solutionsequences.h"
-
-#if defined (_HAVE_M1QN3_) 
-//& !defined(_HAVE_ADOLC_)
-/*m1qn3 prototypes*/
-extern "C" void *ctonbe_; // DIS mode : Conversion
-extern "C" void *ctcabe_; // DIS mode : Conversion
-extern "C" void *euclid_; // Scalar product
-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 [],void* );
-
-/*Use struct to provide arguments*/
-typedef struct{
-	FemModel   * femmodel;
-	IssmPDouble* Jlist;
-	int          M;
-	int          N;
-	int*         i;
-} m1qn3_struct;
-
-/*m1qm3 functions*/
-void simul_starttrace(FemModel* femmodel){/*{{{*/
-
-	/*Retrive ADOLC parameters*/
-	IssmDouble gcTriggerRatio;
-	IssmDouble gcTriggerMaxSize;
-	IssmDouble obufsize;
-	IssmDouble lbufsize;
-	IssmDouble cbufsize;
-	IssmDouble tbufsize;
-	femmodel->parameters->FindParam(&gcTriggerRatio,AutodiffGcTriggerRatioEnum);
-	femmodel->parameters->FindParam(&gcTriggerMaxSize,AutodiffGcTriggerMaxSizeEnum);
-	femmodel->parameters->FindParam(&obufsize,AutodiffObufsizeEnum);
-	femmodel->parameters->FindParam(&lbufsize,AutodiffLbufsizeEnum);
-	femmodel->parameters->FindParam(&cbufsize,AutodiffCbufsizeEnum);
-	femmodel->parameters->FindParam(&tbufsize,AutodiffTbufsizeEnum);
-
-	/*Set garbage collection parameters: */
-	setStoreManagerControl(reCast<IssmPDouble>(gcTriggerRatio),reCast<size_t>(gcTriggerMaxSize));
-
-	/*Start trace: */
-	int skipFileDeletion=1;
-	int keepTaylors=1;
-	int my_rank=IssmComm::GetRank();
-	trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
-}/*}}}*/
-void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
-
-	/*Get rank*/
-	int my_rank=IssmComm::GetRank();
-
-	/*Recover Arguments*/
-	m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
-	FemModel     *femmodel     = input_struct->femmodel;
-	IssmPDouble  *Jlist        = input_struct->Jlist;
-	int           JlistM       = input_struct->M;
-	int           JlistN       = input_struct->N;
-	int          *Jlisti       = input_struct->i;
-	int           intn         = (int)*n;
-
-	/*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 and update controls*/
-	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 i=0;i<numberofvertices;i++){
-		for(int c=0;c<num_controls;c++){
-			int index = num_controls*i+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];
-		}
-	}
-
-	/*Start Tracing*/
-	simul_starttrace(femmodel);
-
-	/*Set X as our new control input abd as INDEPENDENT!!*/
-#ifdef _HAVE_AD_
-	IssmDouble* aX=xNew<IssmDouble>(num_controls*numberofvertices,"t");
-#else
-	IssmDouble* aX=xNew<IssmDouble>(num_controls*numberofvertices);
-#endif
-	if(my_rank==0){
-		for(int i=0;i<intn;i++){
-			aX[i]<<=X[i];
-		}
-	}
-	ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-	SetControlInputsFromVectorx(femmodel,aX);
-	xDelete<IssmDouble>(aX);
-
-	/*Compute solution (forward)*/
-	void (*solutioncore)(FemModel*)=NULL;
-	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
-	solutioncore(femmodel);
-
-	/*Get Dependents*/
-	IssmDouble  output_value;
-	int         num_dependents;
-	IssmPDouble *dependents;
-	DataSet*    dependent_objects=NULL;
-
-	femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
-	femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
-
-	/*Go through our dependent variables, and compute the response:*/
-	dependents=xNew<IssmPDouble>(num_dependents);
-	for(int i=0;i<dependent_objects->Size();i++){
-		DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
-		dep->Responsex(&output_value,femmodel);
-		if (my_rank==0) {
-			output_value>>=dependents[i];
-		}
-	}
-
-	/*Turning off trace tape*/
-	trace_off();
-
-	/*Print tape statistics so that user can kill this run if something is off already:*/
-	if(VerboseAutodiff()){ /*{{{*/
-		size_t  tape_stats[15];
-		tapestats(my_rank,tape_stats); //reading of tape statistics
-		int commSize=IssmComm::GetSize();
-		int *sstats=new int[7];
-		sstats[0]=tape_stats[NUM_OPERATIONS];
-		sstats[1]=tape_stats[OP_FILE_ACCESS];
-		sstats[2]=tape_stats[NUM_LOCATIONS];
-		sstats[3]=tape_stats[LOC_FILE_ACCESS];
-		sstats[4]=tape_stats[NUM_VALUES];
-		sstats[5]=tape_stats[VAL_FILE_ACCESS];
-		sstats[6]=tape_stats[TAY_STACK_SIZE];
-		int *rstats=NULL;
-		if (my_rank==0) rstats=new int[commSize*7];
-		ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
-		if (my_rank==0) {
-			int offset=50;
-			int rOffset=(commSize/10)+1;
-			_printf_("   ADOLC statistics: \n");
-			_printf_("     "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
-			_printf_("     "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
-			_printf_("     "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
-			_printf_("     operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
-			_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
-			for (int r=0;r<commSize;++r)
-			 _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+0] << (rstats[r*7+1]?" ->file":"") << "\n");
-			_printf_("     locations: entry size " << sizeof(locint) << " Bytes\n");
-			_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
-			for (int r=0;r<commSize;++r)
-			 _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+2] << (rstats[r*7+3]?" ->file":"") << "\n");
-			_printf_("     constant values: entry size " << sizeof(double) << " Bytes\n");
-			_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
-			for (int r=0;r<commSize;++r)
-			 _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+4] << (rstats[r*7+5]?" ->file":"") << "\n");
-			_printf_("     Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
-			_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
-			for (int r=0;r<commSize;++r)
-			 _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+6] << (rstats[r*7+6]>tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n");
-			delete []rstats;
-		}
-		delete [] sstats;
-	} /*}}}*/
-
-	/*diverse: */
-	int  dummy;
-	int  num_independents=0;
-
-	/*intermediary: */
-	IssmPDouble *aWeightVector=NULL;
-	IssmPDouble *weightVectorTimesJac=NULL;
-
-	/*output: */
-	IssmPDouble *totalgradient=NULL;
-
-	/*retrieve parameters: */
-	num_independents = intn;
-
-	/*if no dependents, no point in running a driver: */
-	if(!(num_dependents*num_independents)) _error_("this is not allowed");
-
-	/*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/
-	int num_dependents_old   = num_dependents;
-	int num_independents_old = num_independents;
-	if(my_rank!=0){
-		num_dependents   = 0;
-		num_independents = 0;
-	}
-
-	/*get the EDF pointer:*/
-	ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
-
-	/* these are always needed regardless of the interpreter */
-	anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
-	anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
-
-	/* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so
-	 * as to generate num_dependents gradients: */
-	totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
-
-	for(int aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
-
-		/*initialize direction index in the weights vector: */
-		aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
-		if (my_rank==0) aWeightVector[aDepIndex]=1.;
-
-		/*initialize output gradient: */
-		weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
-
-		/*set the forward method function pointer: */
-		#ifdef _HAVE_GSL_
-		anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
-		#endif
-		#ifdef _HAVE_MUMPS_
-		anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
-		#endif
-
-		anEDF_for_solverx_p->dp_U=xNew<IssmPDouble>(anEDF_for_solverx_p->max_m);
-		anEDF_for_solverx_p->dp_Z=xNew<IssmPDouble>(anEDF_for_solverx_p->max_n);
-
-		/*call driver: */
-		fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
-
-		/*Add to totalgradient: */
-		if(my_rank==0) for(int i=0;i<num_independents;i++) totalgradient[i]+=weightVectorTimesJac[i];
-
-		/*free resources :*/
-		xDelete(weightVectorTimesJac);
-		xDelete(aWeightVector);
-	}
-
-	/*Broadcast gradient to other ranks*/
-	ISSM_MPI_Bcast(totalgradient,num_independents_old,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
-
-	/*Check size of Jlist to avoid crashes*/
-	_assert_((*Jlisti)<JlistM);
-	_assert_(JlistN==num_responses+1);
-
-	/*Compute objective function*/
-	IssmDouble* 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] = reCast<IssmPDouble>(Jtemp[i]);
-	Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J);
-	xDelete<IssmDouble>(Jtemp);
-
-	if(*indic==0){
-		/*dry run, no gradient required*/
-
-		/*Retrieve objective functions independently*/
-		_printf0_("            N/A |\n");
-		for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
-		_printf0_("\n");
-
-		*Jlisti = (*Jlisti) +1;
-		xDelete<double>(XU);
-		xDelete<double>(XL);
-		return;
-	}
-
-	/*Compute gradient*/
-	for(long i=0;i<num_independents_old;i++) G[i] = totalgradient[i];
-
-	/*Constrain Gradient*/
-	IssmDouble  Gnorm = 0.;
-	for(int i=0;i<numberofvertices;i++){
-		for(int c=0;c<num_controls;c++){
-			int index = num_controls*i+c;
-			if(X[index]>=XU[index]) G[index]=0.;
-			if(X[index]<=XL[index]) G[index]=0.;
-			G[index] = G[index]*reCast<double>(scaling_factors[c]);
-			X[index] = X[index]/reCast<double>(scaling_factors[c]);
-			Gnorm += G[index]*G[index];
-		}
-	}
-	Gnorm = sqrt(Gnorm);
-
-	/*Print info*/
-	_printf0_("       "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
-	for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
-	_printf0_("\n");
-
-	/*Clean-up and return*/
-	*Jlisti = (*Jlisti) +1;
-	xDelete<double>(XU);
-	xDelete<double>(XL);
-	xDelete<IssmDouble>(scaling_factors);
-	xDelete<IssmPDouble>(totalgradient);
-
-}/*}}}*/
-void controlm1qn3_core(FemModel* femmodel){/*{{{*/
-
-	/*Intermediaries*/
-	long         omode;
-	double		 f;
-	double		 dxmin,gttol; 
-	int          maxsteps,maxiter;
-	int          intn,numberofvertices,num_controls,num_cost_functions,solution_type;
-	IssmDouble	*scaling_factors = NULL;
-	double      *X  = NULL;
-	double      *G  = NULL;
-
-	/*Recover some parameters*/
-	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);
-	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();
-
-	/*Initialize M1QN3 parameters*/
-	if(VerboseControl())_printf0_("   Initialize M1QN3 parameters\n");
-	SimulFunc simul_ptr    = &simul_ad; /*Cost function address*/
-	void**    prosca       = &euclid_;  /*Dot product function (euclid is the default)*/
-	char      normtype[]   = "dfn";     /*Norm type: dfn = scalar product defined by prosca*/
-	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*/
-	long      imode[3]     = {0};       /*scaling and starting mode, 0 by default*/
-	long      indic        = 4;         /*compute f and g*/
-	long      reverse      = 0;         /*reverse or direct mode*/
-	long      io           = 6;         /*Channel number for the output*/
-
-	/*Optimization criterions*/
-	long niter = long(maxsteps); /*Maximum number of iterations*/
-	long nsim  = long(maxiter);/*Maximum number of function calls*/
-
-	/*Get initial guess*/
-	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;
-	_assert_(intn==numberofvertices*num_controls);
-
-	/*Get problem dimension and initialize gradient and initial guess*/
-	long n = long(intn);
-	G = xNew<double>(n);
-
-	/*Scale control for M1QN3*/
-	for(int i=0;i<numberofvertices;i++){
-		for(int c=0;c<num_controls;c++){
-			int index = num_controls*i+c;
-			X[index] = X[index]/reCast<IssmPDouble>(scaling_factors[c]);
-		}
-	}
-
-	/*Allocate m1qn3 working arrays (see documentation)*/
-	long      m   = 100;
-	long      ndz = 4*n+m*(2*n+1);
-	double*   dz  = xNew<double>(ndz);
-
-	if(VerboseControl())_printf0_("   Computing initial solution\n");
-	_printf0_("\n");
-	_printf0_("Cost function f(x)   | Gradient norm |g(x)| |  List of contributions\n");
-	_printf0_("____________________________________________________________________\n");
-
-	/*Prepare structure for m1qn3*/
-	m1qn3_struct mystruct;
-	mystruct.femmodel = femmodel;
-	mystruct.M        = maxiter;
-	mystruct.N        = num_cost_functions+1;
-	mystruct.Jlist    = xNewZeroInit<IssmPDouble>(mystruct.M*mystruct.N);
-	mystruct.i        = xNewZeroInit<int>(1);
-	/*Initialize Gradient and cost function of M1QN3*/
-	indic = 4; /*gradient required*/
-	simul_ad(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
-	/*Estimation of the expected decrease in f during the first iteration*/
-	double df1=f;
-
-	/*Call M1QN3 solver*/
-	m1qn3_(simul_ptr,prosca,&ctonbe_,&ctcabe_,
-				&n,X,&f,G,&dxmin,&df1,
-				&gttol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz,
-				&reverse,&indic,izs,rzs,(void*)&mystruct);
-	switch(int(omode)){
-		case 0:  _printf0_("   Stop requested (indic = 0)\n"); break;
-		case 1:  _printf0_("   Convergence reached (gradient satisfies stopping criterion)\n"); break;
-		case 2:  _printf0_("   Bad initialization\n"); break;
-		case 3:  _printf0_("   Line search failure\n"); break;
-		case 4:  _printf0_("   Maximum number of iterations exceeded\n");break;
-		case 5:  _printf0_("   Maximum number of function calls exceeded\n"); break;
-		case 6:  _printf0_("   stopped on dxmin during line search\n"); break;
-		case 7:  _printf0_("   <g,d> > 0  or  <y,s> <0\n"); break;
-		default: _printf0_("   Unknown end condition\n");
-	}
-	/*Constrain solution vector*/
-	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 i=0;i<numberofvertices;i++){
-		for(int c=0;c<num_controls;c++){
-			int index = num_controls*i+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*/
-	IssmDouble* aX=xNew<IssmDouble>(intn);
-	for(int i=0;i<intn;i++) aX[i] = reCast<IssmDouble>(X[i]);
-	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);
-	xDelete(aG);
-
-	/*Add last cost function to results*/
-	femmodel->OutputControlsx(&femmodel->results);
-	femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,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>(dz);
-	xDelete<double>(XU);
-	xDelete<double>(XL);
-	xDelete<IssmDouble>(scaling_factors);
-	xDelete<IssmPDouble>(mystruct.Jlist);
-	xDelete<int>(mystruct.i);
-}/*}}}*/
-
-#else
-void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");}
-#endif //_HAVE_M1QN3_
Index: issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 22437)
+++ issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 22438)
@@ -11,319 +11,38 @@
 #include "../solutionsequences/solutionsequences.h"
 
-#if defined (_HAVE_M1QN3_) 
-//& !defined(_HAVE_ADOLC_)
+#if defined (_HAVE_M1QN3_) & !defined(_HAVE_ADOLC_)
 /*m1qn3 prototypes*/
 extern "C" void *ctonbe_; // DIS mode : Conversion
 extern "C" void *ctcabe_; // DIS mode : Conversion
 extern "C" void *euclid_; // Scalar product
-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),
+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 *, double [], double *, double [], double*, double *,
+			double *, char [], long *, long *, long *, long *, long *, long *, long [], double [], long *,
 			long *, long *, long [], float [],void* );
+
+/*Cost function prototype*/
+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;
-	IssmPDouble* Jlist;
-	int          M;
-	int          N;
-	int*         i;
+	FemModel  * femmodel;
+	IssmDouble* Jlist;
+	int         M;
+	int         N;
+	int*        i;
 } m1qn3_struct;
 
-/*m1qm3 functions*/
-void simul_starttrace(FemModel* femmodel){/*{{{*/
-
-	/*Retrive ADOLC parameters*/
-	IssmDouble gcTriggerRatio;
-	IssmDouble gcTriggerMaxSize;
-	IssmDouble obufsize;
-	IssmDouble lbufsize;
-	IssmDouble cbufsize;
-	IssmDouble tbufsize;
-	femmodel->parameters->FindParam(&gcTriggerRatio,AutodiffGcTriggerRatioEnum);
-	femmodel->parameters->FindParam(&gcTriggerMaxSize,AutodiffGcTriggerMaxSizeEnum);
-	femmodel->parameters->FindParam(&obufsize,AutodiffObufsizeEnum);
-	femmodel->parameters->FindParam(&lbufsize,AutodiffLbufsizeEnum);
-	femmodel->parameters->FindParam(&cbufsize,AutodiffCbufsizeEnum);
-	femmodel->parameters->FindParam(&tbufsize,AutodiffTbufsizeEnum);
-
-	/*Set garbage collection parameters: */
-	setStoreManagerControl(reCast<IssmPDouble>(gcTriggerRatio),reCast<size_t>(gcTriggerMaxSize));
-
-	/*Start trace: */
-	int skipFileDeletion=1;
-	int keepTaylors=1;
-	int my_rank=IssmComm::GetRank();
-	trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
-}/*}}}*/
-void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
-
-	/*Get rank*/
-	int my_rank=IssmComm::GetRank();
-
-	/*Recover Arguments*/
-	m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
-	FemModel     *femmodel     = input_struct->femmodel;
-	IssmPDouble  *Jlist        = input_struct->Jlist;
-	int           JlistM       = input_struct->M;
-	int           JlistN       = input_struct->N;
-	int          *Jlisti       = input_struct->i;
-	int           intn         = (int)*n;
-
-	/*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 and update controls*/
-	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 i=0;i<numberofvertices;i++){
-		for(int c=0;c<num_controls;c++){
-			int index = num_controls*i+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];
-		}
-	}
-
-	/*Start Tracing*/
-	simul_starttrace(femmodel);
-
-	/*Set X as our new control input abd as INDEPENDENT!!*/
-#ifdef _HAVE_AD_
-	IssmDouble* aX=xNew<IssmDouble>(num_controls*numberofvertices,"t");
-#else
-	IssmDouble* aX=xNew<IssmDouble>(num_controls*numberofvertices);
-#endif
-	if(my_rank==0){
-		for(int i=0;i<intn;i++){
-			aX[i]<<=X[i];
-		}
-	}
-	ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-	SetControlInputsFromVectorx(femmodel,aX);
-	xDelete<IssmDouble>(aX);
-
-	/*Compute solution (forward)*/
-	void (*solutioncore)(FemModel*)=NULL;
-	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
-	solutioncore(femmodel);
-
-	/*Get Dependents*/
-	IssmDouble  output_value;
-	int         num_dependents;
-	IssmPDouble *dependents;
-	DataSet*    dependent_objects=NULL;
-
-	femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
-	femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
-
-	/*Go through our dependent variables, and compute the response:*/
-	dependents=xNew<IssmPDouble>(num_dependents);
-	for(int i=0;i<dependent_objects->Size();i++){
-		DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
-		dep->Responsex(&output_value,femmodel);
-		if (my_rank==0) {
-			output_value>>=dependents[i];
-		}
-	}
-
-	/*Turning off trace tape*/
-	trace_off();
-
-	/*Print tape statistics so that user can kill this run if something is off already:*/
-	if(VerboseAutodiff()){ /*{{{*/
-		size_t  tape_stats[15];
-		tapestats(my_rank,tape_stats); //reading of tape statistics
-		int commSize=IssmComm::GetSize();
-		int *sstats=new int[7];
-		sstats[0]=tape_stats[NUM_OPERATIONS];
-		sstats[1]=tape_stats[OP_FILE_ACCESS];
-		sstats[2]=tape_stats[NUM_LOCATIONS];
-		sstats[3]=tape_stats[LOC_FILE_ACCESS];
-		sstats[4]=tape_stats[NUM_VALUES];
-		sstats[5]=tape_stats[VAL_FILE_ACCESS];
-		sstats[6]=tape_stats[TAY_STACK_SIZE];
-		int *rstats=NULL;
-		if (my_rank==0) rstats=new int[commSize*7];
-		ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
-		if (my_rank==0) {
-			int offset=50;
-			int rOffset=(commSize/10)+1;
-			_printf_("   ADOLC statistics: \n");
-			_printf_("     "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
-			_printf_("     "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
-			_printf_("     "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
-			_printf_("     operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
-			_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
-			for (int r=0;r<commSize;++r)
-			 _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+0] << (rstats[r*7+1]?" ->file":"") << "\n");
-			_printf_("     locations: entry size " << sizeof(locint) << " Bytes\n");
-			_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
-			for (int r=0;r<commSize;++r)
-			 _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+2] << (rstats[r*7+3]?" ->file":"") << "\n");
-			_printf_("     constant values: entry size " << sizeof(double) << " Bytes\n");
-			_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
-			for (int r=0;r<commSize;++r)
-			 _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+4] << (rstats[r*7+5]?" ->file":"") << "\n");
-			_printf_("     Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
-			_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
-			for (int r=0;r<commSize;++r)
-			 _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+6] << (rstats[r*7+6]>tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n");
-			delete []rstats;
-		}
-		delete [] sstats;
-	} /*}}}*/
-
-	/*diverse: */
-	int  dummy;
-	int  num_independents=0;
-
-	/*intermediary: */
-	IssmPDouble *aWeightVector=NULL;
-	IssmPDouble *weightVectorTimesJac=NULL;
-
-	/*output: */
-	IssmPDouble *totalgradient=NULL;
-
-	/*retrieve parameters: */
-	num_independents = intn;
-
-	/*if no dependents, no point in running a driver: */
-	if(!(num_dependents*num_independents)) _error_("this is not allowed");
-
-	/*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/
-	int num_dependents_old   = num_dependents;
-	int num_independents_old = num_independents;
-	if(my_rank!=0){
-		num_dependents   = 0;
-		num_independents = 0;
-	}
-
-	/*get the EDF pointer:*/
-	ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
-
-	/* these are always needed regardless of the interpreter */
-	anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
-	anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
-
-	/* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so
-	 * as to generate num_dependents gradients: */
-	totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
-
-	for(int aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
-
-		/*initialize direction index in the weights vector: */
-		aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
-		if (my_rank==0) aWeightVector[aDepIndex]=1.;
-
-		/*initialize output gradient: */
-		weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
-
-		/*set the forward method function pointer: */
-		#ifdef _HAVE_GSL_
-		anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
-		#endif
-		#ifdef _HAVE_MUMPS_
-		anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
-		#endif
-
-		anEDF_for_solverx_p->dp_U=xNew<IssmPDouble>(anEDF_for_solverx_p->max_m);
-		anEDF_for_solverx_p->dp_Z=xNew<IssmPDouble>(anEDF_for_solverx_p->max_n);
-
-		/*call driver: */
-		fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
-
-		/*Add to totalgradient: */
-		if(my_rank==0) for(int i=0;i<num_independents;i++) totalgradient[i]+=weightVectorTimesJac[i];
-
-		/*free resources :*/
-		xDelete(weightVectorTimesJac);
-		xDelete(aWeightVector);
-	}
-
-	/*Broadcast gradient to other ranks*/
-	ISSM_MPI_Bcast(totalgradient,num_independents_old,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
-
-	/*Check size of Jlist to avoid crashes*/
-	_assert_((*Jlisti)<JlistM);
-	_assert_(JlistN==num_responses+1);
-
-	/*Compute objective function*/
-	IssmDouble* 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] = reCast<IssmPDouble>(Jtemp[i]);
-	Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J);
-	xDelete<IssmDouble>(Jtemp);
-
-	if(*indic==0){
-		/*dry run, no gradient required*/
-
-		/*Retrieve objective functions independently*/
-		_printf0_("            N/A |\n");
-		for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
-		_printf0_("\n");
-
-		*Jlisti = (*Jlisti) +1;
-		xDelete<double>(XU);
-		xDelete<double>(XL);
-		return;
-	}
-
-	/*Compute gradient*/
-	for(long i=0;i<num_independents_old;i++) G[i] = totalgradient[i];
-
-	/*Constrain Gradient*/
-	IssmDouble  Gnorm = 0.;
-	for(int i=0;i<numberofvertices;i++){
-		for(int c=0;c<num_controls;c++){
-			int index = num_controls*i+c;
-			if(X[index]>=XU[index]) G[index]=0.;
-			if(X[index]<=XL[index]) G[index]=0.;
-			G[index] = G[index]*reCast<double>(scaling_factors[c]);
-			X[index] = X[index]/reCast<double>(scaling_factors[c]);
-			Gnorm += G[index]*G[index];
-		}
-	}
-	Gnorm = sqrt(Gnorm);
-
-	/*Print info*/
-	_printf0_("       "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
-	for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
-	_printf0_("\n");
-
-	/*Clean-up and return*/
-	*Jlisti = (*Jlisti) +1;
-	xDelete<double>(XU);
-	xDelete<double>(XL);
-	xDelete<IssmDouble>(scaling_factors);
-	xDelete<IssmPDouble>(totalgradient);
-
-}/*}}}*/
-void controlm1qn3_core(FemModel* femmodel){/*{{{*/
+void controlm1qn3_core(FemModel* femmodel){
 
 	/*Intermediaries*/
 	long         omode;
-	double		 f;
-	double		 dxmin,gttol; 
+	double       f,dxmin,gttol; 
 	int          maxsteps,maxiter;
 	int          intn,numberofvertices,num_controls,num_cost_functions,solution_type;
-	IssmDouble	*scaling_factors = NULL;
-	double      *X  = NULL;
-	double      *G  = NULL;
+	IssmDouble  *scaling_factors = NULL;
+	IssmDouble  *X  = NULL;
+	IssmDouble  *G  = NULL;
 
 	/*Recover some parameters*/
@@ -333,6 +52,6 @@
 	femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum);
 	femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
-	femmodel->parameters->FindParamAndMakePassive(&dxmin,InversionDxminEnum);
-	femmodel->parameters->FindParamAndMakePassive(&gttol,InversionGttolEnum);
+	femmodel->parameters->FindParam(&dxmin,InversionDxminEnum);
+	femmodel->parameters->FindParam(&gttol,InversionGttolEnum);
 	femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
 	femmodel->parameters->SetParam(false,SaveResultsEnum);
@@ -341,5 +60,5 @@
 	/*Initialize M1QN3 parameters*/
 	if(VerboseControl())_printf0_("   Initialize M1QN3 parameters\n");
-	SimulFunc simul_ptr    = &simul_ad; /*Cost function address*/
+	SimulFunc costfuncion  = &simul;    /*Cost function address*/
 	void**    prosca       = &euclid_;  /*Dot product function (euclid is the default)*/
 	char      normtype[]   = "dfn";     /*Norm type: dfn = scalar product defined by prosca*/
@@ -358,6 +77,6 @@
 
 	/*Get initial guess*/
-	Vector<double> *Xpetsc = NULL;
-	GetPassiveVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
+	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);
@@ -373,5 +92,5 @@
 		for(int c=0;c<num_controls;c++){
 			int index = num_controls*i+c;
-			X[index] = X[index]/reCast<IssmPDouble>(scaling_factors[c]);
+			X[index] = X[index]/scaling_factors[c];
 		}
 	}
@@ -392,17 +111,20 @@
 	mystruct.M        = maxiter;
 	mystruct.N        = num_cost_functions+1;
-	mystruct.Jlist    = xNewZeroInit<IssmPDouble>(mystruct.M*mystruct.N);
+	mystruct.Jlist    = xNewZeroInit<IssmDouble>(mystruct.M*mystruct.N);
 	mystruct.i        = xNewZeroInit<int>(1);
+
 	/*Initialize Gradient and cost function of M1QN3*/
 	indic = 4; /*gradient required*/
-	simul_ad(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
+	simul(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
+
 	/*Estimation of the expected decrease in f during the first iteration*/
 	double df1=f;
 
 	/*Call M1QN3 solver*/
-	m1qn3_(simul_ptr,prosca,&ctonbe_,&ctcabe_,
+	m1qn3_(costfuncion,prosca,&ctonbe_,&ctcabe_,
 				&n,X,&f,G,&dxmin,&df1,
 				&gttol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz,
 				&reverse,&indic,izs,rzs,(void*)&mystruct);
+
 	switch(int(omode)){
 		case 0:  _printf0_("   Stop requested (indic = 0)\n"); break;
@@ -416,34 +138,22 @@
 		default: _printf0_("   Unknown end condition\n");
 	}
+
 	/*Constrain solution vector*/
-	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 i=0;i<numberofvertices;i++){
-		for(int c=0;c<num_controls;c++){
-			int index = num_controls*i+c;
-			X[index] = X[index]*reCast<double>(scaling_factors[c]);
+	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");
+	for(int i=0;i<numberofvertices;i++){
+		for(int c=0;c<num_controls;c++){
+			int index = num_controls*i+c;
+			X[index] = X[index]*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*/
-	IssmDouble* aX=xNew<IssmDouble>(intn);
-	for(int i=0;i<intn;i++) aX[i] = reCast<IssmDouble>(X[i]);
-	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);
-	xDelete(aG);
-
-	/*Add last cost function to results*/
+	SetControlInputsFromVectorx(femmodel,X);
+	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
 	femmodel->OutputControlsx(&femmodel->results);
-	femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0));
+	femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0));
 
 	/*Finalize*/
@@ -454,5 +164,4 @@
 	solutioncore(femmodel);
 
-
 	/*Clean-up and return*/
 	xDelete<double>(G);
@@ -461,10 +170,117 @@
 	xDelete<double>(XU);
 	xDelete<double>(XL);
+	xDelete<double>(scaling_factors);
+	xDelete<double>(mystruct.Jlist);
+	xDelete<int>(mystruct.i);
+}
+
+/*Cost function definition*/
+void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){
+
+	/*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 and update controls*/
+	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");
+	for(int i=0;i<numberofvertices;i++){
+		for(int c=0;c<num_controls;c++){
+			int index = num_controls*i+c;
+			X[index] = X[index]*scaling_factors[c];
+			if(X[index]>XU[index]) X[index]=XU[index];
+			if(X[index]<XL[index]) X[index]=XL[index];
+		}
+	}
+	SetControlInputsFromVectorx(femmodel,X);
+
+	/*Compute solution and adjoint*/
+	void (*solutioncore)(FemModel*)=NULL;
+	void (*adjointcore)(FemModel*)=NULL;
+	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
+	solutioncore(femmodel);
+
+	/*Check size of Jlist to avoid crashes*/
+	_assert_((*Jlisti)<JlistM);
+	_assert_(JlistN==num_responses+1);
+
+	/*Compute objective function*/
+	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){
+		/*dry run, no gradient required*/
+
+		/*Retrieve objective functions independently*/
+		_printf0_("            N/A |\n");
+		for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
+		_printf0_("\n");
+
+		*Jlisti = (*Jlisti) +1;
+		xDelete<IssmDouble>(XU);
+		xDelete<IssmDouble>(XL);
+		return;
+	}
+
+	/*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];
+	xDelete<IssmDouble>(G2);
+
+	/*Constrain Gradient*/
+	IssmDouble  Gnorm = 0.;
+	for(int i=0;i<numberofvertices;i++){
+		for(int c=0;c<num_controls;c++){
+			int index = num_controls*i+c;
+			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];
+			Gnorm += G[index]*G[index];
+		}
+	}
+	Gnorm = sqrt(Gnorm);
+
+	/*Print info*/
+	_printf0_("       "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
+	for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
+	_printf0_("\n");
+
+	/*Clean-up and return*/
+	*Jlisti = (*Jlisti) +1;
+	xDelete<IssmDouble>(XU);
+	xDelete<IssmDouble>(XL);
 	xDelete<IssmDouble>(scaling_factors);
-	xDelete<IssmPDouble>(mystruct.Jlist);
-	xDelete<int>(mystruct.i);
-}/*}}}*/
+}
 
 #else
-void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");}
+void controlm1qn3_core(FemModel* femmodel){
+	_error_("M1QN3 not installed");
+}
 #endif //_HAVE_M1QN3_
Index: issm/trunk-jpl/src/c/cores/stressbalance_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/stressbalance_core.cpp	(revision 22437)
+++ issm/trunk-jpl/src/c/cores/stressbalance_core.cpp	(revision 22438)
@@ -79,6 +79,6 @@
 		femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
 	}
-	
-//	if(solution_type==StressbalanceSolutionEnum)femmodel->RequestedDependentsx();
+
+	if(solution_type==StressbalanceSolutionEnum)femmodel->RequestedDependentsx();
 
 	/*Free ressources:*/	
Index: issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp	(revision 22437)
+++ issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp	(revision 22438)
@@ -59,5 +59,5 @@
 
 	/*Retrieve all inputs we will be needing: */
-//	Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+	Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
 	Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
 	Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
@@ -79,5 +79,5 @@
 
 		/*Get all parameters at gaussian point*/
-		//weights_input->GetInputValue(&weight,gauss,SurfaceAbsVelMisfitEnum);
+		weights_input->GetInputValue(&weight,gauss,SurfaceAbsVelMisfitEnum);
 		vx_input->GetInputValue(&vx,gauss);
 		vxobs_input->GetInputValue(&vxobs,gauss);
@@ -98,7 +98,5 @@
 
 		/*Add to cost function*/
-		//Jelem+=misfit*weight*Jdet*gauss->weight;
-		Jelem+=misfit*Jdet*gauss->weight;
-
+		Jelem+=misfit*weight*Jdet*gauss->weight;
 	}
 
Index: issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp	(revision 22437)
+++ issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp	(revision 22438)
@@ -62,5 +62,5 @@
 
 	/*Retrieve all inputs we will be needing: */
-//	Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+	Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
 	Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
 	Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
@@ -82,5 +82,5 @@
 
 		/*Get all parameters at gaussian point*/
-		//weights_input->GetInputValue(&weight,gauss,SurfaceLogVelMisfitEnum);
+		weights_input->GetInputValue(&weight,gauss,SurfaceLogVelMisfitEnum);
 		vx_input->GetInputValue(&vx,gauss);
 		vxobs_input->GetInputValue(&vxobs,gauss);
@@ -106,5 +106,5 @@
 
 		misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
-		weight = 1.;
+
 		/*Add to cost function*/
 		Jelem+=misfit*weight*Jdet*gauss->weight;
Index: issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp	(revision 22437)
+++ issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp	(revision 22438)
@@ -61,5 +61,5 @@
 
 	/*Retrieve all inputs we will be needing: */
-//	Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+	Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
 	Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
 	Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
@@ -81,5 +81,5 @@
 
 		/*Get all parameters at gaussian point*/
-		//weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum);
+		weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum);
 		vx_input->GetInputValue(&vx,gauss);
 		vxobs_input->GetInputValue(&vxobs,gauss);
@@ -106,5 +106,5 @@
 			misfit=0.5*(scalex*pow((vx-vxobs),2));
 		}
-		weight = 1.;
+
 		/*Add to cost function*/
 		Jelem+=misfit*weight*Jdet*gauss->weight;
