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 //    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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include "Thyra_EpetraOperatorWrapper.hpp"
00043 #include "Thyra_EpetraThyraWrappers.hpp"
00044 #include "Thyra_DetachedSpmdVectorView.hpp"
00045 #include "Thyra_DefaultProductVector.hpp"
00046 #include "Thyra_ProductVectorSpaceBase.hpp"
00047 #include "Thyra_SpmdVectorBase.hpp"
00048 #include "Thyra_EpetraLinearOp.hpp"
00049 
00050 #ifdef HAVE_MPI
00051 #  include "Epetra_MpiComm.h"
00052 #endif
00053 #include "Epetra_SerialComm.h"
00054 #include "Epetra_Vector.h"
00055 
00056 #ifdef HAVE_MPI
00057 #  include "Teuchos_DefaultMpiComm.hpp"
00058 #endif
00059 #include "Teuchos_DefaultSerialComm.hpp"
00060 
00061 
00062 namespace Thyra { 
00063 
00064 
00065 // Constructor, utilties
00066 
00067 
00068 EpetraOperatorWrapper::EpetraOperatorWrapper(
00069   const RCP<const LinearOpBase<double> > &thyraOp
00070   )
00071   : useTranspose_(false),
00072     thyraOp_(thyraOp),
00073     range_(thyraOp->range()),
00074     domain_(thyraOp->domain()),
00075     comm_(getEpetraComm(*thyraOp)),
00076     rangeMap_(get_Epetra_Map(*range_, comm_)),
00077     domainMap_(get_Epetra_Map(*domain_, comm_)),
00078     label_(thyraOp->description())
00079 {;}
00080 
00081 
00082 void EpetraOperatorWrapper::copyEpetraIntoThyra(const Epetra_MultiVector& x,
00083   const Ptr<VectorBase<double> > &thyraVec) const
00084 {
00085 
00086   using Teuchos::rcpFromPtr;
00087   using Teuchos::rcp_dynamic_cast;
00088 
00089   const int numVecs = x.NumVectors();
00090 
00091   TEUCHOS_TEST_FOR_EXCEPTION(numVecs != 1, std::runtime_error,
00092     "epetraToThyra does not work with MV dimension != 1");
00093 
00094   const RCP<ProductVectorBase<double> > prodThyraVec =
00095     castOrCreateNonconstProductVectorBase(rcpFromPtr(thyraVec));
00096 
00097   const ArrayView<const double> epetraData(x[0], x.Map().NumMyElements());
00098   // NOTE: I tried using Epetra_MultiVector::operator()(int) to return an
00099   // Epetra_Vector object but it has a defect when Reset(...) is called which
00100   // results in a memory access error (see bug 4700).
00101 
00102   int offset = 0;
00103   const int numBlocks = prodThyraVec->productSpace()->numBlocks();
00104   for (int b = 0; b < numBlocks; ++b) {
00105     const RCP<VectorBase<double> > vec_b = prodThyraVec->getNonconstVectorBlock(b);
00106     const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_b =
00107       rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(vec_b->space(), true);
00108     DetachedSpmdVectorView<double> view(vec_b);
00109     const ArrayRCP<double> thyraData = view.sv().values();
00110     const int localNumElems = spmd_vs_b->localSubDim();
00111     for (int i=0; i < localNumElems; ++i) {
00112       thyraData[i] = epetraData[i+offset];
00113     }
00114     offset += localNumElems;
00115   }
00116 
00117 }
00118 
00119 
00120 void EpetraOperatorWrapper::copyThyraIntoEpetra(
00121   const VectorBase<double>& thyraVec, Epetra_MultiVector& x) const 
00122 {
00123 
00124   using Teuchos::rcpFromRef;
00125   using Teuchos::rcp_dynamic_cast;
00126 
00127   const int numVecs = x.NumVectors();
00128 
00129   TEUCHOS_TEST_FOR_EXCEPTION(numVecs != 1, std::runtime_error,
00130     "epetraToThyra does not work with MV dimension != 1");
00131 
00132   const RCP<const ProductVectorBase<double> > prodThyraVec =
00133     castOrCreateProductVectorBase(rcpFromRef(thyraVec));
00134 
00135   const ArrayView<double> epetraData(x[0], x.Map().NumMyElements());
00136   // NOTE: See above!
00137 
00138   int offset = 0;
00139   const int numBlocks = prodThyraVec->productSpace()->numBlocks();
00140   for (int b = 0; b < numBlocks; ++b) {
00141     const RCP<const VectorBase<double> > vec_b = prodThyraVec->getVectorBlock(b);
00142     const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_b =
00143       rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(vec_b->space(), true);
00144     ConstDetachedSpmdVectorView<double> view(vec_b);
00145     const ArrayRCP<const double> thyraData = view.sv().values();
00146     const int localNumElems = spmd_vs_b->localSubDim();
00147     for (int i=0; i < localNumElems; ++i) {
00148       epetraData[i+offset] = thyraData[i];
00149     }
00150     offset += localNumElems;
00151   }
00152 
00153 }
00154 
00155 
00156 // Overridden from Epetra_Operator
00157 
00158 
00159 int EpetraOperatorWrapper::Apply(const Epetra_MultiVector& X,
00160   Epetra_MultiVector& Y) const
00161 {
00162 
00163   const RCP<const VectorSpaceBase<double> >
00164     opRange = ( !useTranspose_ ? range_ : domain_ ),
00165     opDomain = ( !useTranspose_ ? domain_ : range_ );
00166 
00167   const RCP<VectorBase<double> > tx = createMember(opDomain);
00168   copyEpetraIntoThyra(X, tx.ptr());
00169 
00170   const RCP<VectorBase<double> > ty = createMember(opRange);
00171 
00172   Thyra::apply<double>( *thyraOp_, !useTranspose_ ? NOTRANS : CONJTRANS, *tx, ty.ptr());
00173 
00174   copyThyraIntoEpetra(*ty, Y);
00175 
00176   return 0;
00177 
00178 }
00179 
00180 
00181 int EpetraOperatorWrapper::ApplyInverse(const Epetra_MultiVector& X, 
00182   Epetra_MultiVector& Y) const
00183 {
00184   TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00185     "EpetraOperatorWrapper::ApplyInverse not implemented");
00186   return 1;
00187 }
00188 
00189 
00190 double EpetraOperatorWrapper::NormInf() const 
00191 {
00192   TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00193     "EpetraOperatorWrapper::NormInf not implemated");
00194   return 1.0;
00195 }
00196 
00197 
00198 // private
00199 
00200 
00201 RCP<const Epetra_Comm> 
00202 EpetraOperatorWrapper::getEpetraComm(const LinearOpBase<double>& thyraOp)
00203 {
00204 
00205   using Teuchos::rcp;
00206   using Teuchos::rcp_dynamic_cast;
00207   using Teuchos::SerialComm;
00208 #ifdef HAVE_MPI
00209   using Teuchos::MpiComm;
00210 #endif
00211 
00212   RCP<const VectorSpaceBase<double> > vs = thyraOp.range();
00213   
00214   const RCP<const ProductVectorSpaceBase<double> > prod_vs =
00215     rcp_dynamic_cast<const ProductVectorSpaceBase<double> >(vs);
00216 
00217   if (nonnull(prod_vs))
00218     vs = prod_vs->getBlock(0);
00219 
00220   const RCP<const SpmdVectorSpaceBase<double> > spmd_vs =
00221     rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(vs, true);
00222 
00223   return get_Epetra_Comm(*spmd_vs->getComm());
00224 
00225 }
00226 
00227 
00228 } // namespace Thyra
00229 
00230 
00231 Teuchos::RCP<const Thyra::LinearOpBase<double> > 
00232 Thyra::makeEpetraWrapper(const RCP<const LinearOpBase<double> > &thyraOp)
00233 {
00234   return epetraLinearOp(
00235     Teuchos::rcp(new EpetraOperatorWrapper(thyraOp))
00236     );
00237 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines