Index: /issm/trunk-jpl/src/c/classes/matrix/Vector.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/Vector.cpp	(revision 12909)
+++ /issm/trunk-jpl/src/c/classes/matrix/Vector.cpp	(revision 12910)
@@ -61,4 +61,15 @@
 }
 /*}}}*/
+#ifdef _HAVE_PETSC_
+/*FUNCTION Vector::Vector(Vec petsc_vector){{{*/
+Vector::Vector(Vec petsc_vector){
+
+	this->type=PetscVecType;
+	this->svector=NULL;
+	this->pvector=new PetscVec(petsc_vector);
+
+}
+/*}}}*/
+#endif
 /*FUNCTION Vector::Vector(IssmDouble* serial_vec,int M,int type){{{*/
 Vector::Vector(IssmDouble* serial_vec,int M,int in_type){
Index: /issm/trunk-jpl/src/c/classes/matrix/Vector.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/Vector.h	(revision 12909)
+++ /issm/trunk-jpl/src/c/classes/matrix/Vector.h	(revision 12910)
@@ -36,4 +36,7 @@
 		Vector(int M,bool fromlocalsize=false,int type=PetscVecType);
 		Vector(IssmDouble* serial_vec,int pM,int type=PetscVecType);
+		#ifdef _HAVE_PETSC_
+		Vector(Vec petsc_vector);
+		#endif
 		~Vector();
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/solutions/controltao_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/controltao_core.cpp	(revision 12909)
+++ /issm/trunk-jpl/src/c/solutions/controltao_core.cpp	(revision 12910)
@@ -21,4 +21,5 @@
 typedef struct {
 	FemModel* femmodel;
+	double*   J;
 } AppCtx;
 
@@ -69,9 +70,10 @@
 	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");
-	TaoSetInitialVector(tao,X->vector);
-	TaoSetVariableBounds(tao,XL->vector,XU->vector);
+	TaoSetInitialVector(tao,X->pvector->vector);
+	TaoSetVariableBounds(tao,XL->pvector->vector,XU->pvector->vector);
 	xdelete(&XL);
 	xdelete(&XU);
 
+	user.J=(double*)xcalloc((maxiter+5),sizeof(double));
 	user.femmodel=femmodel;
 	TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user); 
@@ -81,9 +83,10 @@
 	TaoSolve(tao);
 	TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
-	TaoGetSolutionVector(tao,&X->vector);
+	TaoGetSolutionVector(tao,&X->pvector->vector);
 	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
 	for(int i=0;i<num_controls;i++){
 		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_list[i]);
 	}
+	femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,user.J,maxiter+3,1,0));
 
 	/*Finalize*/
@@ -142,5 +145,5 @@
 	/*Compute gradient*/
 	Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-	VecCopy(gradient->vector,G); xdelete(&gradient);
+	VecCopy(gradient->pvector->vector,G); xdelete(&gradient);
 	VecScale(G,-1.);
 
@@ -166,4 +169,5 @@
 	if(its==0) _pprintLine_("-----------------------------------+-----------------------");
 	_pprintString_(setw(4)<<its<<"   "<<setw(12)<<setprecision(7)<<f<<"  "<<setw(12)<<setprecision(7)<<gnorm<<"  | ");
+	user->J[its]=f;
 
 	/*Retrieve objective functions independently*/