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_SpmdVectorBase.hpp"
00031 #ifdef HAVE_MPI
00032 #include "Teuchos_DefaultMpiComm.hpp"
00033 #endif
00034 #include "Teuchos_DefaultSerialComm.hpp"
00035 #include "Thyra_EpetraLinearOp.hpp"
00036 #include "Epetra_SerialComm.h"
00037 #include "Epetra_Vector.h"
00038 #ifdef HAVE_MPI
00039 #include "Epetra_MpiComm.h"
00040 #endif
00041 #include "Thyra_EpetraThyraWrappers.hpp"
00042 
00043 using namespace Thyra;
00044 using namespace Teuchos;
00045 
00046 EpetraOperatorWrapper::EpetraOperatorWrapper(const ConstLinearOperator<double>& thyraOp)
00047   : useTranspose_(false),
00048     thyraOp_(thyraOp),
00049     domain_(thyraOp_.domain()),
00050     range_(thyraOp_.range()),
00051     comm_(getEpetraComm(thyraOp)),
00052     domainMap_(thyraVSToEpetraMap(thyraOp_.domain(), comm_)),
00053     rangeMap_(thyraVSToEpetraMap(thyraOp_.range(), comm_)),
00054     label_(thyraOp_.description())
00055 {;}
00056 
00057 double EpetraOperatorWrapper::NormInf() const 
00058 {
00059   TEST_FOR_EXCEPTION(true, runtime_error,
00060                      "EpetraOperatorWrapper::NormInf not implemated");
00061   return 1.0;
00062 }
00063 
00064 RefCountPtr<Epetra_Map> 
00065 EpetraOperatorWrapper
00066 ::thyraVSToEpetraMap(const VectorSpace<double>& vs,
00067                      const RefCountPtr<Epetra_Comm>& comm) const 
00068 {
00069   int globalDim = vs.dim();
00070   int myLocalElements = 0;
00071   
00072   /* find the number of local elements, summed over all blocks */
00073   if (isSPMD(vs))
00074     {
00075       myLocalElements = numLocalElements(vs);
00076     }
00077   else
00078     {
00079       for (int b=0; b<vs.numBlocks(); b++)
00080         {
00081           TEST_FOR_EXCEPTION(!isSPMD(vs.getBlock(b)), runtime_error, 
00082                              "EpetraOperatorWrapper requires vector space "
00083                              "blocks to be SPMD vectors");
00084           myLocalElements += numLocalElements(vs.getBlock(b));
00085         }
00086     }
00087 
00088   
00089   /* find the GIDs owned by this processor, taken from all blocks */
00090   Array<int> myGIDs(myLocalElements);
00091   
00092   int count=0;
00093   for (int b=0; b<vs.numBlocks(); b++)
00094     {
00095       int lowGIDInBlock = lowestLocallyOwnedIndex(vs.getBlock(b));
00096       int numLocalElementsInBlock = numLocalElements(vs.getBlock(b));
00097       for (int i=0; i<numLocalElementsInBlock; i++, count++)
00098         {
00099           myGIDs[count] = lowGIDInBlock+i;
00100         }
00101     }
00102 
00103   /* create the map */
00104   RefCountPtr<Epetra_Map> rtn 
00105     = rcp(new Epetra_Map(globalDim, myLocalElements, &(myGIDs[0]), 0, *comm));
00106 
00107   return rtn;
00108 }
00109 
00110 
00111 void EpetraOperatorWrapper::copyEpetraIntoThyra(const Epetra_MultiVector& x,
00112                                                 Vector<double> thyraVec) const
00113 {
00114   int numVecs = x.NumVectors();
00115   TEST_FOR_EXCEPTION(numVecs != 1, runtime_error,
00116                      "epetraToThyra does not work with MV dimension != 1");
00117 
00118   double* const epetraData = x[0];
00119 
00120   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00121     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00122   bool verbose = false;
00123   if (verbose)
00124     {
00125       *out << "in ETO::copyEpetraIntoThyra()" << endl;
00126       double normIn=0;
00127       x.Norm2(&normIn);
00128       *out << "before: norm(Epetra vec) = " << normIn << endl;
00129     }
00130 
00131   int offset = 0;
00132   for (int b=0; b<thyraVec.numBlocks(); b++)
00133     {
00134       Vector<double> block = thyraVec.getBlock(b);
00135       TEST_FOR_EXCEPTION(!isSPMD(Thyra::space(block)), runtime_error, "block is not SPMD");
00136       int localOffset = lowestLocallyOwnedIndex(Thyra::space(block));
00137       int localNumElems = numLocalElements(Thyra::space(block));
00138 
00139       RefCountPtr<SpmdVectorBase<double> > spmd 
00140         = rcp_dynamic_cast<SpmdVectorBase<double> >(block.ptr());
00141 
00143       {
00144         /* limit the scope so that it gets destroyed, and thus committed,
00145          * before using the vector */
00146         Teuchos::RefCountPtr<DetachedVectorView<double> > view 
00147           = rcp(new DetachedVectorView<double>(block.ptr(), 
00148                                                Thyra::Range1D(localOffset,
00149                                                        localOffset+localNumElems-1), 
00150                                                true));
00151         double* thyraData = view->values();
00152         for (int i=0; i<localNumElems; i++)
00153           {
00154             thyraData[i] = epetraData[i+offset];
00155           }
00156       }
00157       offset += localNumElems;
00158     }
00159   if (verbose)
00160     {
00161       *out << "after: norm(Thyra vec) = " << norm2(thyraVec) << endl;
00162     }
00163 }
00164 
00165 
00166 void EpetraOperatorWrapper::
00167 copyThyraIntoEpetra(const ConstVector<double>& thyraVec,
00168                     Epetra_MultiVector& v) const 
00169 {
00170   int numVecs = v.NumVectors();
00171   TEST_FOR_EXCEPTION(numVecs != 1, runtime_error,
00172                      "copyThyraIntoEpetra does not work with MV dimension != 1");
00173 
00174   /* get a writable pointer to the contents of the Epetra vector */
00175   double* epetraData = v[0];
00176 
00177   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00178     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00179   bool verbose = false;
00180   if (verbose)
00181     {
00182       *out << "in ETO::copyThyraIntoEpetra()" << endl;
00183       *out << "before: norm(Thyra vec) = " << norm2(thyraVec) << endl;
00184     }
00185 
00186   int offset = 0;
00187   for (int b=0; b<thyraVec.numBlocks(); b++)
00188     {
00189       ConstVector<double> block = thyraVec.getBlock(b);
00190       TEST_FOR_EXCEPTION(!isSPMD(Thyra::space(block)), runtime_error, 
00191                          "block is not SPMD");
00192       int localOffset = lowestLocallyOwnedIndex(Thyra::space(block));
00193       int localNumElems = numLocalElements(Thyra::space(block));
00194       RefCountPtr<const SpmdVectorBase<double> > spmd 
00195         = rcp_dynamic_cast<const SpmdVectorBase<double> >(block.constPtr());
00196 
00198       {
00199         /* limit the scope so that it gets destroyed, and thus committed,
00200          * before using the vector */
00201         Teuchos::RefCountPtr<const ConstDetachedVectorView<double> > view 
00202           = rcp(new ConstDetachedVectorView<double>(block.constPtr(), 
00203                                                Thyra::Range1D(localOffset,
00204                                                               localOffset+localNumElems-1), 
00205                                                true));
00206         const double* thyraData = view->values();
00207         for (int i=0; i<localNumElems; i++)
00208           {
00209             epetraData[i+offset] = thyraData[i];
00210           }
00211       }
00212       offset += localNumElems;
00213     }
00214   
00215   if (verbose)
00216     {
00217       double normOut=0;
00218       v.Norm2(&normOut);
00219       *out << "after: norm(Epetra vec) = " << normOut << endl;
00220     }
00221 }
00222 
00223 
00224 
00225 
00226 int EpetraOperatorWrapper::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00227 {
00228   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00229     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00230   bool verbose = false;
00231   if (verbose)
00232     {
00233       *out << "in ETO::Apply()" << endl;
00234       double normX=0;
00235       double normY=0;
00236       X.Norm2(&normX);
00237       Y.Norm2(&normY);
00238       *out << "before: norm(in) = " << normX << endl;
00239       *out << "before: norm(out) = " << normY << endl;
00240     }
00241 
00242   if (!useTranspose_)
00243     {
00244       Vector<double> tX = domain_.createMember();
00245       copyEpetraIntoThyra(X, tX);
00246       Vector<double> tY = range_.createMember();
00247       zeroOut(tY);
00248       thyraOp_.apply(tX, tY);
00249       copyThyraIntoEpetra(tY, Y);
00250     }
00251   else
00252     {
00253       Vector<double> tX = range_.createMember();
00254       copyEpetraIntoThyra(X, tX);
00255       Vector<double> tY = domain_.createMember();
00256       zeroOut(tY);
00257       thyraOp_.applyTranspose(tX, tY);
00258       copyThyraIntoEpetra(tY, Y);
00259     }
00260 
00261   if (verbose)
00262     {
00263       double normX=0;
00264       double normY=0;
00265       X.Norm2(&normX);
00266       Y.Norm2(&normY);
00267       *out << "after: norm(in) = " << normX << endl;
00268       *out << "after: norm(out) = " << normY << endl;
00269       *out << "leaving ETO::Apply()" << endl;
00270     }
00271   return 0;
00272 }
00273 
00274 int EpetraOperatorWrapper::ApplyInverse(const Epetra_MultiVector& X, 
00275                                       Epetra_MultiVector& Y) const
00276 {
00277   TEST_FOR_EXCEPTION(true, runtime_error,
00278                      "EpetraOperatorWrapper::ApplyInverse not implemented");
00279   return 1;
00280 }
00281 
00282 RefCountPtr<Epetra_Comm> 
00283 EpetraOperatorWrapper::getEpetraComm(const ConstLinearOperator<double>& thyraOp) const
00284 {
00285   RefCountPtr<Epetra_Comm> rtn;
00286   VectorSpace<double> vs = thyraOp.domain().getBlock(0);
00287 
00288   RefCountPtr<const SpmdVectorSpaceBase<double> > spmd 
00289     = rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(vs.constPtr());
00290 
00291   TEST_FOR_EXCEPTION(!isSPMD(vs), runtime_error, 
00292                      "EpetraOperatorWrapper requires vector space "
00293                      "blocks to be SPMD vector spaces");
00294 
00295 
00296   const SerialComm<int>* serialComm 
00297     = dynamic_cast<const SerialComm<int>*>(spmd->getComm().get());
00298 
00299 #ifdef HAVE_MPI
00300   const MpiComm<int>* mpiComm 
00301     = dynamic_cast<const MpiComm<int>*>(spmd->getComm().get());
00302 
00303   TEST_FOR_EXCEPTION(mpiComm==0 && serialComm==0, runtime_error, 
00304                      "SPMD vector space has a communicator that is "
00305                      "neither a serial comm nor an MPI comm");
00306 
00307   if (mpiComm != 0)
00308     {
00309       rtn = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
00310     }
00311   else
00312     {
00313       rtn = rcp(new Epetra_SerialComm());
00314     }
00315 #else
00316   TEST_FOR_EXCEPTION(serialComm==0, runtime_error, 
00317                      "SPMD vector space has a communicator that is "
00318                      "neither a serial comm nor an MPI comm");
00319   rtn = rcp(new Epetra_SerialComm());
00320   
00321 #endif
00322 
00323   TEST_FOR_EXCEPTION(rtn.get()==0, runtime_error, "null communicator created");
00324   return rtn;
00325 }
00326 
00327 
00328 namespace Thyra
00329 {
00330   RefCountPtr<const LinearOpBase<double> > 
00331   makeEpetraWrapper(const ConstLinearOperator<double>& thyraOp)
00332   {
00333     RefCountPtr<const LinearOpBase<double> > rtn 
00334       = rcp(new EpetraLinearOp(rcp(new EpetraOperatorWrapper(thyraOp))));
00335     return rtn;
00336   }
00337 
00338 }

Generated on Thu Sep 18 12:37:48 2008 for Epetra to Thyra Adapter Software by doxygen 1.3.9.1