Belos Version of the Day
BelosTsqrOrthoManager.hpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                  Copyright 2004 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00044 
00045 #ifndef __BelosTsqrOrthoManager_hpp
00046 #define __BelosTsqrOrthoManager_hpp
00047 
00048 #include "BelosTsqrOrthoManagerImpl.hpp"
00049 // Belos doesn't have an SVQB orthomanager implemented yet, so we just
00050 // use a default orthomanager for the case where the inner product
00051 // matrix is nontrivial.
00052 #include "BelosDGKSOrthoManager.hpp"
00053 
00054 
00055 namespace Belos {
00056 
00078   template<class Scalar, class MV>
00079     class OutOfPlaceNormalizerMixin {
00080       public:
00081         typedef Scalar scalar_type;
00082         typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00085         typedef MV multivector_type;
00086 
00087         typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00088         typedef Teuchos::RCP<mat_type>                  mat_ptr;
00089 
00098         virtual int
00099           normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const = 0;
00100 
00119         virtual int
00120           projectAndNormalizeOutOfPlace (MV& X_in,
00121               MV& X_out,
00122               Teuchos::Array<mat_ptr> C,
00123               mat_ptr B,
00124               Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const = 0;
00125 
00127         virtual ~OutOfPlaceNormalizerMixin () {}
00128     };
00129 
00136   template<class Scalar, class MV>
00137     class TsqrOrthoManager :
00138       public OrthoManager<Scalar, MV>,
00139       public OutOfPlaceNormalizerMixin<Scalar, MV>,
00140       public Teuchos::ParameterListAcceptor
00141   {
00142     public:
00143       typedef Scalar scalar_type;
00144       typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00146       typedef MV multivector_type;
00147 
00148       typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00149       typedef Teuchos::RCP<mat_type>                  mat_ptr;
00150 
00151       void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
00152         impl_.setParameterList (params);
00153       }
00154 
00155       Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList () {
00156         return impl_.getNonconstParameterList ();
00157       }
00158 
00159       Teuchos::RCP<Teuchos::ParameterList> unsetParameterList () {
00160         return impl_.unsetParameterList ();
00161       }
00162 
00170       Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
00171         return impl_.getValidParameters();
00172       }
00173 
00183       Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() const {
00184         return impl_.getFastParameters();
00185       }
00186 
00203       TsqrOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params,
00204           const std::string& label = "Belos") :
00205         impl_ (params, label)
00206         {}
00207 
00212       TsqrOrthoManager (const std::string& label) :
00213         impl_ (label)
00214     {}
00215 
00217       virtual ~TsqrOrthoManager() {}
00218 
00240       void
00241         setReorthogonalizationCallback (const Teuchos::RCP<ReorthogonalizationCallback<Scalar> >& callback)
00242         {
00243           impl_.setReorthogonalizationCallback (callback);
00244         }
00245 
00246       void innerProd (const MV &X, const MV &Y, mat_type& Z) const {
00247         return impl_.innerProd (X, Y, Z);
00248       }
00249 
00250       void norm (const MV& X, std::vector<magnitude_type>& normVec) const {
00251         return impl_.norm (X, normVec);
00252       }
00253 
00254       void
00255         project (MV &X,
00256             Teuchos::Array<mat_ptr> C,
00257             Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00258         {
00259           return impl_.project (X, C, Q);
00260         }
00261 
00262       int
00263         normalize (MV &X, mat_ptr B) const
00264         {
00265           return impl_.normalize (X, B);
00266         }
00267 
00268       int
00269         projectAndNormalize (MV &X,
00270             Teuchos::Array<mat_ptr> C,
00271             mat_ptr B,
00272             Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00273         {
00274           return impl_.projectAndNormalize (X, C, B, Q);
00275         }
00276 
00293       int
00294         normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
00295         {
00296           return impl_.normalizeOutOfPlace (X, Q, B);
00297         }
00298 
00319       int
00320         projectAndNormalizeOutOfPlace (MV& X_in,
00321             MV& X_out,
00322             Teuchos::Array<mat_ptr> C,
00323             mat_ptr B,
00324             Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00325         {
00326           return impl_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
00327         }
00328 
00329       magnitude_type orthonormError (const MV &X) const {
00330         return impl_.orthonormError (X);
00331       }
00332 
00333       magnitude_type orthogError (const MV &X1, const MV &X2) const {
00334         return impl_.orthogError (X1, X2);
00335       }
00336 
00344       void setLabel (const std::string& label) {
00345         impl_.setLabel (label);
00346       }
00347 
00348       const std::string& getLabel() const { return impl_.getLabel(); }
00349 
00350     private:
00356       mutable TsqrOrthoManagerImpl<Scalar, MV> impl_;
00357 
00359       std::string label_;
00360   };
00361 
00362 
00376   template<class Scalar, class MV, class OP>
00377     class TsqrMatOrthoManager :
00378       public MatOrthoManager<Scalar, MV, OP>,
00379       public OutOfPlaceNormalizerMixin<Scalar, MV>,
00380       public Teuchos::ParameterListAcceptorDefaultBase
00381   {
00382     public:
00383       typedef Scalar scalar_type;
00384       typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00386       typedef MV multivector_type;
00388       typedef OP operator_type;
00389 
00390       typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00391       typedef Teuchos::RCP<mat_type>                  mat_ptr;
00392 
00393     private:
00402       typedef MatOrthoManager<Scalar, MV, OP> base_type;
00403 
00406       typedef TsqrOrthoManagerImpl<Scalar, MV> tsqr_type;
00407 
00410       typedef DGKSOrthoManager<Scalar, MV, OP> dgks_type;
00411 
00414       typedef MultiVecTraits<Scalar, MV> MVT;
00415 
00416     public:
00437       TsqrMatOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params,
00438           const std::string& label = "Belos",
00439           Teuchos::RCP<const OP> Op = Teuchos::null) :
00440         MatOrthoManager<Scalar, MV, OP> (Op),
00441         tsqr_ (params, label),
00442         pDgks_ (Teuchos::null) // Initialized lazily
00443     {}
00444 
00453       TsqrMatOrthoManager (const std::string& label = "Belos",
00454           Teuchos::RCP<const OP> Op = Teuchos::null) :
00455         MatOrthoManager<Scalar, MV, OP> (Op),
00456         tsqr_ (label),
00457         pDgks_ (Teuchos::null) // Initialized lazily
00458     {}
00459 
00461       virtual ~TsqrMatOrthoManager() {}
00462 
00470       Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
00471         return tsqr_.getValidParameters ();
00472       }
00473 
00483       Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() {
00484         return tsqr_.getFastParameters ();
00485       }
00486 
00487       void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
00488         tsqr_.setParameterList (params);
00489       }
00490 
00491       const std::string& getLabel() const { return tsqr_.getLabel (); }
00492 
00493       void
00494         setOp (Teuchos::RCP<const OP> Op)
00495         {
00496           // We override the base class' setOp() so that the
00497           // DGKSOrthoManager instance gets the new op.
00498           //
00499           // The base_type notation helps C++ resolve the name for a
00500           // member function of a templated base class.
00501           base_type::setOp (Op); // base class gets a copy of the Op too
00502 
00503           if (! Op.is_null()) {
00504             ensureDgksInit (); // Make sure the DGKS object has been initialized
00505             pDgks_->setOp (Op);
00506           }
00507         }
00508 
00509       Teuchos::RCP<const OP> getOp () const {
00510         // The base_type notation helps C++ resolve the name for a
00511         // member function of a templated base class.
00512         return base_type::getOp();
00513       }
00514 
00515       void
00516         project (MV &X,
00517             Teuchos::RCP<MV> MX,
00518             Teuchos::Array<mat_ptr> C,
00519             Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00520         {
00521           if (getOp().is_null()) {
00522             tsqr_.project (X, C, Q);
00523             if (! MX.is_null()) {
00524               // MX gets a copy of X; M is the identity operator.
00525               MVT::Assign (X, *MX);
00526             }
00527           } else {
00528             ensureDgksInit ();
00529             pDgks_->project (X, MX, C, Q);
00530           }
00531         }
00532 
00533       void
00534         project (MV &X,
00535             Teuchos::Array<mat_ptr> C,
00536             Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00537         {
00538           project (X, Teuchos::null, C, Q);
00539         }
00540 
00541       int
00542         normalize (MV& X, Teuchos::RCP<MV> MX, mat_ptr B) const
00543         {
00544           if (getOp().is_null()) {
00545             const int rank = tsqr_.normalize (X, B);
00546             if (! MX.is_null()) {
00547               // MX gets a copy of X; M is the identity operator.
00548               MVT::Assign (X, *MX);
00549             }
00550             return rank;
00551           } else {
00552             ensureDgksInit ();
00553             return pDgks_->normalize (X, MX, B);
00554           }
00555         }
00556 
00557       int normalize (MV& X, mat_ptr B) const {
00558         return normalize (X, Teuchos::null, B);
00559       }
00560 
00561       // Attempted fix for a warning about hiding
00562       // OrthoManager::projectAndNormalize(). Unfortunately, this fix turns out
00563       // to produce a compilation error with cray++, see bug #6129,
00564       // <https://software.sandia.gov/bugzilla/show_bug.cgi?id=6129>.
00565       //using Belos::OrthoManager<Scalar, MV>::projectAndNormalize;
00566 
00567       int
00568         projectAndNormalize (MV &X,
00569             Teuchos::RCP<MV> MX,
00570             Teuchos::Array<mat_ptr> C,
00571             mat_ptr B,
00572             Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00573         {
00574           if (getOp().is_null()) {
00575             const int rank = tsqr_.projectAndNormalize (X, C, B, Q);
00576             if (! MX.is_null()) {
00577               // MX gets a copy of X; M is the identity operator.
00578               MVT::Assign (X, *MX);
00579             }
00580             return rank;
00581           } else {
00582             ensureDgksInit ();
00583             return pDgks_->projectAndNormalize (X, MX, C, B, Q);
00584           }
00585         }
00586 
00587       int
00588         normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
00589         {
00590           if (getOp().is_null()) {
00591             return tsqr_.normalizeOutOfPlace (X, Q, B);
00592           } else {
00593             // DGKS normalizes in place, so we have to copy.
00594             ensureDgksInit ();
00595             const int rank = pDgks_->normalize (X, B);
00596             MVT::Assign (X, Q);
00597             return rank;
00598           }
00599         }
00600 
00601       int
00602         projectAndNormalizeOutOfPlace (MV& X_in,
00603             MV& X_out,
00604             Teuchos::Array<mat_ptr> C,
00605             mat_ptr B,
00606             Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00607         {
00608           using Teuchos::null;
00609 
00610           if (getOp().is_null()) {
00611             return tsqr_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
00612           } else {
00613             // DGKS normalizes in place, so we have to copy.
00614             ensureDgksInit ();
00615             const int rank = pDgks_->projectAndNormalize (X_in, null, C, B, Q);
00616             MVT::Assign (X_in, X_out);
00617             return rank;
00618           }
00619         }
00620 
00621       magnitude_type
00622         orthonormError (const MV &X, Teuchos::RCP<const MV> MX) const
00623         {
00624           if (getOp().is_null()) {
00625             return tsqr_.orthonormError (X); // Ignore MX
00626           } else {
00627             ensureDgksInit ();
00628             return pDgks_->orthonormError (X, MX);
00629           }
00630         }
00631 
00632       magnitude_type orthonormError (const MV &X) const {
00633         return orthonormError (X, Teuchos::null);
00634       }
00635 
00636       magnitude_type orthogError (const MV &X1, const MV &X2) const {
00637         return orthogError (X1, Teuchos::null, X2);
00638       }
00639 
00640       magnitude_type
00641         orthogError (const MV &X1,
00642             Teuchos::RCP<const MV> MX1,
00643             const MV &X2) const
00644         {
00645           if (getOp().is_null()) {
00646             // Ignore MX1, since we don't need to write to it.
00647             return tsqr_.orthogError (X1, X2);
00648           } else {
00649             ensureDgksInit ();
00650             return pDgks_->orthogError (X1, MX1, X2);
00651           }
00652         }
00653 
00654       void
00655         setLabel (const std::string& label)
00656         {
00657           tsqr_.setLabel (label);
00658 
00659           // Make sure DGKS gets the new label, if it's initialized.
00660           // Otherwise, it will get the new label on initialization.
00661           if (! pDgks_.is_null()) {
00662             pDgks_->setLabel (label);
00663           }
00664         }
00665 
00666     private:
00668       void
00669         ensureDgksInit () const
00670         {
00671           // NOTE (mfh 11 Jan 2011) DGKS has a parameter that needs to be
00672           // set.  For now, we just use the default value of the
00673           // parameter.
00674           if (pDgks_.is_null()) {
00675             pDgks_ = Teuchos::rcp (new dgks_type (getLabel(), getOp()));
00676           }
00677         }
00678 
00686       mutable tsqr_type tsqr_;
00687 
00692       mutable Teuchos::RCP<dgks_type> pDgks_;
00693   };
00694 
00695 } // namespace Belos
00696 
00697 #endif // __BelosTsqrOrthoManager_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines