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   public:
00065 
00068 
00073     static Teuchos::RCP< Thyra::MultiVectorBase<ScalarType> > Clone( const  Thyra::MultiVectorBase<ScalarType> & mv, const int numvecs )
00074     { 
00075       Teuchos::RCP< Thyra::MultiVectorBase<ScalarType> > c = Thyra::createMembers( mv.range(), numvecs ); 
00076       return c;
00077     }
00078 
00083     static Teuchos::RCP< Thyra::MultiVectorBase<ScalarType> > CloneCopy( const  Thyra::MultiVectorBase< ScalarType > & mv )
00084     { 
00085       int numvecs = mv.domain()->dim();
00086       // create the new multivector
00087       Teuchos::RCP<  Thyra::MultiVectorBase< ScalarType >  > cc = Thyra::createMembers( mv.range(), numvecs );
00088       // copy the data from the source multivector to the new multivector
00089       Thyra::assign(&*cc, mv);
00090       return cc;
00091     }
00092 
00098     static Teuchos::RCP< Thyra::MultiVectorBase< ScalarType > > CloneCopy( const  Thyra::MultiVectorBase< ScalarType > & mv, const std::vector<int>& index )
00099     { 
00100       int numvecs = index.size();
00101       // create the new multivector
00102       Teuchos::RCP<  Thyra::MultiVectorBase< ScalarType >  > cc = Thyra::createMembers( mv.range(), numvecs );
00103       // create a view to the relevant part of the source multivector
00104       Teuchos::RCP< const  Thyra::MultiVectorBase< ScalarType >  > view = mv.subView( numvecs, &(index[0]) );
00105       // copy the data from the relevant view to the new multivector
00106       Thyra::assign(&*cc, *view);
00107       return cc;
00108     }
00109 
00115     static Teuchos::RCP< Thyra::MultiVectorBase< ScalarType > > CloneViewNonConst(  Thyra::MultiVectorBase< ScalarType > & mv, const std::vector<int>& index )
00116     {
00117       int numvecs = index.size();
00118 
00119       // We do not assume that the indices are sorted, nor do we check that
00120       // index.size() > 0. This code is fail-safe, in the sense that a zero
00121       // length index vector will pass the error on the Thyra.
00122 
00123       // Thyra has two ways to create an indexed View:
00124       // * contiguous (via a range of columns)
00125       // * indexed (via a vector of column indices)
00126       // The former is significantly more efficient than the latter, in terms of
00127       // computations performed with/against the created view.
00128       // We will therefore check to see if the given indices are contiguous, and
00129       // if so, we will use the contiguous view creation method.
00130 
00131       int lb = index[0];
00132       bool contig = true;
00133       for (int i=0; i<numvecs; i++) {
00134         if (lb+i != index[i]) contig = false;
00135       }
00136 
00137       Teuchos::RCP<  Thyra::MultiVectorBase< ScalarType >  > cc;
00138       if (contig) {
00139         const Thyra::Range1D rng(lb,lb+numvecs-1);
00140         // create a contiguous view to the relevant part of the source multivector
00141         cc = mv.subView(rng);
00142       }
00143       else {
00144         // create an indexed view to the relevant part of the source multivector
00145         cc = mv.subView( numvecs, &(index[0]) );
00146       }
00147       return cc;
00148     }
00149 
00155     static Teuchos::RCP<const  Thyra::MultiVectorBase< ScalarType > > CloneView( const  Thyra::MultiVectorBase< ScalarType > & mv, const std::vector<int>& index )
00156     {
00157       int numvecs = index.size();
00158 
00159       // We do not assume that the indices are sorted, nor do we check that
00160       // index.size() > 0. This code is fail-safe, in the sense that a zero
00161       // length index vector will pass the error on the Thyra.
00162 
00163       // Thyra has two ways to create an indexed View:
00164       // * contiguous (via a range of columns)
00165       // * indexed (via a vector of column indices)
00166       // The former is significantly more efficient than the latter, in terms of
00167       // computations performed with/against the created view.
00168       // We will therefore check to see if the given indices are contiguous, and
00169       // if so, we will use the contiguous view creation method.
00170 
00171       int lb = index[0];
00172       bool contig = true;
00173       for (int i=0; i<numvecs; i++) {
00174         if (lb+i != index[i]) contig = false;
00175       }
00176 
00177       Teuchos::RCP< const  Thyra::MultiVectorBase< ScalarType >  > cc;
00178       if (contig) {
00179         const Thyra::Range1D rng(lb,lb+numvecs-1);
00180         // create a contiguous view to the relevant part of the source multivector
00181         cc = mv.subView(rng);
00182       }
00183       else {
00184         // create an indexed view to the relevant part of the source multivector
00185         cc = mv.subView( numvecs, &(index[0]) );
00186       }
00187       return cc;
00188     }
00189 
00191 
00194 
00196     static int GetVecLength( const  Thyra::MultiVectorBase< ScalarType > & mv )
00197     { return mv.range()->dim(); }
00198 
00200     static int GetNumberVecs( const  Thyra::MultiVectorBase< ScalarType > & mv )
00201     { return mv.domain()->dim(); }
00202 
00204 
00207 
00210     static void MvTimesMatAddMv( const ScalarType alpha, const  Thyra::MultiVectorBase< ScalarType > & A, 
00211          const Teuchos::SerialDenseMatrix<int,ScalarType>& B, 
00212          const ScalarType beta,  Thyra::MultiVectorBase< ScalarType > & mv )
00213     {
00214       int m = B.numRows();
00215       int n = B.numCols();
00216       // Create a view of the B object!
00217       Teuchos::RCP< const  Thyra::MultiVectorBase< ScalarType >  >
00218         B_thyra = Thyra::createMembersView(
00219           A.domain()
00220           ,RTOpPack::ConstSubMultiVectorView<ScalarType>(0,m,0,n,&B(0,0),B.stride())
00221           );
00222       // perform the operation via A: mv <- alpha*A*B_thyra + beta*mv
00223       A.apply(Thyra::NONCONJ_ELE,*B_thyra,&mv,alpha,beta);
00224     }
00225 
00228     static void MvAddMv( const ScalarType alpha, const  Thyra::MultiVectorBase< ScalarType > & A, 
00229                          const ScalarType beta,  const  Thyra::MultiVectorBase< ScalarType > & B,  Thyra::MultiVectorBase< ScalarType > & mv )
00230     { 
00231       ScalarType coef[2], zero = Teuchos::ScalarTraits<ScalarType>::zero();
00232       const  Thyra::MultiVectorBase< ScalarType >  * in[2];
00233 
00234       in[0] = &A;
00235       in[1] = &B;
00236       coef[0] = alpha;
00237       coef[1] = beta;
00238 
00239       Thyra::linear_combination(2,coef,in,zero,&mv);
00240     }
00241 
00244     static void MvTransMv( const ScalarType alpha, const  Thyra::MultiVectorBase< ScalarType > & A, const  Thyra::MultiVectorBase< ScalarType > & mv, 
00245                            Teuchos::SerialDenseMatrix<int,ScalarType>& B )
00246     { 
00247       // Create a multivector to hold the result (m by n)
00248       int m = A.domain()->dim();
00249       int n = mv.domain()->dim();
00250       // Create a view of the B object!
00251       Teuchos::RCP<  Thyra::MultiVectorBase< ScalarType >  >
00252         B_thyra = Thyra::createMembersView(
00253           A.domain()
00254           ,RTOpPack::SubMultiVectorView<ScalarType>(0,m,0,n,&B(0,0),B.stride())
00255           );
00256       A.applyTranspose(Thyra::CONJ_ELE,mv,&*B_thyra,alpha,Teuchos::ScalarTraits<ScalarType>::zero());
00257     }
00258 
00262     static void MvDot( const  Thyra::MultiVectorBase< ScalarType > & mv, const  Thyra::MultiVectorBase< ScalarType > & A, std::vector<ScalarType> &b )
00263     { Thyra::dots(mv,A,&(b[0])); }
00264 
00267     static void MvScale (  Thyra::MultiVectorBase< ScalarType > & mv, const ScalarType alpha )
00268     { Thyra::scale(alpha,&mv); }
00269     
00272     static void MvScale (  Thyra::MultiVectorBase< ScalarType > & mv, const std::vector<ScalarType>& alpha ) 
00273     { 
00274       for (unsigned int i=0; i<alpha.size(); i++) {
00275         Thyra::scale(alpha[i],mv.col(i).get());
00276       }
00277     }
00278     
00280 
00283 
00287     static void MvNorm( const  Thyra::MultiVectorBase< ScalarType > & mv, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec )
00288     { Thyra::norms_2(mv,&(normvec[0])); }
00289 
00291 
00294 
00297     static void SetBlock( const  Thyra::MultiVectorBase< ScalarType > & A, const std::vector<int>& index,  Thyra::MultiVectorBase< ScalarType > & mv )
00298     { 
00299       // Extract the "numvecs" columns of mv indicated by the index vector.
00300       int numvecs = index.size();
00301       std::vector<int> indexA(numvecs);
00302       int numAcols = A.domain()->dim();
00303       for (int i=0; i<numvecs; i++) {
00304         indexA[i] = i;
00305       }
00306       // Thyra::assign requires that both arguments have the same number of 
00307       // vectors. Enforce this, by shrinking one to match the other.
00308       if ( numAcols < numvecs ) {
00309         // A does not have enough columns to satisfy index_plus. Shrink
00310         // index_plus.
00311         numvecs = numAcols;
00312       }
00313       else if ( numAcols > numvecs ) {
00314         numAcols = numvecs;
00315         indexA.resize( numAcols );
00316       }
00317       // create a view to the relevant part of the source multivector
00318       Teuchos::RCP< const  Thyra::MultiVectorBase< ScalarType >  > relsource = A.subView( numAcols, &(indexA[0]) );
00319       // create a view to the relevant part of the destination multivector
00320       Teuchos::RCP<  Thyra::MultiVectorBase< ScalarType >  > reldest = mv.subView( numvecs, &(index[0]) );
00321       // copy the data to the destination multivector subview
00322       Thyra::assign(&*reldest, *relsource);
00323     }
00324 
00327     static void MvRandom(  Thyra::MultiVectorBase< ScalarType > & mv )
00328     { 
00329       // Thyra::randomize generates via a uniform distribution on [l,u]
00330       // We will use this to generate on [-1,1]
00331       Thyra::randomize(-Teuchos::ScalarTraits<ScalarType>::one(),
00332                         Teuchos::ScalarTraits<ScalarType>::one(),
00333                        &mv);
00334     }
00335 
00338     static void MvInit(  Thyra::MultiVectorBase< ScalarType > & mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
00339     { Thyra::assign(&mv,alpha); }
00340 
00342 
00345 
00348     static void MvPrint( const  Thyra::MultiVectorBase< ScalarType > & mv, std::ostream& os )
00349     { 
00350        Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&os,false));    
00351        out->setf(std::ios_base::scientific);    
00352        out->precision(16);   
00353        mv.describe(*out,Teuchos::VERB_EXTREME);
00354     }
00355 
00357 
00358   };        
00359 
00361   //
00362   // Implementation of the Anasazi::OperatorTraits for Thyra::LinearOpBase
00363   //
00365 
00375   template <class ScalarType> 
00376   class OperatorTraits < ScalarType, Thyra::MultiVectorBase<ScalarType>, Thyra::LinearOpBase<ScalarType> >
00377   {
00378   public:
00379     
00383     static void Apply ( const Thyra::LinearOpBase< ScalarType >& Op, const  Thyra::MultiVectorBase< ScalarType > & x,  Thyra::MultiVectorBase< ScalarType > & y )
00384     { 
00385       Op.apply(Thyra::NONCONJ_ELE,x,&y);
00386     }
00387     
00388   };
00389   
00390 } // end of Anasazi namespace 
00391 
00392 #endif 
00393 // end of file ANASAZI_THYRA_ADAPTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends