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 "BelosMultiVecTraits.hpp"
00038 #include "BelosOperatorTraits.hpp"
00039 #include "BelosConfigDefs.hpp"
00040 
00041 #include <Thyra_DetachedMultiVectorView.hpp>
00042 #include <Thyra_MultiVectorBase.hpp>
00043 #include <Thyra_MultiVectorStdOps.hpp>
00044 
00045 namespace Belos {
00046   
00048   //
00049   // Implementation of the Belos::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     typedef Teuchos::ScalarTraits<ScalarType> ST;
00067     typedef typename ST::magnitudeType magType;
00068   public:
00069 
00072 
00077     static Teuchos::RCP<TMVB> Clone( const TMVB& mv, const int numvecs )
00078     { 
00079       Teuchos::RCP<TMVB> c = Thyra::createMembers( mv.range(), numvecs ); 
00080       return c;
00081     }
00082 
00087     static Teuchos::RCP<TMVB> CloneCopy( const TMVB& mv )
00088     { 
00089       int numvecs = mv.domain()->dim();
00090       // create the new multivector
00091       Teuchos::RCP< TMVB > cc = Thyra::createMembers( mv.range(), numvecs );
00092       // copy the data from the source multivector to the new multivector
00093       Thyra::assign(&*cc, mv);
00094       return cc;
00095     }
00096 
00102     static Teuchos::RCP<TMVB> CloneCopy( const TMVB& mv, const std::vector<int>& index )
00103     { 
00104       int numvecs = index.size();
00105       // create the new multivector
00106       Teuchos::RCP< TMVB > cc = Thyra::createMembers( mv.range(), numvecs );
00107       // create a view to the relevant part of the source multivector
00108       Teuchos::RCP< const TMVB > view = mv.subView( numvecs, &(index[0]) );
00109       // copy the data from the relevant view to the new multivector
00110       Thyra::assign(&*cc, *view);
00111       return cc;
00112     }
00113 
00119     static Teuchos::RCP<TMVB> CloneViewNonConst( TMVB& mv, const std::vector<int>& index )
00120     {
00121       int numvecs = index.size();
00122 
00123       // We do not assume that the indices are sorted, nor do we check that
00124       // index.size() > 0. This code is fail-safe, in the sense that a zero
00125       // length index std::vector will pass the error on the Thyra.
00126 
00127       // Thyra has two ways to create an indexed View:
00128       // * contiguous (via a range of columns)
00129       // * indexed (via a std::vector of column indices)
00130       // The former is significantly more efficient than the latter, in terms of
00131       // computations performed with/against the created view.
00132       // We will therefore check to see if the given indices are contiguous, and
00133       // if so, we will use the contiguous view creation method.
00134 
00135       int lb = index[0];
00136       bool contig = true;
00137       for (int i=0; i<numvecs; i++) {
00138         if (lb+i != index[i]) contig = false;
00139       }
00140 
00141       Teuchos::RCP< TMVB > cc;
00142       if (contig) {
00143         const Thyra::Range1D rng(lb,lb+numvecs-1);
00144         // create a contiguous view to the relevant part of the source multivector
00145         cc = mv.subView(rng);
00146       }
00147       else {
00148         // create an indexed view to the relevant part of the source multivector
00149         cc = mv.subView( numvecs, &(index[0]) );
00150       }
00151       return cc;
00152     }
00153 
00159     static Teuchos::RCP<const TMVB> CloneView( const TMVB& mv, const std::vector<int>& index )
00160     {
00161       int numvecs = index.size();
00162 
00163       // We do not assume that the indices are sorted, nor do we check that
00164       // index.size() > 0. This code is fail-safe, in the sense that a zero
00165       // length index std::vector will pass the error on the Thyra.
00166 
00167       // Thyra has two ways to create an indexed View:
00168       // * contiguous (via a range of columns)
00169       // * indexed (via a std::vector of column indices)
00170       // The former is significantly more efficient than the latter, in terms of
00171       // computations performed with/against the created view.
00172       // We will therefore check to see if the given indices are contiguous, and
00173       // if so, we will use the contiguous view creation method.
00174 
00175       int lb = index[0];
00176       bool contig = true;
00177       for (int i=0; i<numvecs; i++) {
00178         if (lb+i != index[i]) contig = false;
00179       }
00180 
00181       Teuchos::RCP< const TMVB > cc;
00182       if (contig) {
00183         const Thyra::Range1D rng(lb,lb+numvecs-1);
00184         // create a contiguous view to the relevant part of the source multivector
00185         cc = mv.subView(rng);
00186       }
00187       else {
00188         // create an indexed view to the relevant part of the source multivector
00189         cc = mv.subView( numvecs, &(index[0]) );
00190       }
00191       return cc;
00192     }
00193 
00195 
00198 
00200     static int GetVecLength( const TMVB& mv )
00201     { return mv.range()->dim(); }
00202 
00204     static int GetNumberVecs( const TMVB& mv )
00205     { return mv.domain()->dim(); }
00206 
00208 
00211 
00214     static void MvTimesMatAddMv( const ScalarType alpha, const TMVB& A, 
00215          const Teuchos::SerialDenseMatrix<int,ScalarType>& B, 
00216          const ScalarType beta, TMVB& mv )
00217     {
00218       using Teuchos::arrayView; using Teuchos::arcpFromArrayView;
00219       const int m = B.numRows();
00220       const int n = B.numCols();
00221       // Create a view of the B object!
00222       Teuchos::RCP< const TMVB >
00223         B_thyra = Thyra::createMembersView(
00224           A.domain(),
00225           RTOpPack::ConstSubMultiVectorView<ScalarType>(
00226             0, m, 0, n,
00227             arcpFromArrayView(arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
00228             )
00229           );
00230       // perform the operation via A: mv <- alpha*A*B_thyra + beta*mv
00231       Thyra::apply<ScalarType>(A, Thyra::NOTRANS, *B_thyra, Teuchos::outArg(mv), alpha, beta);
00232     }
00233 
00236     static void MvAddMv( const ScalarType alpha, const TMVB& A, 
00237                          const ScalarType beta,  const TMVB& B, TMVB& mv )
00238     { 
00239       ScalarType coef[2], zero = Teuchos::ScalarTraits<ScalarType>::zero();
00240       const TMVB * in[2];
00241 
00242       in[0] = &A;
00243       in[1] = &B;
00244       coef[0] = alpha;
00245       coef[1] = beta;
00246 
00247       Thyra::linear_combination(2,coef,in,zero,&mv);
00248     }
00249 
00252     static void MvScale ( TMVB& mv, const ScalarType alpha )
00253     { Thyra::scale(alpha,&mv); }
00254 
00257     static void MvScale ( TMVB& mv, const std::vector<ScalarType>& alpha )
00258     {
00259       for (unsigned int i=0; i<alpha.size(); i++) {
00260         Thyra::scale<ScalarType>(alpha[i], mv.col(i).ptr());
00261       }
00262     }
00263 
00266     static void MvTransMv( const ScalarType alpha, const TMVB& A, const TMVB& mv, 
00267                            Teuchos::SerialDenseMatrix<int,ScalarType>& B )
00268     { 
00269       using Teuchos::arrayView; using Teuchos::arcpFromArrayView;
00270       // Create a multivector to hold the result (m by n)
00271       int m = A.domain()->dim();
00272       int n = mv.domain()->dim();
00273       // Create a view of the B object!
00274       Teuchos::RCP< TMVB >
00275         B_thyra = Thyra::createMembersView(
00276           A.domain(),
00277           RTOpPack::SubMultiVectorView<ScalarType>(
00278             0, m, 0, n,
00279             arcpFromArrayView(arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
00280             )
00281           );
00282       Thyra::apply<ScalarType>(A, Thyra::CONJTRANS, mv, B_thyra.ptr(), alpha);
00283     }
00284 
00288     static void MvDot( const TMVB& mv, const TMVB& A, std::vector<ScalarType>& b )
00289     { Thyra::dots(mv,A,&(b[0])); }
00290 
00292 
00295 
00299     static void MvNorm( const TMVB& mv, std::vector<magType>& normvec, NormType type = TwoNorm )
00300     { Thyra::norms_2(mv,&(normvec[0])); }
00301 
00303 
00306 
00309     static void SetBlock( const TMVB& A, const std::vector<int>& index, TMVB& mv )
00310     { 
00311       // Extract the "numvecs" columns of mv indicated by the index std::vector.
00312       int numvecs = index.size();
00313       std::vector<int> indexA(numvecs);
00314       int numAcols = A.domain()->dim();
00315       for (int i=0; i<numvecs; i++) {
00316         indexA[i] = i;
00317       }
00318       // Thyra::assign requires that both arguments have the same number of 
00319       // vectors. Enforce this, by shrinking one to match the other.
00320       if ( numAcols < numvecs ) {
00321         // A does not have enough columns to satisfy index_plus. Shrink
00322         // index_plus.
00323         numvecs = numAcols;
00324       }
00325       else if ( numAcols > numvecs ) {
00326         numAcols = numvecs;
00327         indexA.resize( numAcols );
00328       }
00329       // create a view to the relevant part of the source multivector
00330       Teuchos::RCP< const TMVB > relsource = A.subView( numAcols, &(indexA[0]) );
00331       // create a view to the relevant part of the destination multivector
00332       Teuchos::RCP< TMVB > reldest = mv.subView( numvecs, &(index[0]) );
00333       // copy the data to the destination multivector subview
00334       Thyra::assign(&*reldest, *relsource);
00335     }
00336 
00339     static void MvRandom( TMVB& mv )
00340     { 
00341       // Thyra::randomize generates via a uniform distribution on [l,u]
00342       // We will use this to generate on [-1,1]
00343       Thyra::randomize(-Teuchos::ScalarTraits<ScalarType>::one(),
00344                         Teuchos::ScalarTraits<ScalarType>::one(),
00345                        &mv);
00346     }
00347 
00350     static void MvInit( TMVB& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
00351     { Thyra::assign(&mv,alpha); }
00352 
00354 
00357 
00360     static void MvPrint( const TMVB& mv, std::ostream& os )
00361       { os << describe(mv,Teuchos::VERB_EXTREME); }
00362 
00364 
00365   };        
00366 
00368   //
00369   // Implementation of the Belos::OperatorTraits for Thyra::LinearOpBase
00370   //
00372 
00381   template <class ScalarType> 
00382   class OperatorTraits < ScalarType, Thyra::MultiVectorBase<ScalarType>, Thyra::LinearOpBase<ScalarType> >
00383   {
00384   private:
00385     typedef Thyra::MultiVectorBase<ScalarType> TMVB;
00386     typedef Thyra::LinearOpBase<ScalarType>    TLOB;
00387   public:
00388     
00392     static void Apply ( const TLOB& Op, const TMVB& x, TMVB& y )
00393     { 
00394       Thyra::apply<ScalarType>(Op, Thyra::NOTRANS, x, Teuchos::outArg(y));
00395     }
00396     
00397   };
00398   
00399 } // end of Belos namespace 
00400 
00401 #endif 
00402 // end of file BELOS_THYRA_ADAPTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 10:19:51 2011 for Stratimikos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3