Index: /issm/trunk-jpl/src/c/solutions/controltao_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/controltao_core.cpp	(revision 11296)
+++ /issm/trunk-jpl/src/c/solutions/controltao_core.cpp	(revision 11297)
@@ -25,54 +25,49 @@
 
 	/*TAO*/
-	int       i,n,info;
+	int       ierr,numberofvertices;
+	AppCtx    user;
 	TaoSolver tao;
-	Vec       initial_solution = NULL;
-	AppCtx    user;
+	Vec       X  = NULL;
+	Vec       XL = NULL;
+	Vec       XU = NULL;
 
 	/*Initialize TAO*/
-	int argc; char **args; PetscGetArgs(&argc,&args);
-	int ierr=TaoInitialize(&argc,&args,(char*)0,"");
+	int argc; char **args=NULL;
+	PetscGetArgs(&argc,&args);
+	ierr = TaoInitialize(&argc,&args,(char*)0,"");
 	if(ierr) _error_("Could not initialize Tao");
 
-	info = PetscOptionsSetValue("-tao_monitor",""); if(info) _error_("STOP");
+	/*Initialize TAO*/
+	TaoCreate(PETSC_COMM_WORLD,&tao);
+	PetscOptionsSetValue("-tao_monitor","");
+	TaoSetFromOptions(tao);
+	TaoSetType(tao,"tao_blmvm");
 
-	/*Initialize argument*/
-	user.femmodel=femmodel;
+	/*Prepare all TAO parameters*/
+	TaoSetMaximumFunctionEvaluations(tao,50);
+	TaoSetMaximumIterations(tao,10);
+	TaoSetTolerances(tao,10e-28,0.0000,0.0000,0.0000,10e-28);
 
-	/*Set up and solve TAO*/
-	info = TaoCreate(PETSC_COMM_WORLD,&tao); if(info) _error_("STOP");
-	info = TaoSetType(tao,"tao_blmvm"); if(info) _error_("STOP");
-	//info = TaoSetType(tao,"tao_cg"); if(info) _error_("STOP");
-	TaoSetMaximumFunctionEvaluations(tao,10);
-	TaoSetMaximumIterations(tao,3);
-
-	int size=femmodel->vertices->NumberOfVertices();
-	Vec XL=NULL;
-	Vec XU=NULL;
-	VecCreate(PETSC_COMM_WORLD,&XL);                               VecCreate(PETSC_COMM_WORLD,&XU);
-	VecSetSizes(XL,PetscDetermineLocalSize(size),PETSC_DECIDE);    VecSetSizes(XU,PetscDetermineLocalSize(size),PETSC_DECIDE);
-	VecSetFromOptions(XL);                                         VecSetFromOptions(XU);
-	VecSet(XL,1);                                                  VecSet(XU,200.);
+	numberofvertices = femmodel->vertices->NumberOfVertices();
+	XL=NewVec(numberofvertices); VecSet(XL,1.);
+	XU=NewVec(numberofvertices); VecSet(XU,200.);
 	TaoSetVariableBounds(tao,XL,XU); VecFree(&XL); VecFree(&XU);
 
-	GetVectorFromInputsx(&initial_solution,femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,FrictionCoefficientEnum,VertexEnum);
-	info = TaoSetInitialVector(tao,initial_solution);  if(info) _error_("STOP");
-	VecFree(&initial_solution);
+	GetVectorFromInputsx(&X,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum,VertexEnum);
+	TaoSetInitialVector(tao,X);
 
-	info = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);  if(info) _error_("STOP");
-	info = TaoSetFromOptions(tao);  if(info) _error_("STOP");
-	/* http://www.mcs.anl.gov/research/projects/tao/docs/manpages/taosolver/TaoSetTolerances.html*/
-	/*                          fatol ,frtol ,gatol ,grtol ,gttol*/
-	info = TaoSetTolerances(tao,10e-28,0.0000,0.0000,0.0000,10e-28); if(info) _error_("STOP");
-	info = TaoSolve(tao); if(info) _error_("STOP");
+	user.femmodel=femmodel;
+	TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user); 
 
-	/*Get solution status*/
-	info = TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);  if(info) _error_("STOP");
-
-	/*Clean up*/
-	info = TaoDestroy(&tao);  if(info) _error_("STOP");
+	/*Solver optimization problem*/
+	TaoSolve(tao);
+	TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
+	TaoGetSolutionVector(tao,&X);
+	InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum);
 	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum);
 
-	/* Finalize TAO */
+	/*Clean up and return*/
+	VecFree(&X);
+	TaoDestroy(&tao);
 	TaoFinalize();
 }
@@ -86,14 +81,13 @@
 	int costfunction=SurfaceAbsVelMisfitEnum;
 	femmodel->parameters->SetParam(&costfunction,1,1,StepResponsesEnum);
-	//InputUpdateFromVectorx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,X,MaterialsRheologyBbarEnum,VertexEnum);
+	//InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,MaterialsRheologyBbarEnum,VertexEnum);
 	InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum);
 	adjointdiagnostic_core(user->femmodel);
-	//Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, MaterialsRheologyBbarEnum);
-	Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,FrictionCoefficientEnum);
+	//Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaterialsRheologyBbarEnum);
+	Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum);
 	VecScale(gradient,-1.);
 	VecCopy(gradient,G);
 	VecFree(&gradient);
-	CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
-	_printf_(true,"f(x)=%g\n",*fcn);
+	CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 	return 0;
 }