Anasazi Version of the Day
AnasaziThyraAdapter.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: 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 ANASAZI_THYRA_ADAPTER_HPP
00035 #define ANASAZI_THYRA_ADAPTER_HPP
00036 
00037 #include "AnasaziMultiVecTraits.hpp"
00038 #include "AnasaziOperatorTraits.hpp"
00039 #include "AnasaziConfigDefs.hpp"
00040 
00041 #include <Thyra_DetachedMultiVectorView.hpp>
00042 #include <Thyra_MultiVectorBase.hpp>
00043 #include <Thyra_MultiVectorStdOps.hpp>
00044 
00045 namespace Anasazi {
00046   
00048   //
00049   // Implementation of the Anasazi::MultiVecTraits for Thyra::MultiVectorBase
00050   //
00052 
00061   template<class ScalarType>
00062   class MultiVecTraits< ScalarType, Thyra::MultiVectorBase<ScalarType> >
00063   {
00064   private:
00065     typedef Thyra::MultiVectorBase<ScalarType> TMVB;
00066   public:
00067 
00070 
00075     static Teuchos::RCP< Thyra::MultiVectorBase<ScalarType> > Clone( const  Thyra::MultiVectorBase<ScalarType> & mv, const int numvecs )
00076     { 
00077       Teuchos::RCP< Thyra::MultiVectorBase<ScalarType> > c = Thyra::createMembers( mv.range(), numvecs ); 
00078       return c;
00079     }
00080 
00085     static Teuchos::RCP<Thyra::MultiVectorBase<ScalarType> > 
00086     CloneCopy (const Thyra::MultiVectorBase<ScalarType>& mv)
00087     { 
00088       const int numvecs = mv.domain()->dim();
00089       // create the new multivector
00090       Teuchos::RCP<Thyra::MultiVectorBase<ScalarType> > cc = Thyra::createMembers (mv.range(), numvecs);
00091       // copy the data from the source multivector to the new multivector
00092       Thyra::assign (Teuchos::outArg (*cc), mv);
00093       return cc;
00094     }
00095 
00101     static Teuchos::RCP< Thyra::MultiVectorBase< ScalarType > > CloneCopy( const  Thyra::MultiVectorBase< ScalarType > & mv, const std::vector<int>& index )
00102     { 
00103       const int numvecs = index.size();
00104       // create the new multivector
00105       Teuchos::RCP<Thyra::MultiVectorBase<ScalarType> > cc = Thyra::createMembers (mv.range(), numvecs);
00106       // create a view to the relevant part of the source multivector
00107       Teuchos::RCP<const Thyra::MultiVectorBase<ScalarType> > view = mv.subView (numvecs, &(index[0]));
00108       // copy the data from the relevant view to the new multivector
00109       Thyra::assign (Teuchos::outArg (*cc), *view);
00110       return cc;
00111     }
00112 
00113     static Teuchos::RCP<TMVB>
00114     CloneCopy (const TMVB& mv, const Teuchos::Range1D& index)
00115     { 
00116       const int numVecs = index.size();      
00117       // Create the new multivector
00118       Teuchos::RCP<TMVB> cc = Thyra::createMembers (mv.range(), numVecs);
00119       // Create a view to the relevant part of the source multivector
00120       Teuchos::RCP<const TMVB> view = mv.subView (index);
00121       // Copy the data from the view to the new multivector.
00122       Thyra::assign (Teuchos::outArg (*cc), *view);
00123       return cc;
00124     }
00125 
00131     static Teuchos::RCP< Thyra::MultiVectorBase< ScalarType > > CloneViewNonConst(  Thyra::MultiVectorBase< ScalarType > & mv, const std::vector<int>& index )
00132     {
00133       int numvecs = index.size();
00134 
00135       // We do not assume that the indices are sorted, nor do we check that
00136       // index.size() > 0. This code is fail-safe, in the sense that a zero
00137       // length index vector will pass the error on the Thyra.
00138 
00139       // Thyra has two ways to create an indexed View:
00140       // * contiguous (via a range of columns)
00141       // * indexed (via a vector of column indices)
00142       // The former is significantly more efficient than the latter, in terms of
00143       // computations performed with/against the created view.
00144       // We will therefore check to see if the given indices are contiguous, and
00145       // if so, we will use the contiguous view creation method.
00146 
00147       int lb = index[0];
00148       bool contig = true;
00149       for (int i=0; i<numvecs; i++) {
00150         if (lb+i != index[i]) contig = false;
00151       }
00152 
00153       Teuchos::RCP<  Thyra::MultiVectorBase< ScalarType >  > cc;
00154       if (contig) {
00155         const Thyra::Range1D rng(lb,lb+numvecs-1);
00156         // create a contiguous view to the relevant part of the source multivector
00157         cc = mv.subView(rng);
00158       }
00159       else {
00160         // create an indexed view to the relevant part of the source multivector
00161         cc = mv.subView( numvecs, &(index[0]) );
00162       }
00163       return cc;
00164     }
00165 
00166     static Teuchos::RCP<TMVB> 
00167     CloneViewNonConst (TMVB& mv, const Teuchos::Range1D& index)
00168     {
00169       // We let Thyra be responsible for checking that the index range
00170       // is nonempty.
00171       //
00172       // Create and return a contiguous view to the relevant part of
00173       // the source multivector.
00174       return mv.subView (index);
00175     }
00176 
00182     static Teuchos::RCP<const  Thyra::MultiVectorBase< ScalarType > > CloneView( const  Thyra::MultiVectorBase< ScalarType > & mv, const std::vector<int>& index )
00183     {
00184       int numvecs = index.size();
00185 
00186       // We do not assume that the indices are sorted, nor do we check that
00187       // index.size() > 0. This code is fail-safe, in the sense that a zero
00188       // length index vector will pass the error on the Thyra.
00189 
00190       // Thyra has two ways to create an indexed View:
00191       // * contiguous (via a range of columns)
00192       // * indexed (via a vector of column indices)
00193       // The former is significantly more efficient than the latter, in terms of
00194       // computations performed with/against the created view.
00195       // We will therefore check to see if the given indices are contiguous, and
00196       // if so, we will use the contiguous view creation method.
00197 
00198       int lb = index[0];
00199       bool contig = true;
00200       for (int i=0; i<numvecs; i++) {
00201         if (lb+i != index[i]) contig = false;
00202       }
00203 
00204       Teuchos::RCP< const  Thyra::MultiVectorBase< ScalarType >  > cc;
00205       if (contig) {
00206         const Thyra::Range1D rng(lb,lb+numvecs-1);
00207         // create a contiguous view to the relevant part of the source multivector
00208         cc = mv.subView(rng);
00209       }
00210       else {
00211         // create an indexed view to the relevant part of the source multivector
00212         cc = mv.subView( numvecs, &(index[0]) );
00213       }
00214       return cc;
00215     }
00216 
00217     static Teuchos::RCP<const TMVB> 
00218     CloneView (const TMVB& mv, const Teuchos::Range1D& index)
00219     {
00220       // We let Thyra be responsible for checking that the index range
00221       // is nonempty.
00222       //
00223       // Create and return a contiguous view to the relevant part of
00224       // the source multivector.
00225       return mv.subView (index);
00226     }
00227 
00229 
00232 
00234     static int GetVecLength( const  Thyra::MultiVectorBase< ScalarType > & mv )
00235     { return mv.range()->dim(); }
00236 
00238     static int GetNumberVecs( const  Thyra::MultiVectorBase< ScalarType > & mv )
00239     { return mv.domain()->dim(); }
00240 
00242 
00245 
00248     static void MvTimesMatAddMv( const ScalarType alpha, const  Thyra::MultiVectorBase< ScalarType > & A, 
00249          const Teuchos::SerialDenseMatrix<int,ScalarType>& B, 
00250          const ScalarType beta,  Thyra::MultiVectorBase< ScalarType > & mv )
00251     {
00252       int m = B.numRows();
00253       int n = B.numCols();
00254       // Create a view of the B object!
00255       Teuchos::RCP< const  Thyra::MultiVectorBase< ScalarType >  >
00256         B_thyra = Thyra::createMembersView(
00257           A.domain()
00258           ,RTOpPack::ConstSubMultiVectorView<ScalarType>(0,m,0,n,&B(0,0),B.stride())
00259           );
00260       // perform the operation via A: mv <- alpha*A*B_thyra + beta*mv
00261       A.apply(Thyra::NONCONJ_ELE,*B_thyra,&mv,alpha,beta);
00262     }
00263 
00266     static void MvAddMv( const ScalarType alpha, const  Thyra::MultiVectorBase< ScalarType > & A, 
00267                          const ScalarType beta,  const  Thyra::MultiVectorBase< ScalarType > & B,  Thyra::MultiVectorBase< ScalarType > & mv )
00268     { 
00269       ScalarType coef[2], zero = Teuchos::ScalarTraits<ScalarType>::zero();
00270       const  Thyra::MultiVectorBase< ScalarType >  * in[2];
00271 
00272       in[0] = &A;
00273       in[1] = &B;
00274       coef[0] = alpha;
00275       coef[1] = beta;
00276 
00277       Thyra::linear_combination(2,coef,in,zero,&mv);
00278     }
00279 
00282     static void MvTransMv( const ScalarType alpha, const  Thyra::MultiVectorBase< ScalarType > & A, const  Thyra::MultiVectorBase< ScalarType > & mv, 
00283                            Teuchos::SerialDenseMatrix<int,ScalarType>& B )
00284     { 
00285       // Create a multivector to hold the result (m by n)
00286       int m = A.domain()->dim();
00287       int n = mv.domain()->dim();
00288       // Create a view of the B object!
00289       Teuchos::RCP<  Thyra::MultiVectorBase< ScalarType >  >
00290         B_thyra = Thyra::createMembersView(
00291           A.domain()
00292           ,RTOpPack::SubMultiVectorView<ScalarType>(0,m,0,n,&B(0,0),B.stride())
00293           );
00294       A.applyTranspose(Thyra::CONJ_ELE,mv,&*B_thyra,alpha,Teuchos::ScalarTraits<ScalarType>::zero());
00295     }
00296 
00300     static void MvDot( const  Thyra::MultiVectorBase< ScalarType > & mv, const  Thyra::MultiVectorBase< ScalarType > & A, std::vector<ScalarType> &b )
00301     { Thyra::dots(mv,A,&(b[0])); }
00302 
00305     static void 
00306     MvScale (Thyra::MultiVectorBase<ScalarType>& mv, 
00307        const ScalarType alpha)
00308     { 
00309       Thyra::scale (alpha, Teuchos::inOutArg (mv)); 
00310     }
00311     
00314     static void 
00315     MvScale (Thyra::MultiVectorBase<ScalarType>& mv, 
00316        const std::vector<ScalarType>& alpha) 
00317     { 
00318       for (unsigned int i=0; i<alpha.size(); i++) {
00319         Thyra::scale (alpha[i], mv.col(i).ptr());
00320       }
00321     }
00322     
00324 
00327 
00331     static void MvNorm( const  Thyra::MultiVectorBase< ScalarType > & mv, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec )
00332     { Thyra::norms_2(mv,&(normvec[0])); }
00333 
00335 
00338 
00341     static void SetBlock( const  Thyra::MultiVectorBase< ScalarType > & A, const std::vector<int>& index,  Thyra::MultiVectorBase< ScalarType > & mv )
00342     { 
00343       // Extract the "numvecs" columns of mv indicated by the index vector.
00344       int numvecs = index.size();
00345       std::vector<int> indexA(numvecs);
00346       int numAcols = A.domain()->dim();
00347       for (int i=0; i<numvecs; i++) {
00348         indexA[i] = i;
00349       }
00350       // Thyra::assign requires that both arguments have the same number of 
00351       // vectors. Enforce this, by shrinking one to match the other.
00352       if ( numAcols < numvecs ) {
00353         // A does not have enough columns to satisfy index_plus. Shrink
00354         // index_plus.
00355         numvecs = numAcols;
00356       }
00357       else if ( numAcols > numvecs ) {
00358         numAcols = numvecs;
00359         indexA.resize( numAcols );
00360       }
00361       // create a view to the relevant part of the source multivector
00362       Teuchos::RCP< const  Thyra::MultiVectorBase< ScalarType >  > relsource = A.subView( numAcols, &(indexA[0]) );
00363       // create a view to the relevant part of the destination multivector
00364       Teuchos::RCP<  Thyra::MultiVectorBase< ScalarType >  > reldest = mv.subView( numvecs, &(index[0]) );
00365       // copy the data to the destination multivector subview
00366       Thyra::assign (Teuchos::outArg (*reldest), *relsource);
00367     }
00368 
00369     static void 
00370     SetBlock (const TMVB& A, const Teuchos::Range1D& index, TMVB& mv)
00371     { 
00372       const int numColsA = A.domain()->dim();
00373       const int numColsMv = mv.domain()->dim();
00374       // 'index' indexes into mv; it's the index set of the target.
00375       const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
00376       // We can't take more columns out of A than A has.
00377       const bool validSource = index.size() <= numColsA;
00378 
00379       if (! validIndex || ! validSource)
00380   {
00381     std::ostringstream os;
00382     os << "Anasazi::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar> "
00383       ">::SetBlock(A, [" << index.lbound() << ", " << index.ubound() 
00384        << "], mv): ";
00385     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00386            os.str() << "Range lower bound must be nonnegative.");
00387     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
00388            os.str() << "Range upper bound must be less than "
00389            "the number of columns " << numColsA << " in the "
00390            "'mv' output argument.");
00391     TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
00392            os.str() << "Range must have no more elements than"
00393            " the number of columns " << numColsA << " in the "
00394            "'A' input argument.");
00395     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
00396   }
00397 
00398       // View of the relevant column(s) of the target multivector mv.
00399       // We avoid view creation overhead by only creating a view if
00400       // the index range is different than [0, (# columns in mv) - 1].
00401       Teuchos::RCP<TMVB> mv_view;
00402       if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
00403   mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
00404       else
00405   mv_view = mv.subView (index);
00406 
00407       // View of the relevant column(s) of the source multivector A.
00408       // If A has fewer columns than mv_view, then create a view of
00409       // the first index.size() columns of A.
00410       Teuchos::RCP<const TMVB> A_view;
00411       if (index.size() == numColsA)
00412   A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
00413       else
00414   A_view = A.subView (Teuchos::Range1D(0, index.size()-1));
00415 
00416       // Copy the data to the destination multivector.
00417       Thyra::assign (Teuchos::outArg (*mv_view), *A_view);
00418     }
00419 
00420     static void 
00421     Assign (const TMVB& A, TMVB& mv)
00422     { 
00423       using Teuchos::ptr;
00424       using Teuchos::Range1D;
00425       using Teuchos::RCP;
00426 
00427       const int numColsA = A.domain()->dim();
00428       const int numColsMv = mv.domain()->dim();
00429       if (numColsA > numColsMv) {
00430   const std::string prefix ("Anasazi::MultiVecTraits<Scalar, "
00431           "Thyra::MultiVectorBase<Scalar>"
00432           " >::Assign(A, mv): ");
00433   TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
00434            prefix << "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 (outArg (mv), A);
00442       } 
00443       else {
00444   RCP<TMVB> mv_view = CloneViewNonConst (mv, Range1D(0, numColsA-1));
00445   Thyra::assign (outArg (*mv_view), A);
00446       }
00447     }
00448 
00451     static void MvRandom(  Thyra::MultiVectorBase< ScalarType > & 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(-Teuchos::ScalarTraits<ScalarType>::one(),
00456                         Teuchos::ScalarTraits<ScalarType>::one(),
00457                        &mv);
00458     }
00459 
00462     static void 
00463     MvInit (Thyra::MultiVectorBase< ScalarType >& mv, 
00464       ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero())
00465     { 
00466       Thyra::assign (Teuchos::outArg (mv), alpha); 
00467     }
00468 
00470 
00473 
00476     static void MvPrint( const  Thyra::MultiVectorBase< ScalarType > & mv, std::ostream& os )
00477     { 
00478        Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&os,false));    
00479        out->setf(std::ios_base::scientific);    
00480        out->precision(16);   
00481        mv.describe(*out,Teuchos::VERB_EXTREME);
00482     }
00483 
00485 
00486   };        
00487 
00489   //
00490   // Implementation of the Anasazi::OperatorTraits for Thyra::LinearOpBase
00491   //
00493 
00503   template <class ScalarType> 
00504   class OperatorTraits < ScalarType, Thyra::MultiVectorBase<ScalarType>, Thyra::LinearOpBase<ScalarType> >
00505   {
00506   public:
00507     
00511     static void Apply ( const Thyra::LinearOpBase< ScalarType >& Op, const  Thyra::MultiVectorBase< ScalarType > & x,  Thyra::MultiVectorBase< ScalarType > & y )
00512     { 
00513       Op.apply(Thyra::NONCONJ_ELE,x,&y);
00514     }
00515     
00516   };
00517   
00518 } // end of Anasazi namespace 
00519 
00520 #endif 
00521 // end of file ANASAZI_THYRA_ADAPTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends