Index: ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 18713)
+++ ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 18714)
@@ -31,7 +31,7 @@
 	long         omode;
 	double       f,dxmin,gttol; 
 	int          maxsteps,maxiter;
-	int          intn,num_controls,solution_type;
+	int          intn,numberofvertices,num_controls,solution_type;
 	IssmDouble  *scaling_factors = NULL;
 	IssmDouble  *X  = NULL;
 	IssmDouble  *G  = NULL;
@@ -45,6 +45,7 @@
 	femmodel->parameters->FindParam(&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");
@@ -61,8 +62,8 @@
 	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*/
+	long niter = long(maxsteps); /*Maximum number of iterations*/
+	long nsim  = long(maxiter);/*Maximum number of function calls*/
 
 	/*Get initial guess*/
 	Vector<IssmDouble> *Xpetsc = NULL;
@@ -70,15 +71,18 @@
 	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*/
-	if(num_controls!=1) _error_("not supported yet...");
-	for(long i=0;i<n;i++){
-		X[i] = X[i]/scaling_factors[0];
+	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];
+		}
 	}
 
 	/*Allocate m1qn3 working arrays (see doc)*/
@@ -117,11 +121,13 @@
 	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");
-	if(num_controls!=1) _error_("not supported yet");
-	for(long i=0;i<n;i++){
-		X[i] = X[i]*scaling_factors[0];
-		if(X[i]>XU[i]) X[i]=XU[i];
-		if(X[i]<XL[i]) X[i]=XL[i];
+	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);
 	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
@@ -151,23 +157,25 @@
 	FemModel   *femmodel  = (FemModel*)dzs;
 
 	/*Recover number of cost functions responses*/
-	int         num_responses;
-	int         num_controls;
+	int num_responses,num_controls,numberofvertices;
 	IssmDouble* scaling_factors = NULL;
 	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*/
 	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");
-	if(num_controls!=1) _error_("not supported yet");
-	for(long i=0;i<*n;i++){
-		X[i] = X[i]*scaling_factors[0];
-		if(X[i]>XU[i]) X[i]=XU[i];
-		if(X[i]<XL[i]) X[i]=XL[i];
+	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];
+		}
 	}
 
 	/*Update control input*/
@@ -213,13 +221,15 @@
 
 	/*Constrain Gradient*/
 	IssmDouble  Gnorm = 0.;
-	if(num_controls!=1) _error_("not supported yet");
-	for(long i=0;i<*n;i++){
-		if(X[i]>=XU[i]) G[i]=0.;
-		if(X[i]<=XL[i]) G[i]=0.;
-		G[i] = G[i]*scaling_factors[0];
-		X[i] = X[i]/scaling_factors[0];
-		Gnorm += G[i]*G[i];
+	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);
 
