EpetraExt_MultiPointModelEvaluator.cpp

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

Generated on Wed May 12 21:40:37 2010 for EpetraExt by  doxygen 1.4.7