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_MultiVectorBaseDecl.hpp>
00043 
00044 namespace Anasazi {
00045   
00047   //
00048   // Implementation of the Anasazi::MultiVecTraits for Thyra::MultiVectorBase
00049   //
00051 
00060   template<class ScalarType>
00061   class MultiVecTraits< ScalarType, Thyra::MultiVectorBase<ScalarType> >
00062   {
00063   private:
00064     typedef Thyra::MultiVectorBase<ScalarType> TMVB;
00065   public:
00066 
00069 
00074     static Teuchos::RefCountPtr<TMVB> Clone( const TMVB& mv, const int numvecs )
00075     { 
00076       Teuchos::RefCountPtr<TMVB> c = Thyra::createMembers( mv.range(), numvecs ); 
00077       return c;
00078     }
00079 
00084     static Teuchos::RefCountPtr<TMVB> CloneCopy( const TMVB& mv )
00085     { 
00086       int numvecs = mv.domain()->dim();
00087       // create the new multivector
00088       Teuchos::RefCountPtr< TMVB > cc = Thyra::createMembers( mv.range(), numvecs );
00089       // copy the data from the source multivector to the new multivector
00090       Thyra::assign(&*cc, mv);
00091       return cc;
00092     }
00093 
00099     static Teuchos::RefCountPtr<TMVB> CloneCopy( const TMVB& mv, const std::vector<int>& index )
00100     { 
00101       int numvecs = index.size();
00102       // create the new multivector
00103       Teuchos::RefCountPtr< TMVB > cc = Thyra::createMembers( mv.range(), numvecs );
00104       // create a view to the relevant part of the source multivector
00105       Teuchos::RefCountPtr< const TMVB > view = mv.subView( numvecs, &(index[0]) );
00106       // copy the data from the relevant view to the new multivector
00107       Thyra::assign(&*cc, *view);
00108       return cc;
00109     }
00110 
00116     static Teuchos::RefCountPtr<TMVB> CloneView( TMVB& mv, const std::vector<int>& index )
00117     {
00118       int numvecs = index.size();
00119 
00120       // We do not assume that the indices are sorted, nor do we check that
00121       // index.size() > 0. This code is fail-safe, in the sense that a zero
00122       // length index vector will pass the error on the Thyra.
00123 
00124       // Thyra has two ways to create an indexed View:
00125       // * contiguous (via a range of columns)
00126       // * indexed (via a vector of column indices)
00127       // The former is significantly more efficient than the latter, in terms of
00128       // computations performed with/against the created view.
00129       // We will therefore check to see if the given indices are contiguous, and
00130       // if so, we will use the contiguous view creation method.
00131 
00132       int lb = index[0];
00133       bool contig = true;
00134       for (int i=0; i<numvecs; i++) {
00135         if (lb+i != index[i]) contig = false;
00136       }
00137 
00138       Teuchos::RefCountPtr< TMVB > cc;
00139       if (contig) {
00140         const Thyra::Range1D rng(lb,lb+numvecs-1);
00141         // create a contiguous view to the relevant part of the source multivector
00142         cc = mv.subView(rng);
00143       }
00144       else {
00145         // create an indexed view to the relevant part of the source multivector
00146         cc = mv.subView( numvecs, &(index[0]) );
00147       }
00148       return cc;
00149     }
00150 
00156     static Teuchos::RefCountPtr<const TMVB> CloneView( const TMVB& mv, const std::vector<int>& index )
00157     {
00158       int numvecs = index.size();
00159 
00160       // We do not assume that the indices are sorted, nor do we check that
00161       // index.size() > 0. This code is fail-safe, in the sense that a zero
00162       // length index vector will pass the error on the Thyra.
00163 
00164       // Thyra has two ways to create an indexed View:
00165       // * contiguous (via a range of columns)
00166       // * indexed (via a vector of column indices)
00167       // The former is significantly more efficient than the latter, in terms of
00168       // computations performed with/against the created view.
00169       // We will therefore check to see if the given indices are contiguous, and
00170       // if so, we will use the contiguous view creation method.
00171 
00172       int lb = index[0];
00173       bool contig = true;
00174       for (int i=0; i<numvecs; i++) {
00175         if (lb+i != index[i]) contig = false;
00176       }
00177 
00178       Teuchos::RefCountPtr< const TMVB > cc;
00179       if (contig) {
00180         const Thyra::Range1D rng(lb,lb+numvecs-1);
00181         // create a contiguous view to the relevant part of the source multivector
00182         cc = mv.subView(rng);
00183       }
00184       else {
00185         // create an indexed view to the relevant part of the source multivector
00186         cc = mv.subView( numvecs, &(index[0]) );
00187       }
00188       return cc;
00189     }
00190 
00192 
00195 
00197     static int GetVecLength( const TMVB& mv )
00198     { return mv.range()->dim(); }
00199 
00201     static int GetNumberVecs( const TMVB& mv )
00202     { return mv.domain()->dim(); }
00203 
00205 
00208 
00211     static void MvTimesMatAddMv( const ScalarType alpha, const TMVB& A, 
00212          const Teuchos::SerialDenseMatrix<int,ScalarType>& B, 
00213          const ScalarType beta, TMVB& mv )
00214     {
00215       int m = B.numRows();
00216       int n = B.numCols();
00217       // Create a view of the B object!
00218       Teuchos::RefCountPtr< const TMVB >
00219         B_thyra = Thyra::createMembersView(
00220           A.domain()
00221           ,RTOpPack::ConstSubMultiVectorView<ScalarType>(0,m,0,n,&B(0,0),B.stride())
00222           );
00223       // perform the operation via A: mv <- alpha*A*B_thyra + beta*mv
00224       A.apply(Thyra::NONCONJ_ELE,*B_thyra,&mv,alpha,beta);
00225     }
00226 
00229     static void MvAddMv( const ScalarType alpha, const TMVB& A, 
00230                          const ScalarType beta,  const TMVB& B, TMVB& mv )
00231     { 
00232       ScalarType coef[2], zero = Teuchos::ScalarTraits<ScalarType>::zero();
00233       const TMVB * in[2];
00234 
00235       in[0] = &A;
00236       in[1] = &B;
00237       coef[0] = alpha;
00238       coef[1] = beta;
00239 
00240       Thyra::linear_combination(2,coef,in,zero,&mv);
00241     }
00242 
00245     static void MvTransMv( const ScalarType alpha, const TMVB& A, const TMVB& mv, 
00246                            Teuchos::SerialDenseMatrix<int,ScalarType>& B )
00247     { 
00248       // Create a multivector to hold the result (m by n)
00249       int m = A.domain()->dim();
00250       int n = mv.domain()->dim();
00251       // Create a view of the B object!
00252       Teuchos::RefCountPtr< TMVB >
00253         B_thyra = Thyra::createMembersView(
00254           A.domain()
00255           ,RTOpPack::SubMultiVectorView<ScalarType>(0,m,0,n,&B(0,0),B.stride())
00256           );
00257       A.applyTranspose(Thyra::CONJ_ELE,mv,&*B_thyra,alpha,Teuchos::ScalarTraits<ScalarType>::zero());
00258     }
00259 
00263     static void MvDot( const TMVB& mv, const TMVB& A, std::vector<ScalarType>* b )
00264     { Thyra::dots(mv,A,&((*b)[0])); }
00265 
00268     static void MvScale ( TMVB& mv, const ScalarType alpha )
00269     { Thyra::scale(alpha,&mv); }
00270     
00273     static void MvScale ( TMVB& mv, const std::vector<ScalarType>& alpha ) 
00274     { 
00275       for (unsigned int i=0; i<alpha.size(); i++) {
00276         Thyra::scale(alpha[i],mv.col(i).get());
00277       }
00278     }
00279     
00281 
00284 
00288     static void MvNorm( const TMVB& mv, std::vector<ScalarType>* normvec )
00289     { Thyra::norms_2(mv,&((*normvec)[0])); }
00290 
00292 
00295 
00298     static void SetBlock( const TMVB& A, const std::vector<int>& index, TMVB& mv )
00299     { 
00300       // Extract the "numvecs" columns of mv indicated by the index vector.
00301       int numvecs = index.size();
00302       std::vector<int> indexA(numvecs);
00303       int numAcols = A.domain()->dim();
00304       for (int i=0; i<numvecs; i++) {
00305         indexA[i] = i;
00306       }
00307       // Thyra::assign requires that both arguments have the same number of 
00308       // vectors. Enforce this, by shrinking one to match the other.
00309       if ( numAcols < numvecs ) {
00310         // A does not have enough columns to satisfy index_plus. Shrink
00311         // index_plus.
00312         numvecs = numAcols;
00313       }
00314       else if ( numAcols > numvecs ) {
00315         numAcols = numvecs;
00316         indexA.resize( numAcols );
00317       }
00318       // create a view to the relevant part of the source multivector
00319       Teuchos::RefCountPtr< const TMVB > relsource = A.subView( numAcols, &(indexA[0]) );
00320       // create a view to the relevant part of the destination multivector
00321       Teuchos::RefCountPtr< TMVB > reldest = mv.subView( numvecs, &(index[0]) );
00322       // copy the data to the destination multivector subview
00323       Thyra::assign(&*reldest, *relsource);
00324     }
00325 
00328     static void MvRandom( TMVB& mv )
00329     { 
00330       // Thyra::randomize generates via a uniform distribution on [l,u]
00331       // We will use this to generate on [-1,1]
00332       Thyra::randomize(-Teuchos::ScalarTraits<ScalarType>::one(),
00333                         Teuchos::ScalarTraits<ScalarType>::one(),
00334                        &mv);
00335     }
00336 
00339     static void MvInit( TMVB& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
00340     { Thyra::assign(&mv,alpha); }
00341 
00343 
00346 
00349     static void MvPrint( const TMVB& mv, ostream& os )
00350     { }
00351 
00353 
00354   };        
00355 
00357   //
00358   // Implementation of the Anasazi::OperatorTraits for Thyra::LinearOpBase
00359   //
00361 
00371   template <class ScalarType> 
00372   class OperatorTraits < ScalarType, Thyra::MultiVectorBase<ScalarType>, Thyra::LinearOpBase<ScalarType> >
00373   {
00374   private:
00375     typedef Thyra::MultiVectorBase<ScalarType> TMVB;
00376     typedef Thyra::LinearOpBase<ScalarType>    TLOB;
00377   public:
00378     
00382     static void Apply ( const TLOB& Op, const TMVB& x, TMVB& y )
00383     { 
00384       Op.apply(Thyra::NONCONJ_ELE,x,&y);
00385     }
00386     
00387   };
00388   
00389 } // end of Anasazi namespace 
00390 
00391 #endif 
00392 // end of file ANASAZI_THYRA_ADAPTER_HPP

Generated on Thu Sep 18 12:31:38 2008 for Anasazi by doxygen 1.3.9.1