Index: ../trunk-jpl/src/c/cores/controlad_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/controlad_core.cpp	(revision 18941)
+++ ../trunk-jpl/src/c/cores/controlad_core.cpp	(revision 18942)
@@ -24,33 +24,30 @@
 
 /*Cost function prototype*/
 void simulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs);
+FemModel* presimulad(int* pintn, double** pX, FemModel* femmodel);
+void postsimulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs);
 /*}}}*/
 void controlad_core(FemModel* femmodel){ /*{{{*/
 
 	/*Intermediaries*/
+	FemModel*    femmodelad=NULL;
 	int          i;
 	long         omode;
 	IssmPDouble  f,dxmin,gttol;
 	IssmDouble   dxmind,gttold; 
 	int          maxsteps,maxiter;
-	int          intn,numberofvertices,num_controls,solution_type;
-	IssmDouble  *scaling_factors = NULL;
+	int          intn,solution_type;
 	IssmPDouble  *X  = NULL;
 	IssmDouble   *Xd  = NULL;
-	IssmDouble   *Gd  = NULL;
 	IssmPDouble  *G  = NULL;
-	bool onsid=true;
 
 	/*Recover some parameters*/
 	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
 	femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum);
 	femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
 	femmodel->parameters->FindParam(&dxmind,InversionDxminEnum); dxmin=reCast<IssmPDouble>(dxmind);
 	femmodel->parameters->FindParam(&gttold,InversionGttolEnum); gttol=reCast<IssmPDouble>(gttold);
-	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");
@@ -70,29 +67,13 @@
 	long niter = long(maxsteps); /*Maximum number of iterations*/
 	long nsim  = long(maxiter);/*Maximum number of function calls*/
 
-	/*Get initial guess*/
-	Vector<IssmDouble> *Xad = NULL;
-	GetVectorFromControlInputsx(&Xad,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value",onsid);
-	Xd = Xad->ToMPISerial();
-	Xad->GetSize(&intn);
-	X=xNew<IssmPDouble>(intn);
-	for(i=0;i<intn;i++) X[i]=reCast<IssmPDouble>(Xd[i]);
-	delete Xad;
-	_assert_(intn==numberofvertices*num_controls);
-	
-	/*Get problem dimension and initialize gradient and initial guess*/
+	/*Run the first part of simulad, in order to get things started!:*/
+	femmodelad=presimulad(&intn,&X,femmodel);
+
+	/*Get problem dimension and initialize gradient: */
 	long n = long(intn);
 	G = xNew<IssmPDouble>(n);
-	Gd = xNew<IssmDouble>(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 doc)*/
 	long      m   = 100;
 	long      ndz = 4*n+m*(2*n+1);
@@ -103,8 +84,8 @@
 	_printf0_("Cost function f(x)   | Gradient norm |g(x)| |  List of contributions\n");
 	_printf0_("____________________________________________________________________\n");
 
-	//first run before firing up the control optimization
-	simulad(&indic,&n,X,&f,G,izs,rzs,(void*)femmodel);
+	//run post simular phase, to fire up the control optimization
+	postsimulad(&indic,&n,X,&f,G,izs,rzs,(void*)femmodelad); 
 	double f1=f;
 
 	m1qn3_(costfuncion,prosca,&ctonbe_,&ctcabe_,
@@ -124,20 +105,6 @@
 		default: _printf0_("   Unknown end condition\n");
 	}
 	
-	/*Constrain solution vector*/
-	IssmDouble  *XL = NULL;
-	IssmDouble  *XU = NULL;
-	GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound",onsid);
-	GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound",onsid);
-	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]);
-			if(X[index]>XU[index]) X[index]=reCast<IssmPDouble>(XU[index]);
-			if(X[index]<XL[index]) X[index]=reCast<IssmPDouble>(XL[index]);
-		}
-	}
-
 	/*Save results:*/
 	femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,G,n,1,1,0.0));
 	femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffXpEnum,X,intn,1,1,0.0));
@@ -146,8 +113,125 @@
 	xDelete<double>(G);
 	xDelete<double>(X);
 	xDelete<double>(dz);
-	xDelete<IssmDouble>(scaling_factors);
 }/*}}}*/
+FemModel* presimulad(int* pintn, double** pX, FemModel* femmodel){ /*{{{*/
+
+	/*Intermediaries:*/
+	char* rootpath=NULL;
+	char* inputfilename=NULL;
+	char* outputfilename=NULL;
+	char* toolkitsfilename=NULL;
+	char* lockfilename=NULL;
+	int         solution_type;
+	IssmDouble    pfd;
+	IssmDouble*   Xd=NULL;
+	int           intn;
+	IssmPDouble*   X=NULL;
+	int            i;
+	
+	/*Now things get complicated. The femmodel we recovered did not initialize an AD trace, so we can't compute gradients with it. We are going to recreate 
+	 *a new femmodel, identical in all aspects to the first one, with trace on though, which will allow us to run the forward mode and get the gradient 
+	 in one run of the solution core. So first recover the filenames required for the FemModel constructor, then call a new ad tailored constructor:*/
+	femmodel->parameters->FindParam(&rootpath,RootPathEnum);
+	femmodel->parameters->FindParam(&inputfilename,InputFileNameEnum);
+	femmodel->parameters->FindParam(&outputfilename,OutputFileNameEnum);
+	femmodel->parameters->FindParam(&toolkitsfilename,ToolkitsFileNameEnum);
+	femmodel->parameters->FindParam(&lockfilename,LockFileNameEnum);
+
+	femmodel=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, femmodel->comm, femmodel->solution_type,NULL);
+
+	
+	/*Get initial guess:*/
+	femmodel->parameters->FindParam(&Xd,&intn,AutodiffXpEnum);
+	X=xNew<IssmPDouble>(intn); for(i=0;i<intn;i++) X[i]=reCast<IssmPDouble>(Xd[i]); 
+
+	xDelete<char>(rootpath);
+	xDelete<char>(inputfilename);
+	xDelete<char>(outputfilename);
+	xDelete<char>(toolkitsfilename);
+	xDelete<char>(lockfilename);
+	xDelete<IssmDouble>(Xd);
+
+	*pintn=intn;
+	*pX=X;
+
+	return femmodel;
+
+} /*}}}*/
+void postsimulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){ /*{{{*/
+
+	/*Intermediaries:*/
+	char* rootpath=NULL;
+	char* inputfilename=NULL;
+	char* outputfilename=NULL;
+	char* toolkitsfilename=NULL;
+	char* lockfilename=NULL;
+	IssmPDouble* G2=NULL;
+	int         solution_type;
+	FemModel   *femmodel  =  NULL;
+	IssmDouble    pfd;
+	int            i;
+	
+	/*Recover Femmodel*/
+	femmodel  = (FemModel*)dzs;
+
+	/*Recover number of cost functions responses*/
+	int num_responses;
+	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+
+	/*Recover some parameters*/
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+
+	/*Compute solution:*/
+	void (*solutioncore)(FemModel*)=NULL;
+	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
+	solutioncore(femmodel);
+
+	/*Compute objective function*/
+	IssmDouble* Jlist = NULL;
+	femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast<IssmPDouble>(pfd);
+	_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
+	
+	/*Compute gradient using AD. Gradient is in the results after the ad_core is called*/
+	adgradient_core(femmodel); 
+
+	if(IssmComm::GetRank()==0){
+		IssmPDouble* G_temp=NULL;
+		GenericExternalResult<IssmPDouble*>* gradient=(GenericExternalResult<IssmPDouble*>*)femmodel->results->FindResult(AutodiffJacobianEnum); _assert_(gradient);
+		G_temp=gradient->GetValues();
+		/*copy onto G2, to avoid a leak: */
+		G2=xNew<IssmPDouble>(*n); 
+		xMemCpy<IssmPDouble>(G2,G_temp,*n);
+	}
+	else G2=xNew<IssmPDouble>(*n);
+
+	/*MPI broadcast results:*/
+	ISSM_MPI_Bcast(G2,*n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+	
+	/*Send gradient to m1qn3 core: */
+	for(long i=0;i<*n;i++) G[i] = G2[i];
+	
+	/*Constrain X and G*/
+	IssmDouble  Gnorm = 0.;
+	for(long i=0;i<*n;i++) Gnorm += G[i]*G[i];
+	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[i]);
+	_printf0_("\n");
+
+	/*Clean-up and return*/
+	xDelete<IssmDouble>(Jlist);
+	xDelete<IssmPDouble>(G2);
+	
+	xDelete<char>(rootpath);
+	xDelete<char>(inputfilename);
+	xDelete<char>(outputfilename);
+	xDelete<char>(toolkitsfilename);
+	xDelete<char>(lockfilename);
+
+} /*}}}*/
 void simulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){ /*{{{*/
 
 	/*Intermediaries:*/
@@ -157,9 +241,6 @@
 	char* toolkitsfilename=NULL;
 	char* lockfilename=NULL;
 	IssmPDouble* G2=NULL;
-	bool onsid=true;
-	IssmDouble  *XL = NULL;
-	IssmDouble  *XU = NULL;
 	int         solution_type;
 	FemModel   *femmodel  =  NULL;
 	FemModel   *femmodelad  = NULL;
@@ -170,25 +251,9 @@
 	femmodel  = (FemModel*)dzs;
 
 	/*Recover number of cost functions responses*/
-	int num_responses,num_controls,numberofvertices;
-	IssmDouble* scaling_factors = NULL;
+	int num_responses;
 	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
-	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
-	femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
-	numberofvertices=femmodel->vertices->NumberOfVertices();
 
-	/*Constrain input vector*/
-	GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound",onsid);
-	GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound",onsid);
-	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]);
-			if(X[index]>reCast<IssmPDouble>(XU[index])) X[index]=reCast<IssmPDouble>(XU[index]);
-			if(X[index]<reCast<IssmPDouble>(XL[index])) X[index]=reCast<IssmPDouble>(XL[index]);
-		}
-	}
-
 	/*Now things get complicated. The femmodel we recovered did not initialize an AD trace, so we can't compute gradients with it. We are going to recreate 
 	 *a new femmodel, identical in all aspects to the first one, with trace on though, which will allow us to run the forward mode and get the gradient 
 	 in one run of the solution core. So first recover the filenames required for the FemModel constructor, then call a new ad tailored constructor:*/
@@ -233,18 +298,9 @@
 	/*Send gradient to m1qn3 core: */
 	for(long i=0;i<*n;i++) G[i] = G2[i];
 	
-	/*Constrain X and G*/
+	/*Recover Gnorm: */
 	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<IssmPDouble>(scaling_factors[c]);
-			X[index] = X[index]/reCast<IssmPDouble>(scaling_factors[c]);
-			Gnorm += G[index]*G[index];
-		}
-	}
+	for(int i=0;i<*n;i++) Gnorm += G[i]*G[i];
 	Gnorm = sqrt(Gnorm);
 
 	/*Print info*/
@@ -254,10 +310,7 @@
 
 	/*Clean-up and return*/
 	xDelete<IssmDouble>(Jlist);
-	xDelete<IssmDouble>(XU);
-	xDelete<IssmDouble>(XL);
 	xDelete<IssmPDouble>(G2);
-	xDelete<IssmDouble>(scaling_factors);
 	
 	xDelete<char>(rootpath);
 	xDelete<char>(inputfilename);
