Teko Version of the Day
Teko_EpetraOperatorWrapper.cpp
00001 // @HEADER
00002 // 
00003 // ***********************************************************************
00004 // 
00005 //      Teko: A package for block and physics based preconditioning
00006 //                  Copyright 2010 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 Eric C. Cyr (eccyr@sandia.gov)
00039 // 
00040 // ***********************************************************************
00041 // 
00042 // @HEADER
00043 
00044 
00045 #include "Teko_EpetraOperatorWrapper.hpp"
00046 #include "Thyra_SpmdVectorBase.hpp"
00047 #include "Thyra_MultiVectorStdOps.hpp"
00048 #ifdef HAVE_MPI
00049 #include "Teuchos_DefaultMpiComm.hpp"
00050 #endif
00051 #include "Teuchos_DefaultSerialComm.hpp"
00052 #include "Thyra_EpetraLinearOp.hpp"
00053 #include "Epetra_SerialComm.h"
00054 #include "Epetra_Vector.h"
00055 #ifdef HAVE_MPI
00056 #include "Epetra_MpiComm.h"
00057 #endif
00058 #include "Thyra_EpetraThyraWrappers.hpp"
00059 
00060 // #include "Thyra_LinearOperator.hpp"
00061 #include "Thyra_BlockedLinearOpBase.hpp"
00062 #include "Thyra_ProductVectorSpaceBase.hpp"
00063 #include "Thyra_get_Epetra_Operator.hpp"
00064 
00065 #include "Teko_EpetraThyraConverter.hpp"
00066 #include "Teuchos_Ptr.hpp"
00067 
00068 
00069 namespace Teko { 
00070 namespace Epetra { 
00071 
00072 
00073 using namespace Teuchos;
00074 using namespace Thyra;
00075 
00076 DefaultMappingStrategy::DefaultMappingStrategy(const RCP<const Thyra::LinearOpBase<double> > & thyraOp,const Epetra_Comm & comm)
00077 {
00078    RCP<Epetra_Comm> newComm = rcp(comm.Clone());
00079 
00080    // extract vector spaces from linear operator
00081    domainSpace_ = thyraOp->domain();
00082    rangeSpace_ = thyraOp->range();
00083 
00084    domainMap_ = Teko::Epetra::thyraVSToEpetraMap(*domainSpace_,newComm);
00085    rangeMap_ = Teko::Epetra::thyraVSToEpetraMap(*rangeSpace_,newComm);
00086 }
00087 
00088 void DefaultMappingStrategy::copyEpetraIntoThyra(const Epetra_MultiVector& x, const Ptr<Thyra::MultiVectorBase<double> > & thyraVec) const
00089 {
00090    Teko::Epetra::blockEpetraToThyra(x,thyraVec);
00091 }
00092 
00093 void DefaultMappingStrategy::copyThyraIntoEpetra(const RCP<const Thyra::MultiVectorBase<double> > & thyraVec, Epetra_MultiVector& v) const
00094 {
00095    Teko::Epetra::blockThyraToEpetra(thyraVec,v);
00096 }
00097 
00098 EpetraOperatorWrapper::EpetraOperatorWrapper()
00099 {
00100    useTranspose_ = false;
00101    mapStrategy_ = Teuchos::null;
00102    thyraOp_ = Teuchos::null;
00103    comm_ = Teuchos::null;
00104    label_ = Teuchos::null;
00105 }
00106 
00107 EpetraOperatorWrapper::EpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<double> > & thyraOp)
00108 { 
00109    SetOperator(thyraOp); 
00110 }
00111 
00112 EpetraOperatorWrapper::EpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<double> > & thyraOp,const RCP<const MappingStrategy> & mapStrategy)
00113    : mapStrategy_(mapStrategy)
00114 {
00115    SetOperator(thyraOp);
00116 }
00117 
00118 EpetraOperatorWrapper::EpetraOperatorWrapper(const RCP<const MappingStrategy> & mapStrategy)
00119    : mapStrategy_(mapStrategy)
00120 {
00121    useTranspose_ = false;
00122    thyraOp_ = Teuchos::null;
00123    comm_ = Teuchos::null;
00124    label_ = Teuchos::null;
00125 }
00126 
00127 void EpetraOperatorWrapper::SetOperator(const RCP<const LinearOpBase<double> > & thyraOp,bool buildMap)
00128 {
00129    useTranspose_ = false;
00130    thyraOp_ = thyraOp;
00131    comm_ = getEpetraComm(*thyraOp);
00132    label_ = thyraOp_->description();
00133    if(mapStrategy_==Teuchos::null && buildMap)
00134       mapStrategy_ = Teuchos::rcp(new DefaultMappingStrategy(thyraOp,*comm_));
00135 }
00136 
00137 double EpetraOperatorWrapper::NormInf() const 
00138 {
00139   TEST_FOR_EXCEPTION(true, std::runtime_error,
00140                      "EpetraOperatorWrapper::NormInf not implemated");
00141   return 1.0;
00142 }
00143 
00144 int EpetraOperatorWrapper::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00145 {
00146    if (!useTranspose_)
00147    {
00148        // allocate space for each vector
00149        RCP<Thyra::MultiVectorBase<double> > tX;
00150        RCP<Thyra::MultiVectorBase<double> > tY; 
00151 
00152        tX = Thyra::createMembers(thyraOp_->domain(),X.NumVectors()); 
00153        tY = Thyra::createMembers(thyraOp_->range(),X.NumVectors());
00154 
00155        Thyra::assign(tX.ptr(),0.0);
00156        Thyra::assign(tY.ptr(),0.0);
00157 
00158        // copy epetra X into thyra X
00159        mapStrategy_->copyEpetraIntoThyra(X, tX.ptr());
00160        mapStrategy_->copyEpetraIntoThyra(Y, tY.ptr()); // if this matrix isn't block square, this probably won't work!
00161 
00162        // perform matrix vector multiplication
00163        thyraOp_->apply(Thyra::NOTRANS,*tX,tY.ptr(),1.0,0.0);
00164 
00165        // copy thyra Y into epetra Y
00166        mapStrategy_->copyThyraIntoEpetra(tY, Y);
00167    }
00168    else
00169    {
00170        TEUCHOS_ASSERT(false);
00171    }
00172  
00173    return 0;
00174 }
00175 
00176 
00177 int EpetraOperatorWrapper::ApplyInverse(const Epetra_MultiVector& X, 
00178                                       Epetra_MultiVector& Y) const
00179 {
00180   TEST_FOR_EXCEPTION(true, std::runtime_error,
00181                      "EpetraOperatorWrapper::ApplyInverse not implemented");
00182   return 1;
00183 }
00184 
00185 
00186 RCP<const Epetra_Comm> 
00187 EpetraOperatorWrapper::getEpetraComm(const Thyra::LinearOpBase<double>& inOp) const
00188 {
00189   RCP<const VectorSpaceBase<double> > vs = inOp.domain();
00190 
00191   RCP<const SpmdVectorSpaceBase<double> > spmd;
00192   RCP<const VectorSpaceBase<double> > current = vs;
00193   while(current!=Teuchos::null) {
00194      // try to cast to a product vector space first
00195      RCP<const ProductVectorSpaceBase<double> > prod
00196            = rcp_dynamic_cast<const ProductVectorSpaceBase<double> >(current);
00197 
00198      // figure out what type it is
00199      if(prod==Teuchos::null) {
00200         // hopfully this is a SPMD vector space
00201         spmd = rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(current);
00202 
00203         break;
00204      }
00205      else // get first convenient vector space
00206         current = prod->getBlock(0);
00207   }
00208 
00209   TEST_FOR_EXCEPTION(spmd==Teuchos::null, std::runtime_error, 
00210                      "EpetraOperatorWrapper requires std::vector space "
00211                      "blocks to be SPMD std::vector spaces");
00212 
00213   return Thyra::get_Epetra_Comm(*spmd->getComm());
00214 /*
00215   const Thyra::ConstLinearOperator<double> thyraOp = rcpFromRef(inOp); 
00216 
00217   RCP<Epetra_Comm> rtn;
00218   // VectorSpace<double> vs = thyraOp.domain().getBlock(0);
00219   RCP<const VectorSpaceBase<double> > vs = thyraOp.domain().getBlock(0).constPtr();
00220 
00221   // search for an SpmdVectorSpaceBase object
00222   RCP<const SpmdVectorSpaceBase<double> > spmd;
00223   RCP<const VectorSpaceBase<double> > current = vs;
00224   while(current!=Teuchos::null) {
00225      // try to cast to a product vector space first
00226      RCP<const ProductVectorSpaceBase<double> > prod
00227            = rcp_dynamic_cast<const ProductVectorSpaceBase<double> >(current);
00228 
00229      // figure out what type it is
00230      if(prod==Teuchos::null) {
00231         // hopfully this is a SPMD vector space
00232         spmd = rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(current);
00233 
00234         break;
00235      }
00236      else {
00237         // get first convenient vector space
00238         current = prod->getBlock(0);
00239      }
00240   }
00241 
00242   TEST_FOR_EXCEPTION(spmd==Teuchos::null, std::runtime_error, 
00243                      "EpetraOperatorWrapper requires std::vector space "
00244                      "blocks to be SPMD std::vector spaces");
00245 
00246   const SerialComm<Thyra::Ordinal>* serialComm 
00247     = dynamic_cast<const SerialComm<Thyra::Ordinal>*>(spmd->getComm().get());
00248 
00249 #ifdef HAVE_MPI
00250   const MpiComm<Thyra::Ordinal>* mpiComm 
00251     = dynamic_cast<const MpiComm<Thyra::Ordinal>*>(spmd->getComm().get());
00252 
00253   TEST_FOR_EXCEPTION(mpiComm==0 && serialComm==0, std::runtime_error, 
00254                      "SPMD std::vector space has a communicator that is "
00255                      "neither a serial comm nor an MPI comm");
00256 
00257   if (mpiComm != 0)
00258     {
00259       rtn = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
00260     }
00261   else
00262     {
00263       rtn = rcp(new Epetra_SerialComm());
00264     }
00265 #else
00266   TEST_FOR_EXCEPTION(serialComm==0, std::runtime_error, 
00267                      "SPMD std::vector space has a communicator that is "
00268                      "neither a serial comm nor an MPI comm");
00269   rtn = rcp(new Epetra_SerialComm());
00270   
00271 #endif
00272 
00273   TEST_FOR_EXCEPTION(rtn.get()==0, std::runtime_error, "null communicator created");
00274   return rtn;
00275 */
00276 }
00277 
00278 int EpetraOperatorWrapper::GetBlockRowCount()
00279 {
00280    const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp 
00281          = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(getThyraOp());
00282 
00283    return blkOp->productRange()->numBlocks();
00284 }
00285 
00286 int EpetraOperatorWrapper::GetBlockColCount()
00287 {
00288    const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp 
00289          = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(getThyraOp());
00290 
00291    return blkOp->productDomain()->numBlocks();
00292 }
00293 
00294 Teuchos::RCP<const Epetra_Operator> EpetraOperatorWrapper::GetBlock(int i,int j) const
00295 {
00296    const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp 
00297          = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(getThyraOp());
00298 
00299    return Thyra::get_Epetra_Operator(*blkOp->getBlock(i,j));
00300 }
00301 
00302 } // namespace Epetra
00303 } // namespace Teko
 All Classes Files Functions Variables