Index: /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 18713)
+++ /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 18714)
@@ -32,5 +32,5 @@
 	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;
@@ -46,4 +46,5 @@
 	femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
 	femmodel->parameters->SetParam(false,SaveResultsEnum);
+	numberofvertices=femmodel->vertices->NumberOfVertices();
 
 	/*Initialize M1QN3 parameters*/
@@ -62,6 +63,6 @@
 
 	/*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*/
@@ -71,4 +72,5 @@
 	Xpetsc->GetSize(&intn);
 	delete Xpetsc;
+	_assert_(intn==numberofvertices*num_controls);
 
 	/*Get problem dimension and initialize gradient and initial guess*/
@@ -77,7 +79,9 @@
 
 	/*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];
+		}
 	}
 
@@ -118,9 +122,11 @@
 	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);
@@ -152,10 +158,10 @@
 
 	/*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*/
@@ -164,9 +170,11 @@
 	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];
+		}
 	}
 
@@ -214,11 +222,13 @@
 	/*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);