Thyra_EpetraModelEvaluator.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "Thyra_EpetraModelEvaluator.hpp"
00030 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00031 #include "Thyra_EpetraThyraWrappers.hpp"
00032 #include "Thyra_EpetraLinearOp.hpp"
00033 #include "Teuchos_Time.hpp"
00034 
00035 namespace Thyra {
00036 
00037 // Constructors/initializers/accessors.
00038 
00039 EpetraModelEvaluator::EpetraModelEvaluator()
00040   :finalPointWasSolved_(false)
00041 {}
00042 
00043 EpetraModelEvaluator::EpetraModelEvaluator(
00044   const Teuchos::RefCountPtr<const EpetraExt::ModelEvaluator>         &epetraModel
00045   ,const Teuchos::RefCountPtr<LinearOpWithSolveFactoryBase<double> >  &W_factory
00046   )
00047 {
00048   initialize(epetraModel,W_factory);
00049 }
00050 
00051 void EpetraModelEvaluator::initialize(
00052   const Teuchos::RefCountPtr<const EpetraExt::ModelEvaluator>         &epetraModel
00053   ,const Teuchos::RefCountPtr<LinearOpWithSolveFactoryBase<double> >  &W_factory
00054   )
00055 {
00056   typedef ModelEvaluatorBase MEB;
00057   //
00058   epetraModel_ = epetraModel;
00059   //
00060   W_factory_ = W_factory;
00061   //
00062   x_space_ = create_VectorSpace( x_map_ = epetraModel_->get_x_map() );
00063   f_space_ = create_VectorSpace( f_map_ = epetraModel_->get_f_map() );
00064   //
00065   EpetraExt::ModelEvaluator::InArgs inArgs = epetraModel_->createInArgs();
00066   p_map_.resize(inArgs.Np()); p_space_.resize(inArgs.Np());
00067   for( int l = 0; l < static_cast<int>(p_space_.size()); ++l )
00068     p_space_[l] = create_VectorSpace( p_map_[l] = epetraModel_->get_p_map(l) );
00069   //
00070   EpetraExt::ModelEvaluator::OutArgs outArgs = epetraModel_->createOutArgs();
00071   g_map_.resize(outArgs.Ng()); g_space_.resize(outArgs.Ng());
00072   for( int j = 0; j < static_cast<int>(g_space_.size()); ++j )
00073     g_space_[j] = create_VectorSpace( g_map_[j] = epetraModel_->get_g_map(j) );
00074   //
00075   initialGuess_ = this->createInArgs();
00076   lowerBounds_ = this->createInArgs();
00077   upperBounds_ = this->createInArgs();
00078   if(initialGuess_.supports(MEB::IN_ARG_x)) {
00079     initialGuess_.set_x( create_Vector( epetraModel_->get_x_init(), x_space_ ) );
00080     lowerBounds_.set_x( create_Vector( epetraModel_->get_x_lower_bounds(), x_space_ ) );
00081     upperBounds_.set_x( create_Vector( epetraModel_->get_x_upper_bounds(), x_space_ ) );
00082   }
00083   if(initialGuess_.supports(MEB::IN_ARG_x_dot)) {
00084     initialGuess_.set_x_dot( create_Vector( epetraModel_->get_x_dot_init(), x_space_ ) );
00085   }
00086   for( int l = 0; l < static_cast<int>(p_space_.size()); ++l ) {
00087     initialGuess_.set_p( l, create_Vector( epetraModel_->get_p_init(l), p_space_[l] ) );
00088     lowerBounds_.set_p( l, create_Vector( epetraModel_->get_p_lower_bounds(l), p_space_[l] ) );
00089     upperBounds_.set_p( l, create_Vector( epetraModel_->get_p_upper_bounds(l), p_space_[l] ) );
00090   }
00091   if(initialGuess_.supports(MEB::IN_ARG_t)) {
00092     initialGuess_.set_t(epetraModel_->get_t_init());
00093     lowerBounds_.set_t(epetraModel_->get_t_lower_bound());
00094     upperBounds_.set_t(epetraModel_->get_t_upper_bound());
00095   }
00096   finalPointWasSolved_ = false;
00097 }
00098 
00099 Teuchos::RefCountPtr<const EpetraExt::ModelEvaluator> EpetraModelEvaluator::getEpetraModel() const
00100 {
00101   return epetraModel_;
00102 }
00103 
00104 void EpetraModelEvaluator::setInitialGuess( const ModelEvaluatorBase::InArgs<double>& initialGuess )
00105 {
00106   initialGuess_.setArgs(initialGuess);
00107 }
00108 
00109 void EpetraModelEvaluator::uninitialize(
00110   Teuchos::RefCountPtr<const EpetraExt::ModelEvaluator>         *epetraModel
00111   ,Teuchos::RefCountPtr<LinearOpWithSolveFactoryBase<double> >  *W_factory
00112   )
00113 {
00114   if(epetraModel) *epetraModel = epetraModel_;
00115   if(W_factory) *W_factory = W_factory_;
00116   epetraModel_ = Teuchos::null;
00117   W_factory_ = Teuchos::null;
00118 }
00119 
00120 const ModelEvaluatorBase::InArgs<double>&
00121 EpetraModelEvaluator::getFinalPoint() const
00122 {
00123   return finalPoint_;
00124 }
00125 
00126 bool EpetraModelEvaluator::finalPointWasSolved() const
00127 {
00128   return finalPointWasSolved_;
00129 }
00130 
00131 // Overridden from ModelEvaulator.
00132 
00133 int EpetraModelEvaluator::Np() const
00134 {
00135   return p_space_.size();
00136 }
00137 
00138 int EpetraModelEvaluator::Ng() const
00139 {
00140   return g_space_.size();
00141 }
00142 
00143 Teuchos::RefCountPtr<const VectorSpaceBase<double> >
00144 EpetraModelEvaluator::get_x_space() const
00145 {
00146   return x_space_;
00147 }
00148 
00149 Teuchos::RefCountPtr<const VectorSpaceBase<double> >
00150 EpetraModelEvaluator::get_f_space() const
00151 {
00152   return f_space_;
00153 }
00154 
00155 Teuchos::RefCountPtr<const VectorSpaceBase<double> >
00156 EpetraModelEvaluator::get_p_space(int l) const
00157 {
00158   TEST_FOR_EXCEPT( ! ( 0 <= l && l < this->Np() ) );
00159   return p_space_[l];
00160 }
00161 
00162 Teuchos::RefCountPtr<const VectorSpaceBase<double> >
00163 EpetraModelEvaluator::get_g_space(int j) const
00164 {
00165   TEST_FOR_EXCEPT( ! ( 0 <= j && j < this->Ng() ) );
00166   return g_space_[j];
00167 }
00168 
00169 ModelEvaluatorBase::InArgs<double> EpetraModelEvaluator::getNominalValues() const
00170 {
00171   return initialGuess_;
00172 }
00173 
00174 ModelEvaluatorBase::InArgs<double> EpetraModelEvaluator::getLowerBounds() const
00175 {
00176   return lowerBounds_;
00177 }
00178 
00179 ModelEvaluatorBase::InArgs<double> EpetraModelEvaluator::getUpperBounds() const
00180 {
00181   return upperBounds_;
00182 }
00183 
00184 Teuchos::RefCountPtr<LinearOpWithSolveBase<double> >
00185 EpetraModelEvaluator::create_W() const
00186 {
00187   TEST_FOR_EXCEPTION(
00188     W_factory_.get()==NULL, std::logic_error
00189     ,"Thyra::EpetraModelEvaluator::create_W(): Error, the client did not set a LinearOpWithSolveFactoryBase"
00190     " object for W!"
00191     );
00192   W_factory_->setOStream(this->getOStream());
00193   W_factory_->setVerbLevel(this->getVerbLevel());
00194   return W_factory_->createOp();
00195 }
00196 
00197 Teuchos::RefCountPtr<LinearOpBase<double> >
00198 EpetraModelEvaluator::create_W_op() const
00199 {
00200   return Teuchos::rcp(new Thyra::EpetraLinearOp());
00201 }
00202 
00203 Teuchos::RefCountPtr<LinearOpBase<double> >
00204 EpetraModelEvaluator::create_DfDp_op(int l) const
00205 {
00206   TEST_FOR_EXCEPT(true);
00207   return Teuchos::null;
00208 }
00209 
00210 Teuchos::RefCountPtr<LinearOpBase<double> >
00211 EpetraModelEvaluator::create_DgDx_op(int j) const
00212 {
00213   TEST_FOR_EXCEPT(true);
00214   return Teuchos::null;
00215 }
00216 
00217 Teuchos::RefCountPtr<LinearOpBase<double> >
00218 EpetraModelEvaluator::create_DgDp_op( int j, int l ) const
00219 {
00220   TEST_FOR_EXCEPT(true);
00221   return Teuchos::null;
00222 }
00223 
00224 ModelEvaluatorBase::InArgs<double> EpetraModelEvaluator::createInArgs() const
00225 {
00226   const EpetraExt::ModelEvaluator &epetraModel = *epetraModel_;
00227   InArgsSetup<double> inArgs;
00228   typedef EpetraExt::ModelEvaluator EME;
00229   EME::InArgs epetraInArgs = epetraModel.createInArgs();
00230   inArgs.setModelEvalDescription(this->description());
00231   inArgs.set_Np(epetraInArgs.Np());
00232   inArgs.setSupports(IN_ARG_x_dot,epetraInArgs.supports(EME::IN_ARG_x_dot));
00233   inArgs.setSupports(IN_ARG_x,epetraInArgs.supports(EME::IN_ARG_x));
00234   inArgs.setSupports(IN_ARG_x_dot_poly,
00235          epetraInArgs.supports(EME::IN_ARG_x_dot_poly));
00236   inArgs.setSupports(IN_ARG_x_poly,epetraInArgs.supports(EME::IN_ARG_x_poly));
00237   inArgs.setSupports(IN_ARG_t,epetraInArgs.supports(EME::IN_ARG_t));
00238   inArgs.setSupports(IN_ARG_alpha,epetraInArgs.supports(EME::IN_ARG_alpha));
00239   inArgs.setSupports(IN_ARG_beta,epetraInArgs.supports(EME::IN_ARG_beta));
00240   return inArgs;
00241 }
00242 
00243 ModelEvaluatorBase::OutArgs<double> EpetraModelEvaluator::createOutArgs() const
00244 {
00245   const EpetraExt::ModelEvaluator &epetraModel = *epetraModel_;
00246   OutArgsSetup<double> outArgs;
00247   typedef EpetraExt::ModelEvaluator EME;
00248   EME::InArgs  epetraInArgs  = epetraModel.createInArgs();
00249   EME::OutArgs epetraOutArgs = epetraModel.createOutArgs();
00250   const int Np = epetraOutArgs.Np(), Ng = epetraOutArgs.Ng();
00251   outArgs.setModelEvalDescription(this->description());
00252   outArgs.set_Np_Ng(Np,Ng);
00253   outArgs.setSupports(OUT_ARG_f,epetraOutArgs.supports(EME::OUT_ARG_f));
00254   outArgs.setSupports(OUT_ARG_W,epetraOutArgs.supports(EME::OUT_ARG_W)&&W_factory_.get()!=NULL);
00255   outArgs.setSupports(OUT_ARG_W_op,epetraOutArgs.supports(EME::OUT_ARG_W));
00256   outArgs.set_W_properties(convert(epetraOutArgs.get_W_properties()));
00257   for(int l=0; l<Np; ++l) {
00258     outArgs.setSupports(OUT_ARG_DfDp,l,convert(epetraOutArgs.supports(EME::OUT_ARG_DfDp,l)));
00259     if(!outArgs.supports(OUT_ARG_DfDp,l).none())
00260       outArgs.set_DfDp_properties(l,convert(epetraOutArgs.get_DfDp_properties(l)));
00261   }
00262   for(int j=0; j<Ng; ++j) {
00263     outArgs.setSupports(OUT_ARG_DgDx,j,convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx,j)));
00264     if(!outArgs.supports(OUT_ARG_DgDx,j).none())
00265       outArgs.set_DgDx_properties(j,convert(epetraOutArgs.get_DgDx_properties(j)));
00266   }
00267   for(int j=0; j<Ng; ++j) for(int l=0; l<Np; ++l) {
00268     outArgs.setSupports(OUT_ARG_DgDp,j,l,convert(epetraOutArgs.supports(EME::OUT_ARG_DgDp,j,l)));
00269     if(!outArgs.supports(OUT_ARG_DgDp,j,l).none())
00270       outArgs.set_DgDp_properties(j,l,convert(epetraOutArgs.get_DgDp_properties(j,l)));
00271   }
00272   outArgs.setSupports(OUT_ARG_f_poly,epetraOutArgs.supports(EME::OUT_ARG_f_poly));
00273   return outArgs;
00274 }
00275 
00276 void EpetraModelEvaluator::evalModel( const InArgs<double>& inArgs_in, const OutArgs<double>& outArgs ) const
00277 {
00278 
00279   using Thyra::get_Epetra_Vector;
00280   using Teuchos::RefCountPtr;
00281   using Teuchos::rcp_const_cast;
00282   using Teuchos::rcp_dynamic_cast;
00283   using Teuchos::OSTab;
00284 
00285   typedef EpetraExt::ModelEvaluator EME;
00286 
00287   InArgs<double> inArgs = initialGuess_; // Make sure we grab the initial guess first!
00288   inArgs.setArgs(inArgs_in);
00289 
00290   Teuchos::Time totalTimer(""), timer("");
00291   totalTimer.start(true);
00292 
00293   const Teuchos::RefCountPtr<Teuchos::FancyOStream> out       = this->getOStream();
00294   const Teuchos::EVerbosityLevel                    verbLevel = this->getVerbLevel();
00295   Teuchos::OSTab tab(out);
00296   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00297     *out << "\nEntering Thyra::EpetraModelEvaluator::evalModel(...) ...\n";
00298 
00299   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_EXTREME))
00300     *out
00301       << "\ninArgs =\n" << Teuchos::describe(inArgs,verbLevel)
00302       << "\noutArgs on input =\n" << Teuchos::describe(outArgs,Teuchos::VERB_LOW);
00303   
00304   typedef Teuchos::VerboseObjectTempState<LinearOpWithSolveFactoryBase<double> > VOTSLOWSF;
00305   VOTSLOWSF W_factory_outputTempState(W_factory_,out,verbLevel);
00306   
00307   // InArgs
00308   
00309   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00310     *out << "\nSetting-up/creating input arguments ...\n";
00311   timer.start(true);
00312 
00313   EME::InArgs epetraInArgs = epetraModel_->createInArgs();
00314 
00315   RefCountPtr<const VectorBase<double> > x_dot;
00316   if( inArgs.supports(IN_ARG_x_dot) && (x_dot = inArgs.get_x_dot()).get() )
00317     epetraInArgs.set_x_dot(get_Epetra_Vector(*x_map_,x_dot));
00318 
00319   RefCountPtr<const VectorBase<double> > x;
00320   if( inArgs.supports(IN_ARG_x) && (x = inArgs.get_x()).get() )
00321     epetraInArgs.set_x(get_Epetra_Vector(*x_map_,x));
00322 
00323   if(1) {
00324     RefCountPtr<const VectorBase<double> > p_l;
00325     for(int l = 0;  l < outArgs.Np(); ++l ) {
00326       p_l = inArgs.get_p(l);
00327       if(p_l.get()) epetraInArgs.set_p(l,get_Epetra_Vector(*p_map_[l],p_l));
00328     }
00329   }
00330 
00331   RefCountPtr<const Teuchos::Polynomial< VectorBase<double> > > x_dot_poly;
00332   Teuchos::RefCountPtr<Epetra_Vector> epetra_ptr;
00333   if( inArgs.supports(IN_ARG_x_dot_poly) && \
00334       (x_dot_poly = inArgs.get_x_dot_poly()).get() ) {
00335     RefCountPtr<Teuchos::Polynomial<Epetra_Vector> > epetra_x_dot_poly = 
00336       Teuchos::rcp(new Teuchos::Polynomial<Epetra_Vector>(x_dot_poly->degree()));
00337     for (unsigned int i=0; i<=x_dot_poly->degree(); i++) {
00338       epetra_ptr = 
00339   Teuchos::rcp_const_cast<Epetra_Vector>(get_Epetra_Vector(*x_map_, x_dot_poly->getCoefficient(i)));
00340       epetra_x_dot_poly->setCoefficientPtr(i,epetra_ptr);
00341     }
00342     epetraInArgs.set_x_dot_poly(epetra_x_dot_poly);
00343   }
00344 
00345   RefCountPtr<const Teuchos::Polynomial< VectorBase<double> > > x_poly;
00346   if( inArgs.supports(IN_ARG_x_poly) && \
00347       (x_poly = inArgs.get_x_poly()).get() ) {
00348     RefCountPtr<Teuchos::Polynomial<Epetra_Vector> > epetra_x_poly = 
00349       Teuchos::rcp(new Teuchos::Polynomial<Epetra_Vector>(x_poly->degree()));
00350     for (unsigned int i=0; i<=x_poly->degree(); i++) {
00351       epetra_ptr = 
00352   Teuchos::rcp_const_cast<Epetra_Vector>(get_Epetra_Vector(*x_map_, x_poly->getCoefficient(i)));
00353       epetra_x_poly->setCoefficientPtr(i,epetra_ptr);
00354     }
00355     epetraInArgs.set_x_poly(epetra_x_poly);
00356   }
00357 
00358   if( inArgs.supports(IN_ARG_t) )
00359     epetraInArgs.set_t(inArgs.get_t());
00360 
00361   if( inArgs.supports(IN_ARG_alpha) )
00362     epetraInArgs.set_alpha(inArgs.get_alpha());
00363 
00364   if( inArgs.supports(IN_ARG_beta) )
00365     epetraInArgs.set_beta(inArgs.get_beta());
00366 
00367   timer.stop();
00368   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00369     *OSTab(out).getOStream() << "\nTime to setup InArgs = "<<timer.totalElapsedTime()<<" sec\n";
00370 
00371   // OutArgs
00372   
00373   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00374     *out << "\nSetting-up/creating output arguments ...\n";
00375   timer.start(true);
00376 
00377   EME::OutArgs epetraOutArgs = epetraModel_->createOutArgs();
00378 
00379   RefCountPtr<VectorBase<double> > f;
00380   if( outArgs.supports(OUT_ARG_f) && (f = outArgs.get_f()).get() )
00381     epetraOutArgs.set_f(get_Epetra_Vector(*f_map_,f));
00382 
00383   if(1){
00384     Teuchos::RefCountPtr<VectorBase<double> > g_j;
00385     for(int j = 0;  j < outArgs.Ng(); ++j ) {
00386       g_j = outArgs.get_g(j);
00387       if(g_j.get()) epetraOutArgs.set_g(j,get_Epetra_Vector(*g_map_[j],g_j));
00388     }
00389   }
00390   
00391   RefCountPtr<LinearOpWithSolveBase<double> > W;
00392   RefCountPtr<LinearOpBase<double> >          W_op;
00393   RefCountPtr<const LinearOpBase<double> >    fwdW;
00394   RefCountPtr<EpetraLinearOp>                 efwdW;
00395   if( outArgs.supports(OUT_ARG_W) && (W = outArgs.get_W()).get() ) {
00396     Thyra::uninitializeOp<double>(*W_factory_,&*W,&fwdW);
00397     if(fwdW.get()) {
00398       efwdW = rcp_const_cast<EpetraLinearOp>(rcp_dynamic_cast<const EpetraLinearOp>(fwdW,true));
00399     }
00400     else {
00401       efwdW = Teuchos::rcp(new EpetraLinearOp());
00402       fwdW = efwdW;
00403     }
00404   }
00405   if( outArgs.supports(OUT_ARG_W_op) && (W_op = outArgs.get_W_op()).get() ) {
00406     if( W_op.get() && !efwdW.get() )
00407       efwdW = rcp_const_cast<EpetraLinearOp>(rcp_dynamic_cast<const EpetraLinearOp>(W_op,true));
00408   }
00409   RefCountPtr<Epetra_Operator> eW;
00410   if(efwdW.get()) {
00411     eW = efwdW->epetra_op();
00412     if(!eW.get())
00413       eW = epetraModel_->create_W();
00414     epetraOutArgs.set_W(eW);
00415   }
00416   // NOTE: Above, if both W and W_op are set and have been through at least
00417   // one prior evaluation (and therefore have Epetra_Operator objects embedded
00418   // in them), then we will use the Epetra_Operator embedded in W to pass to
00419   // the EpetraExt::ModelEvaluator object and ignore the Epetra_Operator
00420   // object in W_op.  In the standard use case, these will be the same
00421   // Epetra_Operator objects.  However, it is possible that the client could
00422   // use this interface in such a way that these would have different
00423   // Epetra_Operator objects embedded in them.  In this (very unlikely) case,
00424   // the Epetra_Operator embedded in W_op will be discarded!  This might be
00425   // surprising to a client but it is very unlikely that this will ever be a
00426   // problem, but the issue is duly noted here!  Only dangerous programming
00427   // use of this interface would cause any problem.
00428 
00429   if(1){
00430     Derivative<double> DfDp_l;
00431     for(int l = 0;  l < outArgs.Np(); ++l ) {
00432       if( !outArgs.supports(OUT_ARG_DfDp,l).none() && !(DfDp_l = outArgs.get_DfDp(l)).isEmpty() )
00433         epetraOutArgs.set_DfDp(l,convert(DfDp_l,f_map_,p_map_[l]));
00434     }
00435   }
00436 
00437   if(1){
00438     Derivative<double> DgDx_j;
00439     for(int j = 0;  j < outArgs.Ng(); ++j ) {
00440       if( !outArgs.supports(OUT_ARG_DgDx,j).none() && !(DgDx_j = outArgs.get_DgDx(j)).isEmpty() )
00441         epetraOutArgs.set_DgDx(j,convert(DgDx_j,g_map_[j],x_map_));
00442     }
00443   }
00444 
00445   if(1){
00446     Derivative<double> DgDp_j_l;
00447     for(int j = 0;  j < outArgs.Ng(); ++j ) {
00448       for(int l = 0;  l < outArgs.Np(); ++l ) {
00449         if( !outArgs.supports(OUT_ARG_DgDp,j,l).none() && !(DgDp_j_l = outArgs.get_DgDp(j,l)).isEmpty() )
00450           epetraOutArgs.set_DgDp(j,l,convert(DgDp_j_l,g_map_[j],p_map_[l]));
00451       }
00452     }
00453   }
00454 
00455   RefCountPtr<const Teuchos::Polynomial< VectorBase<double> > > f_poly;
00456   if( outArgs.supports(OUT_ARG_f_poly) && \
00457       (f_poly = outArgs.get_f_poly()).get() ) {
00458     RefCountPtr<Teuchos::Polynomial<Epetra_Vector> > epetra_f_poly = 
00459       Teuchos::rcp(new Teuchos::Polynomial<Epetra_Vector>(f_poly->degree()));
00460     for (unsigned int i=0; i<=f_poly->degree(); i++) {
00461       epetra_ptr = 
00462   Teuchos::rcp_const_cast<Epetra_Vector>(get_Epetra_Vector(*f_map_,
00463                  f_poly->getCoefficient(i)));
00464       epetra_f_poly->setCoefficientPtr(i,epetra_ptr);
00465     }
00466     epetraOutArgs.set_f_poly(epetra_f_poly);
00467   }
00468 
00469   timer.stop();
00470   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00471     *OSTab(out).getOStream() << "\nTime to setup OutArgs = "<<timer.totalElapsedTime()<<" sec\n";
00472 
00473   // Do the evaluation
00474   
00475   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00476     *out << "\nEvaluating the output functions ...\n";
00477   timer.start(true);
00478 
00479   epetraModel_->evalModel(epetraInArgs,epetraOutArgs);
00480 
00481   timer.stop();
00482   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00483     *OSTab(out).getOStream() << "\nTime to evaluate output functions = "<<timer.totalElapsedTime()<<" sec\n";
00484 
00485   // Postprocess arguments
00486   
00487   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00488     *out << "\nPost processing the output objects ...\n";
00489   timer.start(true);
00490 
00491   if(efwdW.get())
00492     efwdW->initialize(eW);  // This will directly update W_op if W.get()==NULL!
00493   
00494   if( W.get() ) {
00495     Thyra::initializeOp<double>(*W_factory_,fwdW,&*W);
00496     W->setOStream(this->getOStream());
00497   }
00498 
00499   if( W_op.get() ) {
00500     if( W_op.shares_resource(efwdW) ) {
00501       // W_op was already updated above since *efwdW is the same object as *W_op
00502     }
00503     else {
00504       rcp_dynamic_cast<EpetraLinearOp>(W_op,true)->initialize(eW);
00505     }
00506   }
00507 
00508   timer.stop();
00509   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00510     *OSTab(out).getOStream() << "\nTime to process output objects = "<<timer.totalElapsedTime()<<" sec\n";
00511 
00512   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_EXTREME))
00513     *out
00514       << "\noutArgs on output =\n" << Teuchos::describe(outArgs,verbLevel);
00515 
00516   totalTimer.stop();
00517   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00518     *out
00519       << "\nTotal evaluation time = "<<totalTimer.totalElapsedTime()<<" sec\n"
00520       << "\nLeaving Thyra::EpetraModelEvaluator::evalModel(...) ...\n";
00521   
00522 }
00523 
00524 void EpetraModelEvaluator::reportFinalPoint(
00525   const ModelEvaluatorBase::InArgs<double>      &finalPoint
00526   ,const bool                                   wasSolved
00527   )
00528 {
00529   finalPoint_ = this->createInArgs();
00530   finalPoint_.setArgs(finalPoint);
00531   finalPointWasSolved_ = wasSolved;
00532 }
00533 
00534 // Public functions overridden from Teuchos::Describable
00535 
00536 std::string EpetraModelEvaluator::description() const
00537 {
00538   std::ostringstream oss;
00539   oss << "Thyra::EpetraModelEvaluator{";
00540   oss << "epetraModel=";
00541   if(epetraModel_.get())
00542     oss << "\'"<<epetraModel_->description()<<"\'";
00543   else
00544     oss << "NULL";
00545   oss << ",W_factory=";
00546   if(W_factory_.get())
00547     oss << "\'"<<W_factory_->description()<<"\'";
00548   else
00549     oss << "NULL";
00550   oss << "}";
00551   return oss.str();
00552 }
00553 
00554 } // namespace Thyra
00555 
00556 //
00557 // Non-member utility functions
00558 //
00559 
00560 Thyra::ModelEvaluatorBase::EDerivativeMultiVectorOrientation
00561 Thyra::convert( const EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation &mvOrientation )
00562 {
00563   switch(mvOrientation) {
00564     case EpetraExt::ModelEvaluator::DERIV_MV_BY_COL :
00565       return ModelEvaluatorBase::DERIV_MV_BY_COL;
00566     case EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW :
00567       return ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW;
00568     default:
00569       TEST_FOR_EXCEPT(true);
00570   }
00571   return ModelEvaluatorBase::DERIV_MV_BY_COL; // Should never be called!
00572 }
00573 
00574 EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation
00575 Thyra::convert( const ModelEvaluatorBase::EDerivativeMultiVectorOrientation &mvOrientation )
00576 {
00577   switch(mvOrientation) {
00578     case ModelEvaluatorBase::DERIV_MV_BY_COL :
00579       return EpetraExt::ModelEvaluator::DERIV_MV_BY_COL;
00580     case ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW :
00581       return EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW;
00582     default:
00583       TEST_FOR_EXCEPT(true);
00584   }
00585   return EpetraExt::ModelEvaluator::DERIV_MV_BY_COL; // Should never be called!
00586 }
00587 
00588 Thyra::ModelEvaluatorBase::DerivativeProperties
00589 Thyra::convert( const EpetraExt::ModelEvaluator::DerivativeProperties &derivativeProperties )
00590 {
00591   ModelEvaluatorBase::EDerivativeLinearity linearity;
00592   switch(derivativeProperties.linearity) {
00593     case EpetraExt::ModelEvaluator::DERIV_LINEARITY_UNKNOWN:
00594       linearity = ModelEvaluatorBase::DERIV_LINEARITY_UNKNOWN;
00595       break;
00596     case EpetraExt::ModelEvaluator::DERIV_LINEARITY_CONST:
00597       linearity = ModelEvaluatorBase::DERIV_LINEARITY_CONST;
00598       break;
00599     case  EpetraExt::ModelEvaluator::DERIV_LINEARITY_NONCONST:
00600       linearity = ModelEvaluatorBase::DERIV_LINEARITY_NONCONST;
00601       break;
00602     default:
00603       TEST_FOR_EXCEPT(true);
00604   }
00605   ModelEvaluatorBase::ERankStatus rank;
00606   switch(derivativeProperties.rank) {
00607     case EpetraExt::ModelEvaluator::DERIV_RANK_UNKNOWN:
00608       rank = ModelEvaluatorBase::DERIV_RANK_UNKNOWN;
00609       break;
00610     case EpetraExt::ModelEvaluator::DERIV_RANK_FULL:
00611       rank = ModelEvaluatorBase::DERIV_RANK_FULL;
00612       break;
00613     case EpetraExt::ModelEvaluator::DERIV_RANK_DEFICIENT:
00614       rank = ModelEvaluatorBase::DERIV_RANK_DEFICIENT;
00615       break;
00616     default:
00617       TEST_FOR_EXCEPT(true);
00618   }
00619   return ModelEvaluatorBase::DerivativeProperties(linearity,rank,derivativeProperties.supportsAdjoint);
00620 }
00621 
00622 Thyra::ModelEvaluatorBase::DerivativeSupport
00623 Thyra::convert( const EpetraExt::ModelEvaluator::DerivativeSupport &derivativeSupport )
00624 {
00625   ModelEvaluatorBase::DerivativeSupport ds;
00626   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_LINEAR_OP))
00627     ds.plus(ModelEvaluatorBase::DERIV_LINEAR_OP);
00628   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_MV_BY_COL))
00629     ds.plus(ModelEvaluatorBase::DERIV_MV_BY_COL);
00630   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW))
00631     ds.plus(ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW);
00632   return ds;
00633 }
00634 
00635 EpetraExt::ModelEvaluator::Derivative
00636 Thyra::convert(
00637   const ModelEvaluatorBase::Derivative<double>        &derivative
00638   ,const Teuchos::RefCountPtr<const Epetra_Map>       &fnc_map
00639   ,const Teuchos::RefCountPtr<const Epetra_Map>       &var_map
00640   )
00641 {
00642   typedef ModelEvaluatorBase MEB;
00643   if(derivative.getLinearOp().get()) {
00644     return EpetraExt::ModelEvaluator::Derivative(
00645       Teuchos::rcp_const_cast<Epetra_Operator>(
00646         Teuchos::dyn_cast<const EpetraLinearOp>(*derivative.getLinearOp()).epetra_op()
00647         )
00648       );
00649   }
00650   else if(derivative.getDerivativeMultiVector().getMultiVector().get()) {
00651     return EpetraExt::ModelEvaluator::Derivative(
00652       EpetraExt::ModelEvaluator::DerivativeMultiVector(
00653         get_Epetra_MultiVector(
00654           derivative.getDerivativeMultiVector().getOrientation() == MEB::DERIV_MV_BY_COL
00655           ? *fnc_map : *var_map
00656           ,derivative.getDerivativeMultiVector().getMultiVector()
00657           )
00658         ,convert(derivative.getDerivativeMultiVector().getOrientation())
00659         )
00660       );
00661   }
00662   return EpetraExt::ModelEvaluator::Derivative();
00663 }
00664 

Generated on Thu Sep 18 12:31:51 2008 for EpetraExt/Thyra Adapters by doxygen 1.3.9.1