Changeset 13513

Timestamp:

10/03/12 07:41:55 (12 years ago)

Author:

utke

Message:

NEW fov_reverse driver

File:

• issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp (modified) (1 diff)

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

@@ -199 +199 @@
                                 xDelete(aWeightVector);
                         }
-                        else if (strcmp(driver,"fov_reverse")==0) {
+                        else if ((strcmp(driver,"fov_reverse")==0) || (strcmp(driver,"fov_reverse_all")==0)){
+
+                                int* depIndices=NULL;
+                                int weightNum;
+                                int dummy;
+                                double **weightsTimesJac=NULL;
+                                double **weights=NULL;
+                                std::set<unsigned int> anIndexSet;
+
+
+                                /*retrieve directions:*/
+                                if (strcmp(driver,"fov_reverse_all")==0){
+                                        weightNum=num_dependents;
+                                        depIndices=xNewZeroInit<int>(weightNum);
+                                        for(i=0;i<num_dependents;i++)depIndices[i]=1;
+                                }
+                                else{
+                                        parameters->FindParam(&depIndices,&weightNum,&dummy,AutodiffFovForwardIndicesEnum);
+                                }
+
+                                /*Some checks: */
+                                if (weightNum<1 || weightNum>num_dependents) _error_("tangentDirNum should be in [1,num_dependents]");
+
+                                /* full Jacobian or Jacobian projection:*/
+                                weightsTimesJac=xNew<double>(weightNum,num_independents);
+
+                                /*set the forward method function pointers: */
+                                anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
+
                                 /*allocate the space for the parameters to invoke the EDF fos_reverse :*/
-                                anEDF_for_solverx_p->dpp_U=xNew<double>(num_dependents,anEDF_for_solverx_p->max_m);
-                                anEDF_for_solverx_p->dpp_Z=xNew<double>(num_dependents,anEDF_for_solverx_p->max_n);
-                                _error_("driver: " << driver << " not yet supported!");
+                                anEDF_for_solverx_p->dpp_U=xNew<double>(weightNum,anEDF_for_solverx_p->max_m);
+                                anEDF_for_solverx_p->dpp_Z=xNew<double>(weightNum,anEDF_for_solverx_p->max_n);
+
+                                /*seed matrix: */
+                                weights=xNewZeroInit<double>(weightNum,num_dependents);
+
+                                /*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/
+                                for (int i=0; i<weightNum; ++i) {
+                                        /* make sure the index is in range*/
+                                        if (depIndices[i]>num_dependents) {
+                                                _error_("depIndices values must be in [0,num_dependents-1]");
+                                        }
+                                        if (anIndexSet.find(depIndices[i])!=anIndexSet.end()) {
+                                                _error_("duplicate depIndices values are not allowed until we implement Jacobian decompression");
+                                        }
+                                        anIndexSet.insert(depIndices[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*/
+                                        weights[depIndices[i]][i]=1.0;
+                                }
+
+                                /*call driver: */
+                                fov_reverse(1,num_dependents,num_independents, weightNum, weights, weightsTimesJac );
+
+                                /*add to results: */
+                                results->AddObject(new GenericExternalResult<IssmPDouble*>(results->Size()+1,AutodiffJacobianEnum,*weightsTimesJac,weightNum*num_independents,1,1,0.0));
+
+                                /*Free resources: */
+                                xDelete(weights);
+                                xDelete(weightsTimesJac);
+                                xDelete(depIndices);
                         }
                         else _error_("driver: " << driver << " not yet supported!");