Changeset 13452

Timestamp:

09/26/12 15:32:42 (12 years ago)

Author:

utke

Message:

CHG allow more forward variants to be used (note that jacobian may call reverse though)

File:

• issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp

@@ -2 +2 @@
  * \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 +21 @@
         IssmDouble *axp                = NULL;
         double     *xp                 = NULL;
-        int         anIndepNum;
+        int         anIndepIndex;
+        unsigned int * indepIndices    = NULL;
+        int         tangentDirNum = 1;  // <----------- set this via config
 
         
@@ -35 +37 @@
 
                         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 +64 @@
                         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 +75 @@
 
                         /* 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:*/