EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_MultiPointModelEvaluator.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ***********************************************************************
00004 //
00005 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00006 //                 Copyright (2011) Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include "EpetraExt_MultiPointModelEvaluator.h"
00045 #include "Epetra_Map.h"
00046 #include "Teuchos_as.hpp"
00047 
00048 EpetraExt::MultiPointModelEvaluator::MultiPointModelEvaluator(
00049     Teuchos::RefCountPtr<EpetraExt::ModelEvaluator> underlyingME_,
00050     const Teuchos::RefCountPtr<EpetraExt::MultiComm> &globalComm_,
00051     const std::vector<Epetra_Vector*> initGuessVec_,
00052     Teuchos::RefCountPtr<std::vector< Teuchos::RefCountPtr<Epetra_Vector> > >  q_vec_,
00053     Teuchos::RefCountPtr<std::vector< Teuchos::RefCountPtr<Epetra_Vector> > >  matching_vec_
00054     ) : 
00055     underlyingME(underlyingME_),
00056     globalComm(globalComm_),
00057     underlyingNg(0),
00058     timeStepsOnTimeDomain(globalComm_->NumTimeStepsOnDomain()),
00059     numTimeDomains(globalComm_->NumSubDomains()),
00060     timeDomain(globalComm_->SubDomainRank()),
00061     rowStencil(0),
00062     rowIndex(0),
00063     q_vec(q_vec_),
00064     matching_vec(matching_vec_)
00065 {
00066   using Teuchos::as;
00067   if (globalComm->MyPID()==0) {
00068      cout  << "----------MultiPoint Partition Info------------"
00069            << "\n\tNumProcs              = " << globalComm->NumProc()
00070            << "\n\tSpatial Decomposition = " << globalComm->SubDomainComm().NumProc()
00071            << "\n\tNumber of Domains     = " << numTimeDomains
00072            << "\n\tSteps on Domain 0     = " << timeStepsOnTimeDomain
00073            << "\n\tTotal Number of Steps = " << globalComm->NumTimeSteps();
00074     cout   << "\n-----------------------------------------------" << endl;
00075     }
00076 
00077    // Construct global block matrix graph from split W and stencil,
00078    // which is just diagonal in this case
00079 
00080    split_W = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(underlyingME->create_W());
00081 
00082    rowStencil = new std::vector< std::vector<int> >(timeStepsOnTimeDomain);
00083    rowIndex = new std::vector<int>;
00084    for (int i=0; i < timeStepsOnTimeDomain; i++) {
00085      (*rowStencil)[i].push_back(0);
00086      (*rowIndex).push_back(i + globalComm->FirstTimeStepOnDomain());
00087    }
00088 
00089    block_W = Teuchos::rcp(new EpetraExt::BlockCrsMatrix(*split_W,
00090                                *rowStencil, *rowIndex, *globalComm));
00091 
00092    // Test for g vector
00093    EpetraExt::ModelEvaluator::OutArgs underlyingOutArgs = underlyingME->createOutArgs();
00094 
00095    underlyingNg = underlyingOutArgs.Ng();
00096    if (underlyingNg) {
00097      if (underlyingOutArgs.supports(OUT_ARG_DgDp,0,0).supports(DERIV_TRANS_MV_BY_ROW))
00098        orientation_DgDp = DERIV_TRANS_MV_BY_ROW;
00099      else
00100        orientation_DgDp = DERIV_MV_BY_COL;
00101    }
00102 
00103    // This code assumes 2 parameter vectors, 1 for opt, second for MultiPoint states
00104    TEUCHOS_TEST_FOR_EXCEPT(underlyingOutArgs.Np()!=2);
00105 
00106    // temporary quantities
00107    const Epetra_Map& split_map = split_W->RowMatrixRowMap();
00108    num_p0 =  underlyingME_->get_p_map(0)->NumMyElements();
00109    if (underlyingNg)  num_g0 = underlyingME_->get_g_map(0)->NumMyElements();
00110    else num_g0 = 0;
00111    num_dg0dp0 = num_g0 * num_p0;
00112 
00113    // Construct global solution vector, residual vector -- local storage
00114    block_x = new EpetraExt::BlockVector(split_map, block_W->RowMap());
00115    block_f = new EpetraExt::BlockVector(*block_x); 
00116    block_DfDp = new EpetraExt::BlockMultiVector(split_map, block_W->RowMap(), num_p0);
00117     if (underlyingNg)  
00118    block_DgDx = new EpetraExt::BlockMultiVector(split_map, block_W->RowMap(), num_g0);
00119 
00120    // Allocate local storage of epetra vectors
00121    split_x = Teuchos::rcp(new Epetra_Vector(split_map));
00122    split_f = Teuchos::rcp(new Epetra_Vector(split_map));
00123    split_DfDp = Teuchos::rcp(new Epetra_MultiVector(split_map, num_p0));
00124    if (underlyingNg)  
00125      split_DgDx = Teuchos::rcp(new Epetra_MultiVector(split_map, num_g0));
00126    if (underlyingNg) { 
00127      if(orientation_DgDp == DERIV_TRANS_MV_BY_ROW)
00128        split_DgDp = Teuchos::rcp(new Epetra_MultiVector(*(underlyingME_->get_p_map(0)), num_g0));
00129      else
00130        split_DgDp = Teuchos::rcp(new Epetra_MultiVector(*(underlyingME_->get_g_map(0)), num_p0));
00131    } 
00132    if (underlyingNg)  
00133      split_g = Teuchos::rcp(new Epetra_Vector(*(underlyingME_->get_g_map(0))));
00134 
00135    // Packaging required for getting multivectors back as Derivatives
00136    derivMV_DfDp = new EpetraExt::ModelEvaluator::DerivativeMultiVector(split_DfDp);
00137    deriv_DfDp = new EpetraExt::ModelEvaluator::Derivative(*derivMV_DfDp);
00138    if (underlyingNg)  {
00139      derivMV_DgDx = new EpetraExt::ModelEvaluator::DerivativeMultiVector(split_DgDx, DERIV_TRANS_MV_BY_ROW);
00140      deriv_DgDx = new EpetraExt::ModelEvaluator::Derivative(*derivMV_DgDx);
00141      derivMV_DgDp = new EpetraExt::ModelEvaluator::DerivativeMultiVector(split_DgDp, orientation_DgDp);
00142      deriv_DgDp = new EpetraExt::ModelEvaluator::Derivative(*derivMV_DgDp);
00143    }
00144 
00145    // For 4D, we will need the overlap vector and importer between them
00146    // Overlap not needed for MultiPoint -- no overlap between blocks
00147    /*   solutionOverlap = new EpetraExt::BlockVector(split_W->RowMatrixRowMap(),
00148                                                      block_W->ColMap());
00149         overlapImporter = new Epetra_Import(solutionOverlap->Map(), block_x->Map());
00150    */
00151 
00152    // Load initial guess into block solution vector
00153    solution_init = Teuchos::rcp(new EpetraExt::BlockVector(*block_x));
00154    for (int i=0; i < timeStepsOnTimeDomain; i++)
00155            solution_init->LoadBlockValues(*(initGuessVec_[i]), (*rowIndex)[i]);
00156 
00157  
00158    //Prepare logic for matching problem
00159    if (Teuchos::is_null(matching_vec))  matchingProblem = false;
00160    else matchingProblem = true;
00161 
00162    if (matchingProblem) {
00163      TEUCHOS_TEST_FOR_EXCEPT(as<int>(matching_vec->size())!=timeStepsOnTimeDomain);
00164      TEUCHOS_TEST_FOR_EXCEPT(!(*matching_vec)[0]->Map().SameAs(*(underlyingME_->get_g_map(0))));
00165      TEUCHOS_TEST_FOR_EXCEPT(num_g0 != 1); //This restriction may be lifted later
00166    }
00167 }
00168 
00169 EpetraExt::MultiPointModelEvaluator::~MultiPointModelEvaluator()
00170 {
00171   delete block_x;
00172   delete block_f;
00173   delete block_DfDp;
00174   if (underlyingNg)  delete block_DgDx;
00175   delete rowStencil;
00176   delete rowIndex;
00177 
00178   delete derivMV_DfDp;
00179   delete deriv_DfDp;
00180   if (underlyingNg) {
00181     delete derivMV_DgDx;
00182     delete deriv_DgDx;
00183     delete derivMV_DgDp;
00184     delete deriv_DgDp;
00185   }
00186 }
00187 
00188 Teuchos::RefCountPtr<const Epetra_Map> EpetraExt::MultiPointModelEvaluator::get_x_map() const
00189 {
00190   return Teuchos::rcp(&(block_W->OperatorDomainMap()), false);
00191 }
00192 
00193 Teuchos::RefCountPtr<const Epetra_Map> EpetraExt::MultiPointModelEvaluator::get_f_map() const
00194 {
00195   return get_x_map();
00196 }
00197 
00198 Teuchos::RefCountPtr<const Epetra_Map> EpetraExt::MultiPointModelEvaluator::get_p_map(int l) const
00199 {
00200   return underlyingME->get_p_map(l);
00201 }
00202 
00203 Teuchos::RefCountPtr<const Epetra_Map> EpetraExt::MultiPointModelEvaluator::get_g_map(int j) const
00204 {
00205   return underlyingME->get_g_map(j);
00206 }
00207 
00208 Teuchos::RefCountPtr<const Epetra_Vector> EpetraExt::MultiPointModelEvaluator::get_x_init() const
00209 {
00210   return solution_init;
00211 }
00212 
00213 Teuchos::RefCountPtr<const Epetra_Vector> EpetraExt::MultiPointModelEvaluator::get_p_init(int l) const
00214 {
00215   return underlyingME->get_p_init(l);
00216 }
00217 
00218 Teuchos::RefCountPtr<Epetra_Operator> EpetraExt::MultiPointModelEvaluator::create_W() const
00219 {
00220   return block_W;
00221 }
00222 
00223 EpetraExt::ModelEvaluator::InArgs EpetraExt::MultiPointModelEvaluator::createInArgs() const
00224 {
00225   //return underlyingME->createInArgs();
00226   InArgsSetup inArgs;
00227   inArgs.setModelEvalDescription(this->description());
00228   inArgs.set_Np(1);
00229   inArgs.setSupports(IN_ARG_x,true);
00230   return inArgs;
00231 }
00232 
00233 EpetraExt::ModelEvaluator::OutArgs EpetraExt::MultiPointModelEvaluator::createOutArgs() const
00234 {
00235   //return underlyingME->createOutArgs();
00236   OutArgsSetup outArgs;
00237   outArgs.setModelEvalDescription(this->description());
00238   outArgs.set_Np_Ng(1, underlyingNg);
00239   outArgs.setSupports(OUT_ARG_f,true);
00240   outArgs.setSupports(OUT_ARG_W,true);
00241   outArgs.set_W_properties(
00242     DerivativeProperties(
00243       DERIV_LINEARITY_NONCONST
00244       ,DERIV_RANK_FULL
00245       ,true // supportsAdjoint
00246       )
00247     );
00248   outArgs.setSupports(OUT_ARG_DfDp,0,DERIV_MV_BY_COL);
00249   outArgs.set_DfDp_properties(
00250     0,DerivativeProperties(
00251       DERIV_LINEARITY_CONST
00252       ,DERIV_RANK_DEFICIENT
00253       ,true // supportsAdjoint
00254       )
00255     );
00256 
00257   if (underlyingNg) {
00258     outArgs.setSupports(OUT_ARG_DgDx,0,DERIV_TRANS_MV_BY_ROW);
00259     outArgs.set_DgDx_properties(
00260       0,DerivativeProperties(
00261         DERIV_LINEARITY_NONCONST
00262         ,DERIV_RANK_DEFICIENT
00263         ,true // supportsAdjoint
00264         )
00265       );
00266     outArgs.setSupports(OUT_ARG_DgDp,0,0, orientation_DgDp);
00267     outArgs.set_DgDp_properties(
00268       0,0,DerivativeProperties(
00269         DERIV_LINEARITY_NONCONST
00270         ,DERIV_RANK_DEFICIENT
00271         ,true // supportsAdjoint
00272         )
00273       );
00274   }
00275   return outArgs;
00276 }
00277 
00278 void EpetraExt::MultiPointModelEvaluator::evalModel( const InArgs& inArgs,
00279                                             const OutArgs& outArgs ) const
00280 {
00281 
00282   EpetraExt::ModelEvaluator::InArgs  underlyingInArgs  = underlyingME->createInArgs();
00283   EpetraExt::ModelEvaluator::OutArgs underlyingOutArgs = underlyingME->createOutArgs();
00284 
00285   //temp code for multipoint param q vec
00286 /*
00287   Teuchos::RefCountPtr<Epetra_Vector> q =
00288     Teuchos::rcp(new Epetra_Vector(*(underlyingME->get_p_map(1))));
00289 */
00290 
00291   // Parse InArgs
00292   Teuchos::RefCountPtr<const Epetra_Vector> p_in = inArgs.get_p(0);
00293   if (p_in.get()) underlyingInArgs.set_p(0, p_in);
00294 
00295   Teuchos::RefCountPtr<const Epetra_Vector> x_in = inArgs.get_x();
00296   block_x->Epetra_Vector::operator=(*x_in); //copy into block vector
00297 
00298   // Parse OutArgs
00299   Teuchos::RefCountPtr<Epetra_Vector> f_out = outArgs.get_f();
00300 
00301   Teuchos::RefCountPtr<Epetra_Operator> W_out = outArgs.get_W();
00302   Teuchos::RefCountPtr<EpetraExt::BlockCrsMatrix> W_block =
00303      Teuchos::rcp_dynamic_cast<EpetraExt::BlockCrsMatrix>(W_out);
00304 
00305   Teuchos::RefCountPtr<Epetra_Vector> g_out;
00306   if (underlyingNg) g_out = outArgs.get_g(0); 
00307   if (g_out.get()) g_out->PutScalar(0.0);
00308 
00309   EpetraExt::ModelEvaluator::Derivative DfDp_out = outArgs.get_DfDp(0);
00310 
00311   EpetraExt::ModelEvaluator::Derivative DgDx_out;
00312   EpetraExt::ModelEvaluator::Derivative DgDp_out;
00313   if (underlyingNg) {
00314     DgDx_out = outArgs.get_DgDx(0);
00315     DgDp_out = outArgs.get_DgDp(0,0);
00316     if (!DgDx_out.isEmpty()) DgDx_out.getMultiVector()->PutScalar(0.0);
00317     if (!DgDp_out.isEmpty()) DgDp_out.getMultiVector()->PutScalar(0.0);
00318   }
00319 
00320   // For mathcingProblems, g is needed to calc DgDx DgDp, so ask for
00321   //  g even if it isn't requested.
00322   bool need_g = g_out.get();
00323   if (matchingProblem)
00324     if ( !DgDx_out.isEmpty() || !DgDp_out.isEmpty() ) need_g = true;
00325 
00326 
00327   // Begin loop over Points (steps) owned on this proc
00328   for (int i=0; i < timeStepsOnTimeDomain; i++) {
00329 
00330     // Set MultiPoint parameter vector
00331     underlyingInArgs.set_p(1, (*q_vec)[i]);
00332 
00333     // Set InArgs
00334     block_x->ExtractBlockValues(*split_x, (*rowIndex)[i]);
00335     underlyingInArgs.set_x(split_x);
00336 
00337     // Set OutArgs
00338     if (f_out.get()) underlyingOutArgs.set_f(split_f);
00339 
00340     if (need_g) underlyingOutArgs.set_g(0, split_g);
00341 
00342     if (W_out.get()) underlyingOutArgs.set_W(split_W);
00343 
00344     if (!DfDp_out.isEmpty()) underlyingOutArgs.set_DfDp(0, *deriv_DfDp);
00345 
00346     if (!DgDx_out.isEmpty()) underlyingOutArgs.set_DgDx(0, *deriv_DgDx);
00347   
00348     if (!DgDp_out.isEmpty()) underlyingOutArgs.set_DgDp(0, 0, *deriv_DgDp);
00349 
00350     //********Eval Model ********/
00351     underlyingME->evalModel(underlyingInArgs, underlyingOutArgs);
00352     //********Eval Model ********/
00353 
00354     // If matchingProblem, modify all g-related quantitites G = 0.5*(g-g*)^2 / g*^2
00355     if (matchingProblem) {
00356       if (need_g) {
00357         double diff = (*split_g)[0] -  (*(*matching_vec)[i])[0];
00358         double nrmlz = fabs((*(*matching_vec)[i])[0]) + 1.0e-6;
00359         (*split_g)[0] = 0.5 * diff * diff/(nrmlz*nrmlz);
00360         if (!DgDx_out.isEmpty()) split_DgDx->Scale(diff/(nrmlz*nrmlz));
00361         if (!DgDp_out.isEmpty()) split_DgDp->Scale(diff/(nrmlz*nrmlz));
00362       }
00363     }
00364 
00365     // Repackage block components into global block matrx/vector/multivector
00366     if (f_out.get()) block_f->LoadBlockValues(*split_f, (*rowIndex)[i]);
00367     if (W_out.get()) W_block->LoadBlock(*split_W, i, 0);
00368         // note: split_DfDp points inside deriv_DfDp
00369     if (!DfDp_out.isEmpty()) block_DfDp->LoadBlockValues(*split_DfDp, (*rowIndex)[i]);
00370     if (!DgDx_out.isEmpty()) block_DgDx->LoadBlockValues(*split_DgDx, (*rowIndex)[i]);
00371 
00372     // Assemble multiple steps on this domain into g and dgdp(0) vectors
00373     if (g_out.get()) g_out->Update(1.0, *split_g, 1.0);
00374 
00375     if (!DgDp_out.isEmpty())
00376       DgDp_out.getMultiVector()->Update(1.0, *split_DgDp, 1.0);
00377 
00378   } // End loop over multiPoint steps on this domain/cluster
00379 
00380   //Copy block vectors into *_out vectors of same size
00381   if (f_out.get()) f_out->operator=(*block_f);
00382   if (!DfDp_out.isEmpty()) 
00383     DfDp_out.getMultiVector()->operator=(*block_DfDp);
00384   if (!DgDx_out.isEmpty()) 
00385     DgDx_out.getMultiVector()->operator=(*block_DgDx);
00386 
00387   //Sum together obj fn contributions from differnt Domains (clusters).
00388   if (numTimeDomains > 1) {
00389     double factorToZeroOutCopies = 0.0;
00390     if (globalComm->SubDomainComm().MyPID()==0) factorToZeroOutCopies = 1.0;
00391     if (g_out.get()) {
00392       (*g_out).Scale(factorToZeroOutCopies);
00393       double* vPtr = &((*g_out)[0]);
00394       Epetra_Vector tmp = *(g_out.get());
00395       globalComm->SumAll( &(tmp[0]), vPtr, num_g0);
00396     }
00397     if (!DgDp_out.isEmpty()) {
00398       DgDp_out.getMultiVector()->Scale(factorToZeroOutCopies);
00399       double* mvPtr = (*DgDp_out.getMultiVector())[0];
00400       Epetra_MultiVector tmp = *(DgDp_out.getMultiVector());
00401       globalComm->SumAll(tmp[0], mvPtr, num_dg0dp0);
00402     }
00403   }
00404 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines