Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_EpetraOperatorWrapper.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //              Meros: Segregated Preconditioning Package
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_EpetraOperatorWrapper.hpp"
00030 #include "Thyra_EpetraThyraWrappers.hpp"
00031 #include "Thyra_DetachedSpmdVectorView.hpp"
00032 #include "Thyra_DefaultProductVector.hpp"
00033 #include "Thyra_ProductVectorSpaceBase.hpp"
00034 #include "Thyra_SpmdVectorBase.hpp"
00035 #include "Thyra_EpetraLinearOp.hpp"
00036 
00037 #ifdef HAVE_MPI
00038 #  include "Epetra_MpiComm.h"
00039 #endif
00040 #include "Epetra_SerialComm.h"
00041 #include "Epetra_Vector.h"
00042 
00043 #ifdef HAVE_MPI
00044 #  include "Teuchos_DefaultMpiComm.hpp"
00045 #endif
00046 #include "Teuchos_DefaultSerialComm.hpp"
00047 
00048 
00049 namespace Thyra { 
00050 
00051 
00052 // Constructor, utilties
00053 
00054 
00055 EpetraOperatorWrapper::EpetraOperatorWrapper(
00056   const RCP<const LinearOpBase<double> > &thyraOp
00057   )
00058   : useTranspose_(false),
00059     thyraOp_(thyraOp),
00060     range_(thyraOp->range()),
00061     domain_(thyraOp->domain()),
00062     comm_(getEpetraComm(*thyraOp)),
00063     rangeMap_(get_Epetra_Map(*range_, comm_)),
00064     domainMap_(get_Epetra_Map(*domain_, comm_)),
00065     label_(thyraOp->description())
00066 {;}
00067 
00068 
00069 void EpetraOperatorWrapper::copyEpetraIntoThyra(const Epetra_MultiVector& x,
00070   const Ptr<VectorBase<double> > &thyraVec) const
00071 {
00072 
00073   using Teuchos::rcpFromPtr;
00074   using Teuchos::rcp_dynamic_cast;
00075 
00076   const int numVecs = x.NumVectors();
00077 
00078   TEST_FOR_EXCEPTION(numVecs != 1, std::runtime_error,
00079     "epetraToThyra does not work with MV dimension != 1");
00080 
00081   const RCP<ProductVectorBase<double> > prodThyraVec =
00082     castOrCreateNonconstProductVectorBase(rcpFromPtr(thyraVec));
00083 
00084   const ArrayView<const double> epetraData(x[0], x.Map().NumMyElements());
00085   // NOTE: I tried using Epetra_MultiVector::operator()(int) to return an
00086   // Epetra_Vector object but it has a defect when Reset(...) is called which
00087   // results in a memory access error (see bug 4700).
00088 
00089   int offset = 0;
00090   const int numBlocks = prodThyraVec->productSpace()->numBlocks();
00091   for (int b = 0; b < numBlocks; ++b) {
00092     const RCP<VectorBase<double> > vec_b = prodThyraVec->getNonconstVectorBlock(b);
00093     const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_b =
00094       rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(vec_b->space(), true);
00095     DetachedSpmdVectorView<double> view(vec_b);
00096     const ArrayRCP<double> thyraData = view.sv().values();
00097     const int localNumElems = spmd_vs_b->localSubDim();
00098     for (int i=0; i < localNumElems; ++i) {
00099       thyraData[i] = epetraData[i+offset];
00100     }
00101     offset += localNumElems;
00102   }
00103 
00104 }
00105 
00106 
00107 void EpetraOperatorWrapper::copyThyraIntoEpetra(
00108   const VectorBase<double>& thyraVec, Epetra_MultiVector& x) const 
00109 {
00110 
00111   using Teuchos::rcpFromRef;
00112   using Teuchos::rcp_dynamic_cast;
00113 
00114   const int numVecs = x.NumVectors();
00115 
00116   TEST_FOR_EXCEPTION(numVecs != 1, std::runtime_error,
00117     "epetraToThyra does not work with MV dimension != 1");
00118 
00119   const RCP<const ProductVectorBase<double> > prodThyraVec =
00120     castOrCreateProductVectorBase(rcpFromRef(thyraVec));
00121 
00122   const ArrayView<double> epetraData(x[0], x.Map().NumMyElements());
00123   // NOTE: See above!
00124 
00125   int offset = 0;
00126   const int numBlocks = prodThyraVec->productSpace()->numBlocks();
00127   for (int b = 0; b < numBlocks; ++b) {
00128     const RCP<const VectorBase<double> > vec_b = prodThyraVec->getVectorBlock(b);
00129     const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_b =
00130       rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(vec_b->space(), true);
00131     ConstDetachedSpmdVectorView<double> view(vec_b);
00132     const ArrayRCP<const double> thyraData = view.sv().values();
00133     const int localNumElems = spmd_vs_b->localSubDim();
00134     for (int i=0; i < localNumElems; ++i) {
00135       epetraData[i+offset] = thyraData[i];
00136     }
00137     offset += localNumElems;
00138   }
00139 
00140 }
00141 
00142 
00143 // Overridden from Epetra_Operator
00144 
00145 
00146 int EpetraOperatorWrapper::Apply(const Epetra_MultiVector& X,
00147   Epetra_MultiVector& Y) const
00148 {
00149 
00150   const RCP<const VectorSpaceBase<double> >
00151     opRange = ( !useTranspose_ ? range_ : domain_ ),
00152     opDomain = ( !useTranspose_ ? domain_ : range_ );
00153 
00154   const RCP<VectorBase<double> > tx = createMember(opDomain);
00155   copyEpetraIntoThyra(X, tx.ptr());
00156 
00157   const RCP<VectorBase<double> > ty = createMember(opRange);
00158 
00159   Thyra::apply<double>( *thyraOp_, !useTranspose_ ? NOTRANS : CONJTRANS, *tx, ty.ptr());
00160 
00161   copyThyraIntoEpetra(*ty, Y);
00162 
00163   return 0;
00164 
00165 }
00166 
00167 
00168 int EpetraOperatorWrapper::ApplyInverse(const Epetra_MultiVector& X, 
00169   Epetra_MultiVector& Y) const
00170 {
00171   TEST_FOR_EXCEPTION(true, std::runtime_error,
00172     "EpetraOperatorWrapper::ApplyInverse not implemented");
00173   return 1;
00174 }
00175 
00176 
00177 double EpetraOperatorWrapper::NormInf() const 
00178 {
00179   TEST_FOR_EXCEPTION(true, std::runtime_error,
00180     "EpetraOperatorWrapper::NormInf not implemated");
00181   return 1.0;
00182 }
00183 
00184 
00185 // private
00186 
00187 
00188 RCP<const Epetra_Comm> 
00189 EpetraOperatorWrapper::getEpetraComm(const LinearOpBase<double>& thyraOp)
00190 {
00191 
00192   using Teuchos::rcp;
00193   using Teuchos::rcp_dynamic_cast;
00194   using Teuchos::SerialComm;
00195 #ifdef HAVE_MPI
00196   using Teuchos::MpiComm;
00197 #endif
00198 
00199   RCP<const VectorSpaceBase<double> > vs = thyraOp.range();
00200   
00201   const RCP<const ProductVectorSpaceBase<double> > prod_vs =
00202     rcp_dynamic_cast<const ProductVectorSpaceBase<double> >(vs);
00203 
00204   if (nonnull(prod_vs))
00205     vs = prod_vs->getBlock(0);
00206 
00207   const RCP<const SpmdVectorSpaceBase<double> > spmd_vs =
00208     rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(vs, true);
00209 
00210   return get_Epetra_Comm(*spmd_vs->getComm());
00211 
00212 }
00213 
00214 
00215 } // namespace Thyra
00216 
00217 
00218 Teuchos::RCP<const Thyra::LinearOpBase<double> > 
00219 Thyra::makeEpetraWrapper(const RCP<const LinearOpBase<double> > &thyraOp)
00220 {
00221   return epetraLinearOp(
00222     Teuchos::rcp(new EpetraOperatorWrapper(thyraOp))
00223     );
00224 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines