Index: /issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp	(revision 13451)
+++ /issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp	(revision 13452)
@@ -2,5 +2,5 @@
  * \brief: compute outputs from the AD mode,  using our dependents and independents, and drivers available in Adolc.
  */
-
+#include <set>
 #include "../../modules/modules.h"
 #include "../../shared/shared.h"
@@ -21,5 +21,7 @@
 	IssmDouble *axp                = NULL;
 	double     *xp                 = NULL;
-	int         anIndepNum;
+	int         anIndepIndex;
+	unsigned int * indepIndices    = NULL;
+	int         tangentDirNum = 1;  // <----------- set this via config
 
 	
@@ -35,4 +37,8 @@
 
 			if(!(num_dependents*num_independents)) return;
+
+			if (tangentDirNum<1 || tangentDirNum>num_independents) {
+			  _error_("tangentDirNum should be in [1,num_dependents]"); // <------------ fix this error message to relate to config
+			}
 
 			/*retrieve state variable: */
@@ -58,8 +64,9 @@
 			anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
 			anEDF_for_solverx_p->dp_X=xNew<double>(anEDF_for_solverx_p->max_n);
-			anEDF_for_solverx_p->dpp_X=xNew<double>(anEDF_for_solverx_p->max_n, num_independents);
+			anEDF_for_solverx_p->dpp_X=xNew<double>(anEDF_for_solverx_p->max_n, tangentDirNum);
 			anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
 			anEDF_for_solverx_p->dp_Y=xNew<double>(anEDF_for_solverx_p->max_m);
-			anEDF_for_solverx_p->dpp_Y=xNew<double>(anEDF_for_solverx_p->max_m, num_independents);
+			anEDF_for_solverx_p->dpp_Y=xNew<double>(anEDF_for_solverx_p->max_m, tangentDirNum);
+                        /*allocate the space for the parameters to invoke the reverse methods:*/
 			anEDF_for_solverx_p->dp_U=xNew<double>(anEDF_for_solverx_p->max_m);
 			anEDF_for_solverx_p->dpp_U=xNew<double>(num_dependents,anEDF_for_solverx_p->max_m);
@@ -68,21 +75,62 @@
 
 			/* Call AD driver:*/
-			// single direction:
-			double *tangentDir=xNewZeroInit<double>(num_independents);
-			parameters->FindParam(&anIndepNum,AutodiffFosForwardIndexEnum);
-			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");
+			if (tangentDirNum==1) {
+			  // single direction:
+			  double *tangentDir=xNewZeroInit<double>(num_independents);
+			  parameters->FindParam(&anIndepIndex,AutodiffFosForwardIndexEnum);
+			  tangentDir[anIndepIndex]=1.0;
+			  double *jacTimesTangentDir=xNew<double>(num_dependents);
+			  double *theOutput=xNew<double>(num_dependents);
+			  if (fos_forward(1,num_dependents,num_independents, 0, xp, tangentDir, theOutput, jacTimesTangentDir ))
+			    _error_("fos_forward returned non-zero error code");
+			  results->AddObject(new GenericExternalResult<IssmPDouble*>(results->Size()+1,AutodiffJacobianEnum,jacTimesTangentDir,num_dependents,1,1,0.0));
+                          xDelete(theOutput);
+                          xDelete(jacTimesTangentDir);
+                          xDelete(tangentDir);
+			}
+			else {
+			  // full Jacobian or Jacobian projection:
+			  double **jacTimesSeed=xNew<double>(num_dependents,tangentDirNum);
+			  if (tangentDirNum>1 && tangentDirNum<num_dependents) {
+			    double **seed=xNewZeroInit<double>(num_independents,tangentDirNum);
 
-			results->AddObject(new GenericExternalResult<IssmPDouble*>(results->Size()+1,AutodiffJacobianEnum,theJacVecProduct,num_dependents,1,1,0.0));
+			    // <<<<<<< from here <<<<<<<<<<<<
+			    unsigned int * indepIndices=xNew<unsigned int>(tangentDirNum);
+			    for(int i =0; i< tangentDirNum;++i) {
+			      indepIndices[i]=i;
+			    }
+			    // <<<<<<<< to here << get this vector of independent indices from the config - should be 0 based and ideally also unsigned
 
-#if 0
-			double **J=xNew<double>(num_dependents,num_independents);
-			jacobian(1,num_dependents,num_independents,xp,J);
-			results->AddObject(new GenericExternalResult<IssmPDouble*>(results->Size()+1,AutodiffJacobianEnum,*J,num_independents,num_dependents,1,1,0.0));
-			xDelete(J);
-#endif
+			    // collect indices in a set to prevent accidental duplicates as long as we don't do compression:
+			    std::set<unsigned int> anIndexSet;
+			    for (int i=0; i<tangentDirNum; ++i) {
+			      // make sure the index is in range
+			      if (indepIndices[i]>num_independents) {
+			        _error_("indepIndices values must be in [0,num_independents-1]");
+			      }
+			      if (anIndexSet.find(indepIndices[i])!=anIndexSet.end()) {
+			        _error_("duplicate indepIndices values are not allowed until we implement Jacobian decompression");
+			      }
+			      anIndexSet.insert(indepIndices[i]);
+	                      // now populate the seed matrix from the set of independent indices;
+			      // simple setup with a single 1.0 per column and at most a single 1.0 per row
+			      seed[indepIndices[i]][i]=1.0;
+			    }
+			    double *theOutput=xNew<double>(num_dependents);
+			    if (fov_forward(1,num_dependents,num_independents, tangentDirNum, xp, seed, theOutput, jacTimesSeed )) {
+			      _error_("fov_forward returned non-zero error code");
+			    }
+			    xDelete(theOutput);
+			    xDelete(indepIndices);
+			    xDelete(seed);
+			  }
+			  else {
+			    if (jacobian(1,num_dependents,num_independents,xp,jacTimesSeed)) {
+	                        _error_("jacobian returned non-zero error code");
+			    }
+			  }
+			  results->AddObject(new GenericExternalResult<IssmPDouble*>(results->Size()+1,AutodiffJacobianEnum,*jacTimesSeed,num_dependents*tangentDirNum,1,1,0.0));
+			  xDelete(jacTimesSeed);
+			}
 
 			/* delete the allocated space for the parameters:*/