Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
BelosThyraAdapter.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Belos: Block Eigensolvers 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 
00034 #ifndef BELOS_THYRA_ADAPTER_HPP
00035 #define BELOS_THYRA_ADAPTER_HPP
00036 
00037 #include "BelosConfigDefs.hpp"
00038 #include "BelosMultiVecTraits.hpp"
00039 #include "BelosOperatorTraits.hpp"
00040 
00041 #include <Thyra_DetachedMultiVectorView.hpp>
00042 #include <Thyra_MultiVectorBase.hpp>
00043 #include <Thyra_MultiVectorStdOps.hpp>
00044 #ifdef HAVE_BELOS_TSQR
00045 #  include <Thyra_TsqrAdaptor.hpp>
00046 #endif // HAVE_BELOS_TSQR
00047 
00048 namespace Belos {
00049   
00051   //
00052   // Implementation of the Belos::MultiVecTraits for Thyra::MultiVectorBase
00053   //
00055 
00062   template<class ScalarType>
00063   class MultiVecTraits< ScalarType, Thyra::MultiVectorBase<ScalarType> >
00064   {
00065   private:
00066     typedef Thyra::MultiVectorBase<ScalarType> TMVB;
00067     typedef Teuchos::ScalarTraits<ScalarType> ST;
00068     typedef typename ST::magnitudeType magType;
00069 
00070   public:
00071 
00074 
00079     static Teuchos::RCP<TMVB> Clone( const TMVB& mv, const int numvecs )
00080     { 
00081       Teuchos::RCP<TMVB> c = Thyra::createMembers( mv.range(), numvecs ); 
00082       return c;
00083     }
00084 
00089     static Teuchos::RCP<TMVB> CloneCopy( const TMVB& mv )
00090     { 
00091       int numvecs = mv.domain()->dim();
00092       // create the new multivector
00093       Teuchos::RCP< TMVB > cc = Thyra::createMembers( mv.range(), numvecs );
00094       // copy the data from the source multivector to the new multivector
00095       Thyra::assign(cc.ptr(), mv);
00096       return cc;
00097     }
00098 
00104     static Teuchos::RCP<TMVB> CloneCopy( const TMVB& mv, const std::vector<int>& index )
00105     { 
00106       int numvecs = index.size();
00107       // create the new multivector
00108       Teuchos::RCP<TMVB> cc = Thyra::createMembers( mv.range(), numvecs );
00109       // create a view to the relevant part of the source multivector
00110       Teuchos::RCP<const TMVB> view = mv.subView(index);
00111       // copy the data from the relevant view to the new multivector
00112       Thyra::assign(cc.ptr(), *view);
00113       return cc;
00114     }
00115 
00116     static Teuchos::RCP<TMVB> 
00117     CloneCopy (const TMVB& mv, const Teuchos::Range1D& index)
00118     { 
00119       const int numVecs = index.size();      
00120       // Create the new multivector
00121       Teuchos::RCP<TMVB> cc = Thyra::createMembers (mv.range(), numVecs);
00122       // Create a view to the relevant part of the source multivector
00123       Teuchos::RCP<const TMVB> view = mv.subView (index);
00124       // Copy the data from the view to the new multivector.
00125       Thyra::assign (cc.ptr(), *view);
00126       return cc;
00127     }
00128 
00134     static Teuchos::RCP<TMVB> CloneViewNonConst( TMVB& mv, const std::vector<int>& index )
00135     {
00136       int numvecs = index.size();
00137 
00138       // We do not assume that the indices are sorted, nor do we check that
00139       // index.size() > 0. This code is fail-safe, in the sense that a zero
00140       // length index std::vector will pass the error on the Thyra.
00141 
00142       // Thyra has two ways to create an indexed View:
00143       // * contiguous (via a range of columns)
00144       // * indexed (via a std::vector of column indices)
00145       // The former is significantly more efficient than the latter, in terms of
00146       // computations performed with/against the created view.
00147       // We will therefore check to see if the given indices are contiguous, and
00148       // if so, we will use the contiguous view creation method.
00149 
00150       int lb = index[0];
00151       bool contig = true;
00152       for (int i=0; i<numvecs; i++) {
00153         if (lb+i != index[i]) contig = false;
00154       }
00155 
00156       Teuchos::RCP< TMVB > cc;
00157       if (contig) {
00158         const Thyra::Range1D rng(lb,lb+numvecs-1);
00159         // create a contiguous view to the relevant part of the source multivector
00160         cc = mv.subView(rng);
00161       }
00162       else {
00163         // create an indexed view to the relevant part of the source multivector
00164         cc = mv.subView(index);
00165       }
00166       return cc;
00167     }
00168 
00169     static Teuchos::RCP<TMVB> 
00170     CloneViewNonConst (TMVB& mv, const Teuchos::Range1D& index)
00171     {
00172       // We let Thyra be responsible for checking that the index range
00173       // is nonempty.
00174       //
00175       // Create and return a contiguous view to the relevant part of
00176       // the source multivector.
00177       return mv.subView (index);
00178     }
00179 
00180 
00186     static Teuchos::RCP<const TMVB> CloneView( const TMVB& mv, const std::vector<int>& index )
00187     {
00188       int numvecs = index.size();
00189 
00190       // We do not assume that the indices are sorted, nor do we check that
00191       // index.size() > 0. This code is fail-safe, in the sense that a zero
00192       // length index std::vector will pass the error on the Thyra.
00193 
00194       // Thyra has two ways to create an indexed View:
00195       // * contiguous (via a range of columns)
00196       // * indexed (via a std::vector of column indices)
00197       // The former is significantly more efficient than the latter, in terms of
00198       // computations performed with/against the created view.
00199       // We will therefore check to see if the given indices are contiguous, and
00200       // if so, we will use the contiguous view creation method.
00201 
00202       int lb = index[0];
00203       bool contig = true;
00204       for (int i=0; i<numvecs; i++) {
00205         if (lb+i != index[i]) contig = false;
00206       }
00207 
00208       Teuchos::RCP< const TMVB > cc;
00209       if (contig) {
00210         const Thyra::Range1D rng(lb,lb+numvecs-1);
00211         // create a contiguous view to the relevant part of the source multivector
00212         cc = mv.subView(rng);
00213       }
00214       else {
00215         // create an indexed view to the relevant part of the source multivector
00216         cc = mv.subView(index);
00217       }
00218       return cc;
00219     }
00220 
00221     static Teuchos::RCP<const TMVB> 
00222     CloneView (const TMVB& mv, const Teuchos::Range1D& index)
00223     {
00224       // We let Thyra be responsible for checking that the index range
00225       // is nonempty.
00226       //
00227       // Create and return a contiguous view to the relevant part of
00228       // the source multivector.
00229       return mv.subView (index);
00230     }
00231 
00233 
00236 
00238     static int GetVecLength( const TMVB& mv )
00239     { return mv.range()->dim(); }
00240 
00242     static int GetNumberVecs( const TMVB& mv )
00243     { return mv.domain()->dim(); }
00244 
00246 
00249 
00252     static void MvTimesMatAddMv( const ScalarType alpha, const TMVB& A, 
00253          const Teuchos::SerialDenseMatrix<int,ScalarType>& B, 
00254          const ScalarType beta, TMVB& mv )
00255     {
00256       using Teuchos::arrayView; using Teuchos::arcpFromArrayView;
00257       const int m = B.numRows();
00258       const int n = B.numCols();
00259       // Create a view of the B object!
00260       Teuchos::RCP< const TMVB >
00261         B_thyra = Thyra::createMembersView(
00262           A.domain(),
00263           RTOpPack::ConstSubMultiVectorView<ScalarType>(
00264             0, m, 0, n,
00265             arcpFromArrayView(arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
00266             )
00267           );
00268       // perform the operation via A: mv <- alpha*A*B_thyra + beta*mv
00269       Thyra::apply<ScalarType>(A, Thyra::NOTRANS, *B_thyra, Teuchos::outArg(mv), alpha, beta);
00270     }
00271 
00274     static void MvAddMv( const ScalarType alpha, const TMVB& A, 
00275                          const ScalarType beta,  const TMVB& B, TMVB& mv )
00276     { 
00277       typedef Teuchos::ScalarTraits<ScalarType> ST;
00278       using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::inoutArg;
00279 
00280       Thyra::linear_combination<ScalarType>(
00281         tuple(alpha, beta)(), tuple(ptrInArg(A), ptrInArg(B))(), ST::zero(), inoutArg(mv));
00282     }
00283 
00286     static void MvScale ( TMVB& mv, const ScalarType alpha )
00287       { Thyra::scale(alpha, Teuchos::inoutArg(mv)); }
00288 
00291     static void MvScale (TMVB& mv, const std::vector<ScalarType>& alpha)
00292     {
00293       for (unsigned int i=0; i<alpha.size(); i++) {
00294         Thyra::scale<ScalarType> (alpha[i], mv.col(i).ptr());
00295       }
00296     }
00297 
00300     static void MvTransMv( const ScalarType alpha, const TMVB& A, const TMVB& mv, 
00301       Teuchos::SerialDenseMatrix<int,ScalarType>& B )
00302     { 
00303       using Teuchos::arrayView; using Teuchos::arcpFromArrayView;
00304       // Create a multivector to hold the result (m by n)
00305       int m = A.domain()->dim();
00306       int n = mv.domain()->dim();
00307       // Create a view of the B object!
00308       Teuchos::RCP< TMVB >
00309         B_thyra = Thyra::createMembersView(
00310           A.domain(),
00311           RTOpPack::SubMultiVectorView<ScalarType>(
00312             0, m, 0, n,
00313             arcpFromArrayView(arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
00314             )
00315           );
00316       Thyra::apply<ScalarType>(A, Thyra::CONJTRANS, mv, B_thyra.ptr(), alpha);
00317     }
00318 
00322     static void MvDot( const TMVB& mv, const TMVB& A, std::vector<ScalarType>& b )
00323       { Thyra::dots(mv, A, Teuchos::arrayViewFromVector(b)); }
00324 
00326 
00329 
00333     static void MvNorm( const TMVB& mv, std::vector<magType>& normvec,
00334       NormType type = TwoNorm )
00335       { Thyra::norms_2(mv, Teuchos::arrayViewFromVector(normvec)); }
00336 
00338 
00341 
00344     static void SetBlock( const TMVB& A, const std::vector<int>& index, TMVB& mv )
00345     { 
00346       // Extract the "numvecs" columns of mv indicated by the index std::vector.
00347       int numvecs = index.size();
00348       std::vector<int> indexA(numvecs);
00349       int numAcols = A.domain()->dim();
00350       for (int i=0; i<numvecs; i++) {
00351         indexA[i] = i;
00352       }
00353       // Thyra::assign requires that both arguments have the same number of 
00354       // vectors. Enforce this, by shrinking one to match the other.
00355       if ( numAcols < numvecs ) {
00356         // A does not have enough columns to satisfy index_plus. Shrink
00357         // index_plus.
00358         numvecs = numAcols;
00359       }
00360       else if ( numAcols > numvecs ) {
00361         numAcols = numvecs;
00362         indexA.resize( numAcols );
00363       }
00364       // create a view to the relevant part of the source multivector
00365       Teuchos::RCP< const TMVB > relsource = A.subView(indexA);
00366       // create a view to the relevant part of the destination multivector
00367       Teuchos::RCP< TMVB > reldest = mv.subView(index);
00368       // copy the data to the destination multivector subview
00369       Thyra::assign(reldest.ptr(), *relsource);
00370     }
00371 
00372     static void 
00373     SetBlock (const TMVB& A, const Teuchos::Range1D& index, TMVB& mv)
00374     { 
00375       const int numColsA = A.domain()->dim();
00376       const int numColsMv = mv.domain()->dim();
00377       // 'index' indexes into mv; it's the index set of the target.
00378       const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
00379       // We can't take more columns out of A than A has.
00380       const bool validSource = index.size() <= numColsA;
00381 
00382       if (! validIndex || ! validSource)
00383   {
00384     std::ostringstream os;
00385     os << "Belos::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar> "
00386       ">::SetBlock(A, [" << index.lbound() << ", " << index.ubound() 
00387        << "], mv): ";
00388     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00389            os.str() << "Range lower bound must be nonnegative.");
00390     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
00391            os.str() << "Range upper bound must be less than "
00392            "the number of columns " << numColsA << " in the "
00393            "'mv' output argument.");
00394     TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
00395            os.str() << "Range must have no more elements than"
00396            " the number of columns " << numColsA << " in the "
00397            "'A' input argument.");
00398     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
00399   }
00400 
00401       // View of the relevant column(s) of the target multivector mv.
00402       // We avoid view creation overhead by only creating a view if
00403       // the index range is different than [0, (# columns in mv) - 1].
00404       Teuchos::RCP<TMVB> mv_view;
00405       if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
00406   mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
00407       else
00408   mv_view = mv.subView (index);
00409 
00410       // View of the relevant column(s) of the source multivector A.
00411       // If A has fewer columns than mv_view, then create a view of
00412       // the first index.size() columns of A.
00413       Teuchos::RCP<const TMVB> A_view;
00414       if (index.size() == numColsA)
00415   A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
00416       else
00417   A_view = A.subView (Teuchos::Range1D(0, index.size()-1));
00418 
00419       // Copy the data to the destination multivector.
00420       Thyra::assign(mv_view.ptr(), *A_view);
00421     }
00422 
00423     static void 
00424     Assign (const TMVB& A, TMVB& mv)
00425     { 
00426       const int numColsA = A.domain()->dim();
00427       const int numColsMv = mv.domain()->dim();
00428       if (numColsA > numColsMv)
00429   {
00430     std::ostringstream os;
00431     os << "Belos::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar>"
00432       " >::Assign(A, mv): ";
00433     TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
00434            os.str() << "Input multivector 'A' has " 
00435            << numColsA << " columns, but output multivector "
00436            "'mv' has only " << numColsMv << " columns.");
00437     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
00438   }
00439       // Copy the data to the destination multivector.
00440       if (numColsA == numColsMv) {
00441   Thyra::assign (Teuchos::outArg (mv), A);
00442       } else {
00443   Teuchos::RCP<TMVB> mv_view = 
00444     CloneViewNonConst (mv, Teuchos::Range1D(0, numColsA-1));
00445   Thyra::assign (mv_view.ptr(), A);
00446       }
00447     }
00448 
00451     static void MvRandom( TMVB& mv )
00452     { 
00453       // Thyra::randomize generates via a uniform distribution on [l,u]
00454       // We will use this to generate on [-1,1]
00455       Thyra::randomize<ScalarType>(
00456         -Teuchos::ScalarTraits<ScalarType>::one(),
00457         Teuchos::ScalarTraits<ScalarType>::one(),
00458         Teuchos::outArg(mv));
00459     }
00460 
00462     static void 
00463     MvInit (TMVB& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero())
00464     { 
00465       Thyra::assign (Teuchos::outArg (mv), alpha); 
00466     }
00467 
00469 
00472 
00475     static void MvPrint( const TMVB& mv, std::ostream& os )
00476       { os << describe(mv,Teuchos::VERB_EXTREME); }
00477 
00479 
00480 #ifdef HAVE_BELOS_TSQR
00481 
00482 
00483 
00484 
00485 
00486     typedef Thyra::TsqrAdaptor< ScalarType > tsqr_adaptor_type;
00487 #endif // HAVE_BELOS_TSQR
00488   };        
00489 
00491   //
00492   // Implementation of the Belos::OperatorTraits for Thyra::LinearOpBase
00493   //
00495 
00503   template<class ScalarType> 
00504   class OperatorTraits <ScalarType, 
00505       Thyra::MultiVectorBase<ScalarType>, 
00506       Thyra::LinearOpBase<ScalarType> >
00507   {
00508   private:
00509     typedef Thyra::MultiVectorBase<ScalarType> TMVB;
00510     typedef Thyra::LinearOpBase<ScalarType>    TLOB;
00511 
00512   public:
00528     static void 
00529     Apply (const TLOB& Op, 
00530      const TMVB& x, 
00531      TMVB& y,
00532      ETrans trans = NOTRANS)
00533     { 
00534       Thyra::EOpTransp whichOp;
00535 
00536       // We don't check here whether the operator implements the
00537       // requested operation.  Call HasApplyTranspose() to check.
00538       // Thyra::LinearOpBase implementations are not required to
00539       // implement NOTRANS.  However, Belos needs NOTRANS
00540       // (obviously!), so we assume that Op implements NOTRANS.
00541       if (trans == NOTRANS)
00542   whichOp = Thyra::NOTRANS;
00543       else if (trans == TRANS)
00544   whichOp = Thyra::TRANS;
00545       else if (trans == CONJTRANS)
00546   whichOp = Thyra::CONJTRANS;
00547       else
00548   TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
00549          "Belos::OperatorTraits::Apply (Thyra specialization): "
00550          "'trans' argument must be neither NOTRANS=" << NOTRANS 
00551          << ", TRANS=" << TRANS << ", or CONJTRANS=" << CONJTRANS
00552          << ", but instead has an invalid value of " << trans << ".");
00553       Thyra::apply<ScalarType>(Op, whichOp, x, Teuchos::outArg(y));
00554     }
00555 
00557     static bool HasApplyTranspose (const TLOB& Op)
00558     {
00559       typedef Teuchos::ScalarTraits<ScalarType> STS;
00560 
00561       // Thyra::LinearOpBase's interface lets you check whether the
00562       // operator implements any of all four possible combinations of
00563       // conjugation and transpose.  Belos only needs transpose
00564       // (TRANS) if the operator is real; in that case, Apply() does
00565       // the same thing with trans = CONJTRANS or TRANS.  If the
00566       // operator is complex, Belos needs both transpose and conjugate
00567       // transpose (CONJTRANS) if the operator is complex.  
00568       return Op.opSupported (Thyra::TRANS) && 
00569   (! STS::isComplex || Op.opSupported (Thyra::CONJTRANS));
00570     }
00571   };
00572 
00573 } // end of Belos namespace 
00574 
00575 #endif 
00576 // end of file BELOS_THYRA_ADAPTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines