EpetraExt_MultiPointModelEvaluator.cpp

Go to the documentation of this file.
00001 #include "EpetraExt_MultiPointModelEvaluator.h"
00002 #include "Epetra_Map.h"
00003 #include "Teuchos_as.hpp"
00004 
00005 EpetraExt::MultiPointModelEvaluator::MultiPointModelEvaluator(
00006     Teuchos::RefCountPtr<EpetraExt::ModelEvaluator> underlyingME_,
00007     const Teuchos::RefCountPtr<EpetraExt::MultiComm> &globalComm_,
00008     const std::vector<Epetra_Vector*> initGuessVec_,
00009     Teuchos::RefCountPtr<std::vector< Teuchos::RefCountPtr<Epetra_Vector> > >  q_vec_,
00010     Teuchos::RefCountPtr<std::vector< Teuchos::RefCountPtr<Epetra_Vector> > >  matching_vec_
00011     ) : 
00012     underlyingME(underlyingME_),
00013     globalComm(globalComm_),
00014     timeStepsOnTimeDomain(globalComm_->NumTimeStepsOnDomain()),
00015     numTimeDomains(globalComm_->NumSubDomains()),
00016     timeDomain(globalComm_->SubDomainRank()),
00017     rowStencil(0),
00018     rowIndex(0),
00019     underlyingNg(0),
00020     q_vec(q_vec_),
00021     matching_vec(matching_vec_)
00022 {
00023   using Teuchos::as;
00024   if (globalComm->MyPID()==0) {
00025      cout  << "----------MultiPoint Partition Info------------"
00026            << "\n\tNumProcs              = " << globalComm->NumProc()
00027            << "\n\tSpatial Decomposition = " << globalComm->SubDomainComm().NumProc()
00028            << "\n\tNumber of Domains     = " << numTimeDomains
00029            << "\n\tSteps on Domain 0     = " << timeStepsOnTimeDomain
00030            << "\n\tTotal Number of Steps = " << globalComm->NumTimeSteps();
00031     cout   << "\n-----------------------------------------------" << endl;
00032     }
00033 
00034    // Construct global block matrix graph from split W and stencil,
00035    // which is just diagonal in this case
00036 
00037    split_W = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(underlyingME->create_W());
00038 
00039    rowStencil = new std::vector< std::vector<int> >(timeStepsOnTimeDomain);
00040    rowIndex = new std::vector<int>;
00041    for (int i=0; i < timeStepsOnTimeDomain; i++) {
00042      (*rowStencil)[i].push_back(0);
00043      (*rowIndex).push_back(i + globalComm->FirstTimeStepOnDomain());
00044    }
00045 
00046    block_W = Teuchos::rcp(new EpetraExt::BlockCrsMatrix(*split_W,
00047                                *rowStencil, *rowIndex, *globalComm));
00048 
00049    // Test for g vector
00050    EpetraExt::ModelEvaluator::OutArgs underlyingOutArgs = underlyingME->createOutArgs();
00051 
00052    underlyingNg = underlyingOutArgs.Ng();
00053    if (underlyingNg) {
00054      if (underlyingOutArgs.supports(OUT_ARG_DgDp,0,0).supports(DERIV_TRANS_MV_BY_ROW))
00055        orientation_DgDp = DERIV_TRANS_MV_BY_ROW;
00056      else
00057        orientation_DgDp = DERIV_MV_BY_COL;
00058    }
00059 
00060    // This code assumes 2 parameter vectors, 1 for opt, second for MultiPoint states
00061    TEST_FOR_EXCEPT(underlyingOutArgs.Np()!=2);
00062 
00063    // temporary quantities
00064    const Epetra_Map& split_map = split_W->RowMatrixRowMap();
00065    num_p0 =  underlyingME_->get_p_map(0)->NumMyElements();
00066    if (underlyingNg)  num_g0 = underlyingME_->get_g_map(0)->NumMyElements();
00067    else num_g0 = 0;
00068    num_dg0dp0 = num_g0 * num_p0;
00069 
00070    // Construct global solution vector, residual vector -- local storage
00071    block_x = new EpetraExt::BlockVector(split_map, block_W->RowMap());
00072    block_f = new EpetraExt::BlockVector(*block_x); 
00073    block_DfDp = new EpetraExt::BlockMultiVector(split_map, block_W->RowMap(), num_p0);
00074     if (underlyingNg)  
00075    block_DgDx = new EpetraExt::BlockMultiVector(split_map, block_W->RowMap(), num_g0);
00076 
00077    // Allocate local storage of epetra vectors
00078    split_x = Teuchos::rcp(new Epetra_Vector(split_map));
00079    split_f = Teuchos::rcp(new Epetra_Vector(split_map));
00080    split_DfDp = Teuchos::rcp(new Epetra_MultiVector(split_map, num_p0));
00081    if (underlyingNg)  
00082      split_DgDx = Teuchos::rcp(new Epetra_MultiVector(split_map, num_g0));
00083    if (underlyingNg) { 
00084      if(orientation_DgDp == DERIV_TRANS_MV_BY_ROW)
00085        split_DgDp = Teuchos::rcp(new Epetra_MultiVector(*(underlyingME_->get_p_map(0)), num_g0));
00086      else
00087        split_DgDp = Teuchos::rcp(new Epetra_MultiVector(*(underlyingME_->get_g_map(0)), num_p0));
00088    } 
00089    if (underlyingNg)  
00090      split_g = Teuchos::rcp(new Epetra_Vector(*(underlyingME_->get_g_map(0))));
00091 
00092    // Packaging required for getting multivectors back as Derivatives
00093    derivMV_DfDp = new EpetraExt::ModelEvaluator::DerivativeMultiVector(split_DfDp);
00094    deriv_DfDp = new EpetraExt::ModelEvaluator::Derivative(*derivMV_DfDp);
00095    if (underlyingNg)  {
00096      derivMV_DgDx = new EpetraExt::ModelEvaluator::DerivativeMultiVector(split_DgDx, DERIV_TRANS_MV_BY_ROW);
00097      deriv_DgDx = new EpetraExt::ModelEvaluator::Derivative(*derivMV_DgDx);
00098      derivMV_DgDp = new EpetraExt::ModelEvaluator::DerivativeMultiVector(split_DgDp, orientation_DgDp);
00099      deriv_DgDp = new EpetraExt::ModelEvaluator::Derivative(*derivMV_DgDp);
00100    }
00101 
00102    // For 4D, we will need the overlap vector and importer between them
00103    // Overlap not needed for MultiPoint -- no overlap between blocks
00104    /*   solutionOverlap = new EpetraExt::BlockVector(split_W->RowMatrixRowMap(),
00105                                                      block_W->ColMap());
00106         overlapImporter = new Epetra_Import(solutionOverlap->Map(), block_x->Map());
00107    */
00108 
00109    // Load initial guess into block solution vector
00110    solution_init = Teuchos::rcp(new EpetraExt::BlockVector(*block_x));
00111    for (int i=0; i < timeStepsOnTimeDomain; i++)
00112            solution_init->LoadBlockValues(*(initGuessVec_[i]), (*rowIndex)[i]);
00113 
00114  
00115    //Prepare logic for matching problem
00116    if (Teuchos::is_null(matching_vec))  matchingProblem = false;
00117    else matchingProblem = true;
00118 
00119    if (matchingProblem) {
00120      TEST_FOR_EXCEPT(as<int>(matching_vec->size())!=timeStepsOnTimeDomain);
00121      TEST_FOR_EXCEPT(!(*matching_vec)[0]->Map().SameAs(*(underlyingME_->get_g_map(0))));
00122      TEST_FOR_EXCEPT(num_g0 != 1); //This restriction may be lifted later
00123    }
00124 }
00125 
00126 EpetraExt::MultiPointModelEvaluator::~MultiPointModelEvaluator()
00127 {
00128   delete block_x;
00129   delete block_f;
00130   delete block_DfDp;
00131   if (underlyingNg)  delete block_DgDx;
00132   delete rowStencil;
00133   delete rowIndex;
00134 
00135   delete derivMV_DfDp;
00136   delete deriv_DfDp;
00137   if (underlyingNg) {
00138     delete derivMV_DgDx;
00139     delete deriv_DgDx;
00140     delete derivMV_DgDp;
00141     delete deriv_DgDp;
00142   }
00143 }
00144 
00145 Teuchos::RefCountPtr<const Epetra_Map> EpetraExt::MultiPointModelEvaluator::get_x_map() const
00146 {
00147   return Teuchos::rcp(&(block_W->OperatorDomainMap()), false);
00148 }
00149 
00150 Teuchos::RefCountPtr<const Epetra_Map> EpetraExt::MultiPointModelEvaluator::get_f_map() const
00151 {
00152   return get_x_map();
00153 }
00154 
00155 Teuchos::RefCountPtr<const Epetra_Map> EpetraExt::MultiPointModelEvaluator::get_p_map(int l) const
00156 {
00157   return underlyingME->get_p_map(l);
00158 }
00159 
00160 Teuchos::RefCountPtr<const Epetra_Map> EpetraExt::MultiPointModelEvaluator::get_g_map(int j) const
00161 {
00162   return underlyingME->get_g_map(j);
00163 }
00164 
00165 Teuchos::RefCountPtr<const Epetra_Vector> EpetraExt::MultiPointModelEvaluator::get_x_init() const
00166 {
00167   return solution_init;
00168 }
00169 
00170 Teuchos::RefCountPtr<const Epetra_Vector> EpetraExt::MultiPointModelEvaluator::get_p_init(int l) const
00171 {
00172   return underlyingME->get_p_init(l);
00173 }
00174 
00175 Teuchos::RefCountPtr<Epetra_Operator> EpetraExt::MultiPointModelEvaluator::create_W() const
00176 {
00177   return block_W;
00178 }
00179 
00180 EpetraExt::ModelEvaluator::InArgs EpetraExt::MultiPointModelEvaluator::createInArgs() const
00181 {
00182   //return underlyingME->createInArgs();
00183   InArgsSetup inArgs;
00184   inArgs.setModelEvalDescription(this->description());
00185   inArgs.set_Np(1);
00186   inArgs.setSupports(IN_ARG_x,true);
00187   return inArgs;
00188 }
00189 
00190 EpetraExt::ModelEvaluator::OutArgs EpetraExt::MultiPointModelEvaluator::createOutArgs() const
00191 {
00192   //return underlyingME->createOutArgs();
00193   OutArgsSetup outArgs;
00194   outArgs.setModelEvalDescription(this->description());
00195   outArgs.set_Np_Ng(1, underlyingNg);
00196   outArgs.setSupports(OUT_ARG_f,true);
00197   outArgs.setSupports(OUT_ARG_W,true);
00198   outArgs.set_W_properties(
00199     DerivativeProperties(
00200       DERIV_LINEARITY_NONCONST
00201       ,DERIV_RANK_FULL
00202       ,true // supportsAdjoint
00203       )
00204     );
00205   outArgs.setSupports(OUT_ARG_DfDp,0,DERIV_MV_BY_COL);
00206   outArgs.set_DfDp_properties(
00207     0,DerivativeProperties(
00208       DERIV_LINEARITY_CONST
00209       ,DERIV_RANK_DEFICIENT
00210       ,true // supportsAdjoint
00211       )
00212     );
00213 
00214   if (underlyingNg) {
00215     outArgs.setSupports(OUT_ARG_DgDx,0,DERIV_TRANS_MV_BY_ROW);
00216     outArgs.set_DgDx_properties(
00217       0,DerivativeProperties(
00218         DERIV_LINEARITY_NONCONST
00219         ,DERIV_RANK_DEFICIENT
00220         ,true // supportsAdjoint
00221         )
00222       );
00223     outArgs.setSupports(OUT_ARG_DgDp,0,0, orientation_DgDp);
00224     outArgs.set_DgDp_properties(
00225       0,0,DerivativeProperties(
00226         DERIV_LINEARITY_NONCONST
00227         ,DERIV_RANK_DEFICIENT
00228         ,true // supportsAdjoint
00229         )
00230       );
00231   }
00232   return outArgs;
00233 }
00234 
00235 void EpetraExt::MultiPointModelEvaluator::evalModel( const InArgs& inArgs,
00236                                             const OutArgs& outArgs ) const
00237 {
00238 
00239   EpetraExt::ModelEvaluator::InArgs  underlyingInArgs  = underlyingME->createInArgs();
00240   EpetraExt::ModelEvaluator::OutArgs underlyingOutArgs = underlyingME->createOutArgs();
00241 
00242   //temp code for multipoint param q vec
00243 /*
00244   Teuchos::RefCountPtr<Epetra_Vector> q =
00245     Teuchos::rcp(new Epetra_Vector(*(underlyingME->get_p_map(1))));
00246 */
00247 
00248   // Parse InArgs
00249   Teuchos::RefCountPtr<const Epetra_Vector> p_in = inArgs.get_p(0);
00250   if (p_in.get()) underlyingInArgs.set_p(0, p_in);
00251 
00252   Teuchos::RefCountPtr<const Epetra_Vector> x_in = inArgs.get_x();
00253   block_x->Epetra_Vector::operator=(*x_in); //copy into block vector
00254 
00255   // Parse OutArgs
00256   Teuchos::RefCountPtr<Epetra_Vector> f_out = outArgs.get_f();
00257 
00258   Teuchos::RefCountPtr<Epetra_Operator> W_out = outArgs.get_W();
00259   Teuchos::RefCountPtr<EpetraExt::BlockCrsMatrix> W_block =
00260      Teuchos::rcp_dynamic_cast<EpetraExt::BlockCrsMatrix>(W_out);
00261 
00262   Teuchos::RefCountPtr<Epetra_Vector> g_out;
00263   if (underlyingNg) g_out = outArgs.get_g(0); 
00264   if (g_out.get()) g_out->PutScalar(0.0);
00265 
00266   EpetraExt::ModelEvaluator::Derivative DfDp_out = outArgs.get_DfDp(0);
00267 
00268   EpetraExt::ModelEvaluator::Derivative DgDx_out;
00269   EpetraExt::ModelEvaluator::Derivative DgDp_out;
00270   if (underlyingNg) {
00271     DgDx_out = outArgs.get_DgDx(0);
00272     DgDp_out = outArgs.get_DgDp(0,0);
00273     if (!DgDx_out.isEmpty()) DgDx_out.getMultiVector()->PutScalar(0.0);
00274     if (!DgDp_out.isEmpty()) DgDp_out.getMultiVector()->PutScalar(0.0);
00275   }
00276 
00277   // For mathcingProblems, g is needed to calc DgDx DgDp, so ask for
00278   //  g even if it isn't requested.
00279   bool need_g = g_out.get();
00280   if (matchingProblem)
00281     if ( !DgDx_out.isEmpty() || !DgDp_out.isEmpty() ) need_g = true;
00282 
00283 
00284   // Begin loop over Points (steps) owned on this proc
00285   for (int i=0; i < timeStepsOnTimeDomain; i++) {
00286 
00287     // Set MultiPoint parameter vector
00288     underlyingInArgs.set_p(1, (*q_vec)[i]);
00289 
00290     // Set InArgs
00291     block_x->ExtractBlockValues(*split_x, (*rowIndex)[i]);
00292     underlyingInArgs.set_x(split_x);
00293 
00294     // Set OutArgs
00295     if (f_out.get()) underlyingOutArgs.set_f(split_f);
00296 
00297     if (need_g) underlyingOutArgs.set_g(0, split_g);
00298 
00299     if (W_out.get()) underlyingOutArgs.set_W(split_W);
00300 
00301     if (!DfDp_out.isEmpty()) underlyingOutArgs.set_DfDp(0, *deriv_DfDp);
00302 
00303     if (!DgDx_out.isEmpty()) underlyingOutArgs.set_DgDx(0, *deriv_DgDx);
00304   
00305     if (!DgDp_out.isEmpty()) underlyingOutArgs.set_DgDp(0, 0, *deriv_DgDp);
00306 
00307     //********Eval Model ********/
00308     underlyingME->evalModel(underlyingInArgs, underlyingOutArgs);
00309     //********Eval Model ********/
00310 
00311     // If matchingProblem, modify all g-related quantitites G = 0.5*(g-g*)^2 / g*^2
00312     if (matchingProblem) {
00313       if (need_g) {
00314         double diff = (*split_g)[0] -  (*(*matching_vec)[i])[0];
00315         double nrmlz = fabs((*(*matching_vec)[i])[0]) + 1.0e-6;
00316         (*split_g)[0] = 0.5 * diff * diff/(nrmlz*nrmlz);
00317         if (!DgDx_out.isEmpty()) split_DgDx->Scale(diff/(nrmlz*nrmlz));
00318         if (!DgDp_out.isEmpty()) split_DgDp->Scale(diff/(nrmlz*nrmlz));
00319       }
00320     }
00321 
00322     // Repackage block components into global block matrx/vector/multivector
00323     if (f_out.get()) block_f->LoadBlockValues(*split_f, (*rowIndex)[i]);
00324     if (W_out.get()) W_block->LoadBlock(*split_W, i, 0);
00325         // note: split_DfDp points inside deriv_DfDp
00326     if (!DfDp_out.isEmpty()) block_DfDp->LoadBlockValues(*split_DfDp, (*rowIndex)[i]);
00327     if (!DgDx_out.isEmpty()) block_DgDx->LoadBlockValues(*split_DgDx, (*rowIndex)[i]);
00328 
00329     // Assemble multiple steps on this domain into g and dgdp(0) vectors
00330     if (g_out.get()) g_out->Update(1.0, *split_g, 1.0);
00331 
00332     if (!DgDp_out.isEmpty())
00333       DgDp_out.getMultiVector()->Update(1.0, *split_DgDp, 1.0);
00334 
00335   } // End loop over multiPoint steps on this domain/cluster
00336 
00337   //Copy block vectors into *_out vectors of same size
00338   if (f_out.get()) f_out->operator=(*block_f);
00339   if (!DfDp_out.isEmpty()) 
00340     DfDp_out.getMultiVector()->operator=(*block_DfDp);
00341   if (!DgDx_out.isEmpty()) 
00342     DgDx_out.getMultiVector()->operator=(*block_DgDx);
00343 
00344   //Sum together obj fn contributions from differnt Domains (clusters).
00345   if (numTimeDomains > 1) {
00346     double factorToZeroOutCopies = 0.0;
00347     if (globalComm->SubDomainComm().MyPID()==0) factorToZeroOutCopies = 1.0;
00348     if (g_out.get()) {
00349       (*g_out).Scale(factorToZeroOutCopies);
00350       double* vPtr = &((*g_out)[0]);
00351       Epetra_Vector tmp = *(g_out.get());
00352       globalComm->SumAll( &(tmp[0]), vPtr, num_g0);
00353     }
00354     if (!DgDp_out.isEmpty()) {
00355       DgDp_out.getMultiVector()->Scale(factorToZeroOutCopies);
00356       double* mvPtr = (*DgDp_out.getMultiVector())[0];
00357       Epetra_MultiVector tmp = *(DgDp_out.getMultiVector());
00358       globalComm->SumAll(tmp[0], mvPtr, num_dg0dp0);
00359     }
00360   }
00361 }

Generated on Wed May 12 21:24:46 2010 for EpetraExt by  doxygen 1.4.7