Index: /issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp	(revision 13313)
+++ /issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp	(revision 13314)
@@ -12,49 +12,29 @@
 
 void AutodiffDriversx(Elements* elements,Nodes* nodes,Vertices* vertices,Loads* loads,Materials* materials,Parameters* parameters,Results* results){
-
-	int  i                  ,j;
 	bool isautodiff       = false;
-	int  num_dependents;
-	int  num_independents;
-	int  numberofvertices;
-
-	/*outputs: */
-	IssmDouble  *matJ = NULL;
-	IssmDouble  *axp  = NULL;
-	double     **J    = NULL;
-	double      *xp   = NULL;
-		
-	/*diverse: */
-	int configuration_type;
-	int solveSize;
-	int edf_n,edf_m;
-
 	/*AD mode on?: */
 	parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
-
 	if(isautodiff){
-		#ifdef _HAVE_ADOLC_
-		parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
-		parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
-		parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum);
-
-		if(!(num_dependents*num_independents)) return;
-		
+                #ifdef _HAVE_ADOLC_
+	        int  num_dependents;
+                parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
+	        int  num_independents;
+                parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
+                if(!(num_dependents*num_independents)) return;
+	        int  numberofvertices;
+                parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum);
 		/*retrieve state variable: */
+                IssmDouble  *axp  = NULL;
 		parameters->FindParam(&axp,&num_independents,AutodiffXpEnum);
-
-		/*allocate driver results: */
-		J=xNew<double*>(num_independents);
-		xp=xNew<double>(num_independents); 
-
-		for(i=0;i<num_independents;i++){
-			J[i]=xNew<double>(num_dependents);
+		/* driver argument */
+		double *xp=xNew<double>(num_independents);
+		for(int i=0;i<num_independents;i++){
 			xp[i]=reCast<double,IssmDouble>(axp[i]);
 		}
-
 		/* get the dimension for the solverx arguments*/
+                int configuration_type;
 		parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-		solveSize=nodes->NumberOfDofs(configuration_type,FsetEnum);
-		edf_n=solveSize*(solveSize+1), edf_m=solveSize;
+		int solveSize=nodes->NumberOfDofs(configuration_type,FsetEnum);
+		int edf_n=solveSize*(solveSize+1), edf_m=solveSize;
 
 		/*get the EDF pointer:*/
@@ -62,7 +42,9 @@
 		
 		/*set the forward method function pointers: */
+		anEDF_for_solverx_p->zos_forward=EDF_for_solverx;
 		anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx;
 		anEDF_for_solverx_p->fov_forward=EDF_fov_forward_for_solverx;
-
+                anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
+                // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
 		/*allocate the space for the parameters to invoke the forward methods:*/
 		anEDF_for_solverx_p->dp_x=xNew<double>(edf_n);
@@ -72,31 +54,47 @@
 		anEDF_for_solverx_p->dp_Y=xNew<double>(edf_m);
 		anEDF_for_solverx_p->dpp_Y=xNew<double>(edf_m, num_independents);
+		anEDF_for_solverx_p->dp_U=xNew<double>(edf_m);
+		anEDF_for_solverx_p->dpp_U=xNew<double>(num_dependents,edf_m);
+                anEDF_for_solverx_p->dp_Z=xNew<double>(edf_n);
+                anEDF_for_solverx_p->dpp_Z=xNew<double>(num_dependents,edf_n);
+		/* Call AD driver:*/
+                // single direction:
+                double *tangentDir=xNewZeroInit<double>(num_independents);
+                unsigned int anIndepNum=55; // <---------- make this selectable via config
+                tangentDir[anIndepNum]=1.0;
+                double *theJacVecProduct=xNew<double>(num_dependents);
+                double *theOutput=xNew<double>(num_dependents);
+                if (fos_forward(1,num_dependents,num_independents, 0, xp, tangentDir, theOutput, theJacVecProduct ))
+                  _error_("fos_forward returned non-zero error code");
+                IssmDouble *theJacVecProductA=xNew<IssmDouble>(num_dependents);
+                for (int i=0;i<num_dependents;++i)
+                  theJacVecProductA[i]=theJacVecProduct[i];
+                results->AddObject(new DoubleVecExternalResult(results->Size()+1,AutodiffJacobianEnum,theJacVecProductA,num_dependents,1,0.0));
 
-		/* Call AD driver:*/
+#if 0
+                double **J=xNew<double>(num_dependents,num_independents);
 		jacobian(1,num_dependents,num_independents,xp,J);
+		IssmDouble *matJ=xNew<IssmDouble>(num_dependents*num_independents);
+                for (int i=0;i<num_dependents*num_independents;++i)
+                  matJ[i]=J[i];
+                results->AddObject(new DoubleMatExternalResult(results->Size()+1,AutodiffJacobianEnum,matJ,num_independents,num_dependents,1,0.0));
+                xDelete(matJ);
+                xDelete(J);
+#endif
 
-		/* delete the allocated space for the parameters:*/
-		xDelete(anEDF_for_solverx_p->dp_x);
-		xDelete(anEDF_for_solverx_p->dp_X);
-		xDelete(anEDF_for_solverx_p->dpp_X,edf_n);
-		xDelete(anEDF_for_solverx_p->dp_y);
-		xDelete(anEDF_for_solverx_p->dp_Y);
-		xDelete(anEDF_for_solverx_p->dpp_Y,edf_m);
-		
-		/*Add jacobian matrix to results: */
-		matJ=xNew<IssmDouble>(num_dependents*num_independents);
-		for (i=0;i<num_independents;++i){
-			double* rowi=J[i];
-			for (j=0;j<num_dependents;++j){
-				*(matJ+num_dependents*i+j)=rowi[j];
-			}
-		}
-		results->AddObject(new DoubleMatExternalResult(results->Size()+1,AutodiffJacobianEnum,matJ,num_independents,num_dependents,1,0.0));
-
+                /* delete the allocated space for the parameters:*/
+                xDelete(anEDF_for_solverx_p->dp_x);
+                xDelete(anEDF_for_solverx_p->dp_X);
+                xDelete(anEDF_for_solverx_p->dpp_X);
+                xDelete(anEDF_for_solverx_p->dp_y);
+                xDelete(anEDF_for_solverx_p->dp_Y);
+                xDelete(anEDF_for_solverx_p->dpp_Y);
+                xDelete(anEDF_for_solverx_p->dp_U);
+                xDelete(anEDF_for_solverx_p->dpp_U);
+                xDelete(anEDF_for_solverx_p->dp_Z);
+                xDelete(anEDF_for_solverx_p->dpp_Z);
 		/*Free ressources: */
-		xDelete<double*>(J);
-		xDelete<double>(xp);
-		xDelete<IssmDouble>(matJ);
-		xDelete<IssmDouble>(axp);
+		xDelete(xp);
+		xDelete(axp); // did we allocate this?
 		#else
 		_error_("Should not be requesting AD drivers when an AD library is not available!");