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     // Fix warning about hiding OrthoManager::projectAndNormalize()
00562     using Belos::OrthoManager<Scalar, MV>::projectAndNormalize;
00563 
00564     int 
00565     projectAndNormalize (MV &X, 
00566        Teuchos::RCP<MV> MX,
00567        Teuchos::Array<mat_ptr> C,
00568        mat_ptr B, 
00569        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00570     {
00571       if (getOp().is_null()) {
00572   const int rank = tsqr_.projectAndNormalize (X, C, B, Q); 
00573   if (! MX.is_null()) {
00574     // MX gets a copy of X; M is the identity operator.
00575     MVT::Assign (X, *MX);
00576   }
00577   return rank;
00578       } else {
00579   ensureDgksInit ();
00580   return pDgks_->projectAndNormalize (X, MX, C, B, Q);
00581       }
00582     }
00583 
00584     int 
00585     normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
00586     {
00587       if (getOp().is_null()) {
00588   return tsqr_.normalizeOutOfPlace (X, Q, B);
00589       } else {
00590   // DGKS normalizes in place, so we have to copy.
00591   ensureDgksInit ();
00592   const int rank = pDgks_->normalize (X, B);
00593   MVT::Assign (X, Q);
00594   return rank;
00595       }
00596     }
00597 
00598     int 
00599     projectAndNormalizeOutOfPlace (MV& X_in, 
00600            MV& X_out,
00601            Teuchos::Array<mat_ptr> C,
00602            mat_ptr B,
00603            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00604     {
00605       using Teuchos::null;
00606 
00607       if (getOp().is_null()) {
00608   return tsqr_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
00609       } else {
00610   // DGKS normalizes in place, so we have to copy.
00611   ensureDgksInit ();
00612   const int rank = pDgks_->projectAndNormalize (X_in, null, C, B, Q);
00613   MVT::Assign (X_in, X_out);
00614   return rank;
00615       }
00616     }
00617 
00618     magnitude_type 
00619     orthonormError (const MV &X, Teuchos::RCP<const MV> MX) const
00620     {
00621       if (getOp().is_null()) {
00622   return tsqr_.orthonormError (X); // Ignore MX
00623       } else {
00624   ensureDgksInit ();
00625   return pDgks_->orthonormError (X, MX);
00626       }
00627     }
00628 
00629     magnitude_type orthonormError (const MV &X) const {
00630       return orthonormError (X, Teuchos::null);
00631     }
00632 
00633     magnitude_type orthogError (const MV &X1, const MV &X2) const {
00634       return orthogError (X1, Teuchos::null, X2);
00635     }
00636 
00637     magnitude_type
00638     orthogError (const MV &X1, 
00639      Teuchos::RCP<const MV> MX1,
00640      const MV &X2) const
00641     {
00642       if (getOp().is_null()) {
00643   // Ignore MX1, since we don't need to write to it.
00644   return tsqr_.orthogError (X1, X2);
00645       } else {
00646   ensureDgksInit ();
00647   return pDgks_->orthogError (X1, MX1, X2);
00648       }
00649     }
00650 
00651     void 
00652     setLabel (const std::string& label) 
00653     {
00654       tsqr_.setLabel (label);
00655 
00656       // Make sure DGKS gets the new label, if it's initialized.
00657       // Otherwise, it will get the new label on initialization.
00658       if (! pDgks_.is_null()) {
00659   pDgks_->setLabel (label);
00660       }
00661     }
00662 
00663   private:
00665     void 
00666     ensureDgksInit () const
00667     {
00668       // NOTE (mfh 11 Jan 2011) DGKS has a parameter that needs to be
00669       // set.  For now, we just use the default value of the
00670       // parameter.
00671       if (pDgks_.is_null()) {
00672   pDgks_ = Teuchos::rcp (new dgks_type (getLabel(), getOp()));
00673       }
00674     }
00675 
00683     mutable tsqr_type tsqr_;
00684 
00689     mutable Teuchos::RCP<dgks_type> pDgks_;
00690   };
00691 
00692 } // namespace Belos
00693 
00694 #endif // __BelosTsqrOrthoManager_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines