Index: issm/trunk/src/c/Makefile.am
===================================================================
--- issm/trunk/src/c/Makefile.am	(revision 8277)
+++ issm/trunk/src/c/Makefile.am	(revision 8279)
@@ -1,3 +1,3 @@
-INCLUDES = @DAKOTAINCL@ @PETSCINCL@ @SLEPCINCL@ @MPIINCL@ @MATLABINCL@  @METISINCL@  @CHACOINCL@ @SCOTCHINCL@ @PLAPACKINCL@  @BLASLAPACKINCL@ @MKLINCL@ @MUMPSINCL@  @TRIANGLEINCL@ @HYPREINCL@ @MLINCL@ 
+INCLUDES = @DAKOTAINCL@ @PETSCINCL@ @SLEPCINCL@ @MPIINCL@ @MATLABINCL@  @METISINCL@  @CHACOINCL@ @SCOTCHINCL@ @PLAPACKINCL@  @BLASLAPACKINCL@ @MKLINCL@ @MUMPSINCL@  @TRIANGLEINCL@ @HYPREINCL@ @MLINCL@ @TAOINCL@
 
 #Compile serial library, and then try and compile parallel library
@@ -375,4 +375,5 @@
 					./toolkits/petsc/patches/MatMultPatch.cpp\
 					./toolkits/petsc/petscincludes.h\
+					./toolkits/tao/taoincludes.h\
 					./toolkits/mpi/mpiincludes.h\
 					./toolkits/mpi/patches/mpipatches.h\
@@ -1015,4 +1016,5 @@
 					./toolkits/petsc/patches/MatMultPatch.cpp\
 					./toolkits/petsc/petscincludes.h\
+					./toolkits/tao/taoincludes.h\
 					./toolkits/mpi/mpiincludes.h\
 					./toolkits/mpi/patches/mpipatches.h\
@@ -1285,4 +1287,5 @@
 					./solutions/WriteLockFile.cpp\
 					./solutions/control_core.cpp\
+					./solutions/controltao_core.cpp\
 					./solutions/controlrestart.cpp\
 					./solutions/controlconvergence.cpp\
@@ -1325,5 +1328,5 @@
 endif
 
-LDADD =    ./libpISSM.a $(PETSCLIB)     $(FLIBS)  $(PLAPACKLIB)  $(MUMPSLIB) $(SCALAPACKLIB)  $(BLACSLIB) $(HYPRELIB) $(MLLIB)   $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB)  $(MKLLIB) $(MPILIB) $(MATHLIB) $(FORTRANLIB) $(GRAPHICSLIB) libOverload.a $(MULTITHREADINGLIB)
+LDADD =    ./libpISSM.a $(PETSCLIB) $(TAOLIB) $(FLIBS) $(PLAPACKLIB)  $(MUMPSLIB) $(SCALAPACKLIB)  $(BLACSLIB) $(HYPRELIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB)  $(MKLLIB) $(MPILIB) $(MATHLIB) $(FORTRANLIB) $(GRAPHICSLIB) libOverload.a $(MULTITHREADINGLIB)
 
 issm_exe_SOURCES = solutions/issm.cpp
Index: issm/trunk/src/c/shared/Matlab/ModuleBoot.cpp
===================================================================
--- issm/trunk/src/c/shared/Matlab/ModuleBoot.cpp	(revision 8277)
+++ issm/trunk/src/c/shared/Matlab/ModuleBoot.cpp	(revision 8279)
@@ -24,4 +24,9 @@
 	PetscInitializeNoArguments();
 
+	/*Initialize Tao*/
+	#ifdef _HAVE_TAO_
+	TaoInitialize(0,0,0,0);
+	#endif
+
 	return 1;
 }
Index: issm/trunk/src/c/solutions/controltao_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/controltao_core.cpp	(revision 8279)
+++ issm/trunk/src/c/solutions/controltao_core.cpp	(revision 8279)
@@ -0,0 +1,143 @@
+/*!\file: control_core.cpp
+ * \brief: core of the control solution 
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "./solutions.h"
+#include "../modules/modules.h"
+#include "../include/include.h"
+#include "../solvers/solvers.h"
+
+
+void controltao_core(FemModel* femmodel){
+
+	int     i,n;
+	
+	/*parameters: */
+	int     num_controls;
+	int     nsteps;
+	double  eps_cm;
+	double  tolx;
+	bool    cm_gradient;
+	int     dim;
+	int     solution_type;
+	bool    isstokes;
+	bool    qmu_analysis=false;
+
+	int*    control_type = NULL;
+	double* responses=NULL;
+	double* maxiter=NULL;
+	double* cm_jump=NULL;
+
+		
+	/*intermediary: */
+	double  search_scalar=0;
+	OptArgs optargs;
+	OptPars optpars;
+
+	/*Solution and Adjoint core pointer*/
+	void (*solutioncore)(FemModel*)=NULL;
+	void (*adjointcore)(FemModel*)=NULL;
+
+	/*output: */
+	double* J=NULL;
+
+	printf("I am in TAO!!\n");
+
+	/*Recover parameters used throughout the solution:{{{1*/
+	femmodel->parameters->FindParam(&nsteps,NStepsEnum);
+	femmodel->parameters->FindParam(&num_controls,NumControlsEnum);
+	femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
+	femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum);
+	femmodel->parameters->FindParam(&maxiter,NULL,MaxIterEnum);
+	femmodel->parameters->FindParam(&cm_jump,NULL,CmJumpEnum);
+	femmodel->parameters->FindParam(&eps_cm,EpsCmEnum);
+	femmodel->parameters->FindParam(&tolx,TolXEnum);
+	femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum);
+	femmodel->parameters->FindParam(&dim,DimEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
+	femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
+	/*}}}*/
+
+	/*out of solution_type, figure out solution core and adjoint function pointer*/
+	CorePointerFromSolutionEnum(&solutioncore,solution_type);
+	AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
+
+	/*Launch once a complete solution to set up all inputs*/
+	_printf_(VerboseControl(),"%s\n","   preparing initial solution");
+	if (isstokes) solutioncore(femmodel);
+
+	/*Initialize responses: */
+	J=(double*)xmalloc(nsteps*sizeof(double));
+		
+	/*Initialize some of the BrentSearch arguments: */
+	optargs.femmodel=femmodel;
+	optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx ;
+	
+	/*Start looping: */
+	for(n=0;n<nsteps;n++){
+
+		/*Display info*/
+		_printf_(VerboseControl(),"\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
+		InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)responses[n],CmResponseEnum);
+		
+		/*In case we are running a steady state control method, compute new temperature field using new parameter distribution: */
+		if (solution_type==SteadystateSolutionEnum) solutioncore(femmodel);
+
+		_printf_(VerboseControl(),"%s\n","   compute adjoint state:");
+		adjointcore(femmodel);
+	
+		gradient_core(femmodel,n,search_scalar);
+
+		/*Return gradient if asked: */
+		if (cm_gradient){
+			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GradientEnum);
+			goto cleanup_and_return;
+		}
+
+		_printf_(VerboseControl(),"%s\n","   optimizing along gradient direction");
+		optargs.n=n; optpars.maxiter=(int)maxiter[n]; optpars.cm_jump=cm_jump[n];
+		BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
+		//OptimalSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
+
+		_printf_(VerboseControl(),"%s\n","   updating parameter using optimized search scalar"); //true means update save controls
+		InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true);
+		
+		if(controlconvergence(J,responses,eps_cm,n)) break;
+
+		/*Temporary saving every 5 control steps: */
+		/*if (((n+1)%5)==0){
+			_printf_(VerboseControl(),"%s\n","   saving temporary results");
+			controlrestart(femmodel,J);
+		}*/
+	}
+
+	_printf_(VerboseControl(),"%s\n","   preparing final solution");
+	femmodel->parameters->SetParam(false,ControlAnalysisEnum); //needed to turn control result output in solutioncore
+	solutioncore(femmodel);
+
+	/*some results not computed by steadystate_core or diagnostic_core: */
+	if(!qmu_analysis){ //do not save this if we are running the control core from a qmu run!
+		for(i=0;i<num_controls;i++) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]);
+		femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0));
+		//femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumToStringx(control_type),1,0));
+	}
+
+	cleanup_and_return:
+	/*Free ressources: */
+	xfree((void**)&control_type);
+	xfree((void**)&responses);
+	xfree((void**)&maxiter);
+	xfree((void**)&cm_jump);
+	xfree((void**)&J);
+	
+	/*control_core might be used in Qmu, so leave everything similar to where it started: */
+	femmodel->parameters->SetParam(true,ControlAnalysisEnum);
+
+	/* Finalize TAO */
+	TaoFinalize();
+}
Index: issm/trunk/src/c/solutions/issm.cpp
===================================================================
--- issm/trunk/src/c/solutions/issm.cpp	(revision 8277)
+++ issm/trunk/src/c/solutions/issm.cpp	(revision 8279)
@@ -79,5 +79,11 @@
 
 	/*if control is being run on top of a solution, change core: */
-	if(control_analysis)solutioncore=&control_core;
+	if(control_analysis){
+		#ifdef _HAVE_TAO_
+		solutioncore=&controltao_core;
+		#else
+		solutioncore=&control_core;
+		#endif
+	}
 
 	/*are we running the solution sequence, or a qmu wrapper around it? : */
Index: issm/trunk/src/c/solutions/solutions.h
===================================================================
--- issm/trunk/src/c/solutions/solutions.h	(revision 8277)
+++ issm/trunk/src/c/solutions/solutions.h	(revision 8279)
@@ -24,4 +24,5 @@
 void bedslope_core(FemModel* femmodel);
 void control_core(FemModel* femmodel);
+void controltao_core(FemModel* femmodel);
 void prognostic_core(FemModel* femmodel);
 void balancedthickness_core(FemModel* femmodel);
Index: issm/trunk/src/c/toolkits/tao/taoincludes.h
===================================================================
--- issm/trunk/src/c/toolkits/tao/taoincludes.h	(revision 8279)
+++ issm/trunk/src/c/toolkits/tao/taoincludes.h	(revision 8279)
@@ -0,0 +1,17 @@
+/* \file toolkits.h
+ * \brief: this API allows use of external packages, provides patches, etc ...
+ */
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#ifndef _TAOINCLUDES_H_
+#define _TOOINCLUDES_H_
+
+#ifdef _HAVE_TAO_
+#include "tao.h"
+#endif
+
+#endif
Index: issm/trunk/src/c/toolkits/toolkits.h
===================================================================
--- issm/trunk/src/c/toolkits/toolkits.h	(revision 8277)
+++ issm/trunk/src/c/toolkits/toolkits.h	(revision 8279)
@@ -10,4 +10,5 @@
 #include "./metis/metisincludes.h"
 #include "./triangle/triangleincludes.h"
+#include "./tao/taoincludes.h"
 
 #endif
