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