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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
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 #include <Thyra_VectorStdOps.hpp>
00045 
00046 #include <Teuchos_Ptr.hpp>
00047 #include <Teuchos_ArrayRCP.hpp>
00048 #include <Teuchos_ArrayView.hpp>
00049 
00050 namespace Anasazi {
00051   
00053   //
00054   // Implementation of the Anasazi::MultiVecTraits for Thyra::MultiVectorBase
00055   //
00057 
00066   template<class ScalarType>
00067   class MultiVecTraits< ScalarType, Thyra::MultiVectorBase<ScalarType> >
00068   {
00069   private:
00070     typedef Thyra::MultiVectorBase<ScalarType> TMVB;
00071     typedef Teuchos::ScalarTraits<ScalarType> ST;
00072     typedef typename ST::magnitudeType magType;
00073 
00074   public:
00075 
00078 
00083     static Teuchos::RCP<TMVB> Clone( const TMVB & mv, const int numvecs )
00084     { 
00085       Teuchos::RCP<TMVB> c = Thyra::createMembers( mv.range(), numvecs ); 
00086       return c;
00087     }
00088 
00093     static Teuchos::RCP<TMVB> 
00094     CloneCopy (const TMVB& mv)
00095     { 
00096       const int numvecs = mv.domain()->dim();
00097       // create the new multivector
00098       Teuchos::RCP< TMVB > cc = Thyra::createMembers (mv.range(), numvecs);
00099       // copy the data from the source multivector to the new multivector
00100       Thyra::assign (Teuchos::outArg (*cc), mv);
00101       return cc;
00102     }
00103 
00109     static Teuchos::RCP< TMVB > CloneCopy( const TMVB & mv, const std::vector<int>& index )
00110     { 
00111       const int numvecs = index.size();
00112       // create the new multivector
00113       Teuchos::RCP<TMVB > cc = Thyra::createMembers (mv.range(), numvecs);
00114       // create a view to the relevant part of the source multivector
00115       Teuchos::RCP<const TMVB > view = mv.subView ( Teuchos::arrayViewFromVector( index ) );
00116       // copy the data from the relevant view to the new multivector
00117       Thyra::assign (Teuchos::outArg (*cc), *view);
00118       return cc;
00119     }
00120 
00121     static Teuchos::RCP<TMVB>
00122     CloneCopy (const TMVB& mv, const Teuchos::Range1D& index)
00123     { 
00124       const int numVecs = index.size();      
00125       // Create the new multivector
00126       Teuchos::RCP<TMVB> cc = Thyra::createMembers (mv.range(), numVecs);
00127       // Create a view to the relevant part of the source multivector
00128       Teuchos::RCP<const TMVB> view = mv.subView (index);
00129       // Copy the data from the view to the new multivector.
00130       Thyra::assign (Teuchos::outArg (*cc), *view);
00131       return cc;
00132     }
00133 
00139     static Teuchos::RCP< TMVB > CloneViewNonConst( TMVB & mv, const std::vector<int>& index )
00140     {
00141       int numvecs = index.size();
00142 
00143       // We do not assume that the indices are sorted, nor do we check that
00144       // index.size() > 0. This code is fail-safe, in the sense that a zero
00145       // length index vector will pass the error on the Thyra.
00146 
00147       // Thyra has two ways to create an indexed View:
00148       // * contiguous (via a range of columns)
00149       // * indexed (via a vector of column indices)
00150       // The former is significantly more efficient than the latter, in terms of
00151       // computations performed with/against the created view.
00152       // We will therefore check to see if the given indices are contiguous, and
00153       // if so, we will use the contiguous view creation method.
00154 
00155       int lb = index[0];
00156       bool contig = true;
00157       for (int i=0; i<numvecs; i++) {
00158         if (lb+i != index[i]) contig = false;
00159       }
00160 
00161       Teuchos::RCP< TMVB > cc;
00162       if (contig) {
00163         const Thyra::Range1D rng(lb,lb+numvecs-1);
00164         // create a contiguous view to the relevant part of the source multivector
00165         cc = mv.subView(rng);
00166       }
00167       else {
00168         // create an indexed view to the relevant part of the source multivector
00169         cc = mv.subView( Teuchos::arrayViewFromVector( index ) );
00170       }
00171       return cc;
00172     }
00173 
00174     static Teuchos::RCP<TMVB> 
00175     CloneViewNonConst (TMVB& mv, const Teuchos::Range1D& index)
00176     {
00177       // We let Thyra be responsible for checking that the index range
00178       // is nonempty.
00179       //
00180       // Create and return a contiguous view to the relevant part of
00181       // the source multivector.
00182       return mv.subView (index);
00183     }
00184 
00190     static Teuchos::RCP<const TMVB > CloneView( const TMVB & mv, const std::vector<int>& index )
00191     {
00192       int numvecs = index.size();
00193 
00194       // We do not assume that the indices are sorted, nor do we check that
00195       // index.size() > 0. This code is fail-safe, in the sense that a zero
00196       // length index vector will pass the error on the Thyra.
00197 
00198       // Thyra has two ways to create an indexed View:
00199       // * contiguous (via a range of columns)
00200       // * indexed (via a vector of column indices)
00201       // The former is significantly more efficient than the latter, in terms of
00202       // computations performed with/against the created view.
00203       // We will therefore check to see if the given indices are contiguous, and
00204       // if so, we will use the contiguous view creation method.
00205 
00206       int lb = index[0];
00207       bool contig = true;
00208       for (int i=0; i<numvecs; i++) {
00209         if (lb+i != index[i]) contig = false;
00210       }
00211 
00212       Teuchos::RCP< const TMVB > cc;
00213       if (contig) {
00214         const Thyra::Range1D rng(lb,lb+numvecs-1);
00215         // create a contiguous view to the relevant part of the source multivector
00216         cc = mv.subView(rng);
00217       }
00218       else {
00219         // create an indexed view to the relevant part of the source multivector
00220         cc = mv.subView(Teuchos::arrayViewFromVector( index ) );
00221       }
00222       return cc;
00223     }
00224 
00225     static Teuchos::RCP<const TMVB> 
00226     CloneView (const TMVB& mv, const Teuchos::Range1D& index)
00227     {
00228       // We let Thyra be responsible for checking that the index range
00229       // is nonempty.
00230       //
00231       // Create and return a contiguous view to the relevant part of
00232       // the source multivector.
00233       return mv.subView (index);
00234     }
00235 
00237 
00240 
00242     static int GetVecLength( const TMVB & mv )
00243     { return mv.range()->dim(); }
00244 
00246     static int GetNumberVecs( const TMVB & mv )
00247     { return mv.domain()->dim(); }
00248 
00250 
00253 
00256     static void MvTimesMatAddMv( const ScalarType alpha, const TMVB & A, 
00257          const Teuchos::SerialDenseMatrix<int,ScalarType>& B, 
00258          const ScalarType beta, TMVB & mv )
00259     {
00260       int m = B.numRows();
00261       int n = B.numCols();
00262       // Create a view of the B object!
00263       Teuchos::RCP< const TMVB >
00264         B_thyra = Thyra::createMembersView(
00265           A.domain(),
00266           RTOpPack::ConstSubMultiVectorView<ScalarType>(
00267             0, m, 0, n,
00268             Teuchos::arcpFromArrayView(Teuchos::arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
00269             )
00270           );
00271       // perform the operation via A: mv <- alpha*A*B_thyra + beta*mv
00272       A.apply(Thyra::NOTRANS,*B_thyra,Teuchos::outArg (mv),alpha,beta);
00273     }
00274 
00277     static void MvAddMv( const ScalarType alpha, const TMVB & A, 
00278                          const ScalarType beta,  const TMVB & B, TMVB & mv )
00279     { 
00280       using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::inoutArg;
00281 
00282       Thyra::linear_combination<ScalarType>(
00283         tuple(alpha, beta)(), tuple(ptrInArg(A), ptrInArg(B))(), ST::zero(), inoutArg(mv));
00284     }
00285 
00288     static void MvTransMv( const ScalarType alpha, const TMVB & A, const TMVB & mv, 
00289                            Teuchos::SerialDenseMatrix<int,ScalarType>& B )
00290     { 
00291       // Create a multivector to hold the result (m by n)
00292       int m = A.domain()->dim();
00293       int n = mv.domain()->dim();
00294       // Create a view of the B object!
00295       Teuchos::RCP< TMVB >
00296         B_thyra = Thyra::createMembersView(
00297           A.domain(),
00298           RTOpPack::SubMultiVectorView<ScalarType>(
00299             0, m, 0, n,
00300             Teuchos::arcpFromArrayView(Teuchos::arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
00301             )
00302           );
00303       A.apply(Thyra::CONJTRANS,mv,B_thyra.ptr(),alpha,Teuchos::ScalarTraits<ScalarType>::zero());
00304     }
00305 
00309     static void MvDot( const TMVB & mv, const TMVB & A, std::vector<ScalarType> &b )
00310     { Thyra::dots(mv,A,Teuchos::arrayViewFromVector( b )); }
00311 
00314     static void 
00315     MvScale (TMVB& mv, 
00316        const ScalarType alpha)
00317     { 
00318       Thyra::scale (alpha, Teuchos::inOutArg (mv)); 
00319     }
00320     
00323     static void 
00324     MvScale (TMVB& mv, 
00325        const std::vector<ScalarType>& alpha) 
00326     { 
00327       for (unsigned int i=0; i<alpha.size(); i++) {
00328         Thyra::scale (alpha[i], mv.col(i).ptr());
00329       }
00330     }
00331     
00333 
00336 
00340     static void MvNorm( const TMVB & mv, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec )
00341     { Thyra::norms_2(mv,Teuchos::arrayViewFromVector( normvec )); }
00342 
00344 
00347 
00350     static void SetBlock( const TMVB & A, const std::vector<int>& index, TMVB & mv )
00351     { 
00352       // Extract the "numvecs" columns of mv indicated by the index vector.
00353       int numvecs = index.size();
00354       std::vector<int> indexA(numvecs);
00355       int numAcols = A.domain()->dim();
00356       for (int i=0; i<numvecs; i++) {
00357         indexA[i] = i;
00358       }
00359       // Thyra::assign requires that both arguments have the same number of 
00360       // vectors. Enforce this, by shrinking one to match the other.
00361       if ( numAcols < numvecs ) {
00362         // A does not have enough columns to satisfy index_plus. Shrink
00363         // index_plus.
00364         numvecs = numAcols;
00365       }
00366       else if ( numAcols > numvecs ) {
00367         numAcols = numvecs;
00368         indexA.resize( numAcols );
00369       }
00370       // create a view to the relevant part of the source multivector
00371       Teuchos::RCP< const TMVB > relsource = A.subView( Teuchos::arrayViewFromVector( indexA ) );
00372       // create a view to the relevant part of the destination multivector
00373       Teuchos::RCP< TMVB > reldest = mv.subView( Teuchos::arrayViewFromVector( index ) );
00374       // copy the data to the destination multivector subview
00375       Thyra::assign (Teuchos::outArg (*reldest), *relsource);
00376     }
00377 
00378     static void 
00379     SetBlock (const TMVB& A, const Teuchos::Range1D& index, TMVB& mv)
00380     { 
00381       const int numColsA = A.domain()->dim();
00382       const int numColsMv = mv.domain()->dim();
00383       // 'index' indexes into mv; it's the index set of the target.
00384       const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
00385       // We can't take more columns out of A than A has.
00386       const bool validSource = index.size() <= numColsA;
00387 
00388       if (! validIndex || ! validSource)
00389   {
00390     std::ostringstream os;
00391     os << "Anasazi::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar> "
00392       ">::SetBlock(A, [" << index.lbound() << ", " << index.ubound() 
00393        << "], mv): ";
00394     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00395            os.str() << "Range lower bound must be nonnegative.");
00396     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
00397            os.str() << "Range upper bound must be less than "
00398            "the number of columns " << numColsA << " in the "
00399            "'mv' output argument.");
00400     TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
00401            os.str() << "Range must have no more elements than"
00402            " the number of columns " << numColsA << " in the "
00403            "'A' input argument.");
00404     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
00405   }
00406 
00407       // View of the relevant column(s) of the target multivector mv.
00408       // We avoid view creation overhead by only creating a view if
00409       // the index range is different than [0, (# columns in mv) - 1].
00410       Teuchos::RCP<TMVB> mv_view;
00411       if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
00412   mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
00413       else
00414   mv_view = mv.subView (index);
00415 
00416       // View of the relevant column(s) of the source multivector A.
00417       // If A has fewer columns than mv_view, then create a view of
00418       // the first index.size() columns of A.
00419       Teuchos::RCP<const TMVB> A_view;
00420       if (index.size() == numColsA)
00421   A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
00422       else
00423   A_view = A.subView (Teuchos::Range1D(0, index.size()-1));
00424 
00425       // Copy the data to the destination multivector.
00426       Thyra::assign (Teuchos::outArg (*mv_view), *A_view);
00427     }
00428 
00429     static void 
00430     Assign (const TMVB& A, TMVB& mv)
00431     { 
00432       using Teuchos::ptr;
00433       using Teuchos::Range1D;
00434       using Teuchos::RCP;
00435 
00436       const int numColsA = A.domain()->dim();
00437       const int numColsMv = mv.domain()->dim();
00438       if (numColsA > numColsMv) {
00439   const std::string prefix ("Anasazi::MultiVecTraits<Scalar, "
00440           "Thyra::MultiVectorBase<Scalar>"
00441           " >::Assign(A, mv): ");
00442   TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
00443            prefix << "Input multivector 'A' has " 
00444            << numColsA << " columns, but output multivector "
00445            "'mv' has only " << numColsMv << " columns.");
00446   TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
00447       }
00448       // Copy the data to the destination multivector.
00449       if (numColsA == numColsMv) {
00450   Thyra::assign (outArg (mv), A);
00451       } 
00452       else {
00453   RCP<TMVB> mv_view = CloneViewNonConst (mv, Range1D(0, numColsA-1));
00454   Thyra::assign (outArg (*mv_view), A);
00455       }
00456     }
00457 
00460     static void MvRandom( TMVB & mv )
00461     { 
00462       // Thyra::randomize generates via a uniform distribution on [l,u]
00463       // We will use this to generate on [-1,1]
00464       Thyra::randomize(-Teuchos::ScalarTraits<ScalarType>::one(),
00465                         Teuchos::ScalarTraits<ScalarType>::one(),
00466                         Teuchos::outArg (mv));
00467     }
00468 
00471     static void 
00472     MvInit (TMVB& mv, 
00473       ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero())
00474     { 
00475       Thyra::assign (Teuchos::outArg (mv), alpha); 
00476     }
00477 
00479 
00482 
00485     static void MvPrint( const TMVB & mv, std::ostream& os )
00486     { 
00487        Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&os,false));    
00488        out->setf(std::ios_base::scientific);    
00489        out->precision(16);   
00490        mv.describe(*out,Teuchos::VERB_EXTREME);
00491     }
00492 
00494 
00495   };        
00496 
00498   //
00499   // Implementation of the Anasazi::OperatorTraits for Thyra::LinearOpBase
00500   //
00502 
00512   template <class ScalarType> 
00513   class OperatorTraits < ScalarType, Thyra::MultiVectorBase<ScalarType>, Thyra::LinearOpBase<ScalarType> >
00514   {
00515   public:
00516     
00520     static void Apply ( const Thyra::LinearOpBase< ScalarType >& Op, const  Thyra::MultiVectorBase< ScalarType > & x,  Thyra::MultiVectorBase< ScalarType > & y )
00521     { 
00522       Op.apply(Thyra::NOTRANS,x,Teuchos::outArg (y), Teuchos::ScalarTraits<ScalarType>::one(),Teuchos::ScalarTraits<ScalarType>::zero());
00523     }
00524     
00525   };
00526   
00527 } // end of Anasazi namespace 
00528 
00529 #endif 
00530 // end of file ANASAZI_THYRA_ADAPTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends