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 
00056 
00057 namespace Belos {
00058 
00076   template<class Scalar, class MV>
00077   class OutOfPlaceNormalizerMixin {
00078   public:
00079     typedef Scalar scalar_type;
00080     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00082     typedef MV multivector_type;
00083 
00084     typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00085     typedef Teuchos::RCP<mat_type>                  mat_ptr;
00086 
00097     virtual int 
00098     normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const = 0;
00099 
00118     virtual int 
00119     projectAndNormalizeOutOfPlace (MV& X_in, 
00120            MV& X_out, 
00121            Teuchos::Array<mat_ptr> C,
00122            mat_ptr B,
00123            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const = 0;
00124 
00126     virtual ~OutOfPlaceNormalizerMixin () {}
00127   };
00128 
00134   template<class Scalar, class MV>
00135   class TsqrOrthoManager : 
00136     public OrthoManager<Scalar, MV>, 
00137     public OutOfPlaceNormalizerMixin<Scalar, MV> {
00138   public:
00139     typedef Scalar scalar_type;
00140     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00142     typedef MV multivector_type;
00143 
00144     typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00145     typedef Teuchos::RCP<mat_type>                  mat_ptr;
00146 
00163     static Teuchos::RCP<const Teuchos::ParameterList> getDefaultParameters() {
00164       return TsqrOrthoManagerImpl<Scalar, MV>::getDefaultParameters();
00165     }
00166 
00181     static Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() {
00182       return TsqrOrthoManagerImpl<Scalar, MV>::getFastParameters();
00183     }
00184 
00196     TsqrOrthoManager (const Teuchos::RCP<const Teuchos::ParameterList>& params, 
00197           const std::string& label = "Belos") :
00198       impl_ (params, label)
00199     {}
00200 
00202     virtual ~TsqrOrthoManager() {}
00203 
00205     void innerProd (const MV &X, const MV &Y, mat_type& Z) const {
00206       return impl_.innerProd (X, Y, Z);
00207     }
00208 
00216     void norm (const MV& X, std::vector<magnitude_type>& normVec) const {
00217       return impl_.norm (X, normVec);
00218     }
00219 
00221     void 
00222     project (MV &X, 
00223        Teuchos::Array<mat_ptr> C,
00224        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00225     {
00226       return impl_.project (X, C, Q);
00227     }
00228 
00232     int 
00233     normalize (MV &X, mat_ptr B) const
00234     {
00235       return impl_.normalize (X, B);
00236     }
00237 
00248     int 
00249     projectAndNormalize (MV &X, 
00250        Teuchos::Array<mat_ptr> C,
00251        mat_ptr B,
00252        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00253     {
00254       return impl_.projectAndNormalize (X, C, B, Q);
00255     }
00256 
00275     int 
00276     normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
00277     {
00278       return impl_.normalizeOutOfPlace (X, Q, B);
00279     }
00280 
00300     int 
00301     projectAndNormalizeOutOfPlace (MV& X_in, 
00302            MV& X_out,
00303            Teuchos::Array<mat_ptr> C,
00304            mat_ptr B,
00305            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00306     {
00307       return impl_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
00308     }
00309 
00311     magnitude_type orthonormError (const MV &X) const {
00312       return impl_.orthonormError (X);
00313     }
00314 
00315     magnitude_type orthogError (const MV &X1, const MV &X2) const {
00316       return impl_.orthogError (X1, X2);
00317     }
00318 
00326     void setLabel (const std::string& label) { impl_.setLabel (label); }
00327 
00329     const std::string& getLabel() const { return impl_.getLabel(); }
00330 
00331   private:
00335     mutable TsqrOrthoManagerImpl<Scalar, MV> impl_;
00336 
00338     std::string label_;
00339   };
00340 
00341 
00353   template<class Scalar, class MV, class OP>
00354   class TsqrMatOrthoManager : 
00355     public MatOrthoManager<Scalar, MV, OP>,
00356     public OutOfPlaceNormalizerMixin<Scalar, MV>
00357   {
00358   public:
00359     typedef Scalar scalar_type;
00360     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00362     typedef MV multivector_type;
00364     typedef OP operator_type;
00365 
00366     typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00367     typedef Teuchos::RCP<mat_type>                  mat_ptr;
00368 
00369   private:
00378     typedef MatOrthoManager<Scalar, MV, OP> base_type;
00379 
00382     typedef TsqrOrthoManagerImpl<Scalar, MV> tsqr_type;
00383 
00386     typedef DGKSOrthoManager<Scalar, MV, OP> dgks_type;
00387 
00390     typedef MultiVecTraits<Scalar, MV> MVT;
00391 
00392   public:
00409     static Teuchos::RCP<const Teuchos::ParameterList> getDefaultParameters() {
00410       return TsqrOrthoManagerImpl<Scalar, MV>::getDefaultParameters();
00411     }
00412 
00427     static Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() {
00428       return TsqrOrthoManagerImpl<Scalar, MV>::getFastParameters();
00429     }
00430 
00432     TsqrMatOrthoManager () :
00433       MatOrthoManager<Scalar, MV, OP>(Teuchos::null),
00434       pTsqr_ (Teuchos::null), // Lazy initialization
00435       pDgks_ (Teuchos::null)  // Lazy initialization
00436     {}
00437 
00442     void setLabel (const std::string& label) {
00443       label_ = label; 
00444     }
00446     const std::string& getLabel() const { return label_; }
00447 
00459     TsqrMatOrthoManager (const Teuchos::RCP<const Teuchos::ParameterList>& params,
00460        const std::string& label = "Belos",
00461        Teuchos::RCP< const OP > Op = Teuchos::null) :
00462       MatOrthoManager<Scalar, MV, OP>(Op),
00463       params_ (params),
00464       label_ (label),
00465       pTsqr_ (Teuchos::null), // Lazy initialization
00466       pDgks_ (Teuchos::null)  // Lazy initialization
00467     {}
00468 
00470     virtual ~TsqrMatOrthoManager() {}
00471 
00480     void 
00481     setOp (Teuchos::RCP<const OP> Op) 
00482     {
00483       // We use this notation to help C++ resolve the name.
00484       // Otherwise, it won't know where to find setOp(), since this is
00485       // a member function of the base class which does not depend on
00486       // the template parameters.
00487       base_type::setOp (Op); // base class gets a copy of the Op too
00488       ensureDgksInit (); // Make sure the DGKS object has been initialized
00489       pDgks_->setOp (Op);
00490     }
00491 
00496     Teuchos::RCP<const OP> getOp () const { 
00497       return base_type::getOp(); 
00498     }
00499 
00505     void 
00506     project (MV &X, 
00507        Teuchos::RCP<MV> MX,
00508        Teuchos::Array<mat_ptr> C,
00509        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00510     {
00511       if (getOp().is_null())
00512   {
00513     ensureTsqrInit ();
00514     pTsqr_->project (X, C, Q);
00515     if (! MX.is_null())
00516       // MX gets a copy of X; M is the identity operator.
00517       MVT::Assign (X, *MX);
00518   }
00519       else
00520   {
00521     ensureDgksInit ();
00522     pDgks_->project (X, MX, C, Q);
00523   }
00524     }
00525 
00528     void 
00529     project (MV &X, 
00530        Teuchos::Array<mat_ptr> C,
00531        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00532       project (X, Teuchos::null, C, Q);
00533     }
00534 
00535     int 
00536     normalize (MV& X, Teuchos::RCP<MV> MX, mat_ptr B) const 
00537     {
00538       if (getOp().is_null())
00539   {
00540     ensureTsqrInit ();
00541     const int rank = pTsqr_->normalize (X, B);
00542     if (! MX.is_null())
00543       // MX gets a copy of X; M is the identity operator.
00544       MVT::Assign (X, *MX);
00545     return rank;
00546   }
00547       else
00548   {
00549     ensureDgksInit ();
00550     return pDgks_->normalize (X, MX, B);
00551   }
00552     }
00553 
00554     int normalize (MV& X, mat_ptr B) const {
00555       return normalize (X, Teuchos::null, B);
00556     }
00557 
00558     int 
00559     projectAndNormalize (MV &X, 
00560        Teuchos::RCP<MV> MX,
00561        Teuchos::Array<mat_ptr> C,
00562        mat_ptr B, 
00563        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00564     {
00565       if (getOp().is_null())
00566   {
00567     ensureTsqrInit ();
00568     const int rank = pTsqr_->projectAndNormalize (X, C, B, Q); 
00569     if (! MX.is_null())
00570       // MX gets a copy of X; M is the identity operator.
00571       MVT::Assign (X, *MX);
00572     return rank;
00573   }
00574       else
00575   {
00576     ensureDgksInit ();
00577     return pDgks_->projectAndNormalize (X, MX, C, B, Q);
00578   }
00579     }
00580 
00599     int 
00600     normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
00601     {
00602       if (getOp().is_null())
00603   {
00604     ensureTsqrInit ();
00605     return pTsqr_->normalizeOutOfPlace (X, Q, B);
00606   }
00607       else
00608   {
00609     // DGKS normalizes in place, so we have to copy.
00610     ensureDgksInit ();
00611     const int rank = pDgks_->normalize (X, B);
00612     MVT::Assign (X, Q);
00613     return rank;
00614   }
00615     }
00616 
00636     int 
00637     projectAndNormalizeOutOfPlace (MV& X_in, 
00638            MV& X_out,
00639            Teuchos::Array<mat_ptr> C,
00640            mat_ptr B,
00641            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00642     {
00643       if (getOp().is_null())
00644   {
00645     ensureTsqrInit ();
00646     return pTsqr_->projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
00647   }
00648       else
00649   {
00650     // DGKS normalizes in place, so we have to copy.
00651     ensureDgksInit ();
00652     const int rank = pDgks_->projectAndNormalize (X_in, Teuchos::null, C, B, Q);
00653     MVT::Assign (X_in, X_out);
00654     return rank;
00655   }
00656     }
00657 
00658     magnitude_type 
00659     orthonormError (const MV &X, Teuchos::RCP<const MV> MX) const
00660     {
00661       if (getOp().is_null())
00662   {
00663     ensureTsqrInit ();
00664     return pTsqr_->orthonormError (X); // Ignore MX
00665   }
00666       else
00667   {
00668     ensureDgksInit ();
00669     return pDgks_->orthonormError (X, MX);
00670   }
00671     }
00672 
00673     magnitude_type orthonormError (const MV &X) const {
00674       return orthonormError (X, Teuchos::null);
00675     }
00676 
00677     magnitude_type orthogError (const MV &X1, const MV &X2) const {
00678       return orthogError (X1, Teuchos::null, X2);
00679     }
00680 
00681     magnitude_type
00682     orthogError (const MV &X1, 
00683      Teuchos::RCP<const MV> MX1,
00684      const MV &X2) const
00685     {
00686       if (getOp().is_null())
00687   {
00688     ensureTsqrInit ();
00689     // Ignore MX1, since we don't need to write to it.
00690     return pTsqr_->orthogError (X1, X2);
00691   }
00692       else
00693   {
00694     ensureDgksInit ();
00695     return pDgks_->orthogError (X1, MX1, X2);
00696   }
00697     }
00698 
00699   private:
00700     void
00701     ensureTsqrInit () const
00702     {
00703       if (pTsqr_.is_null())
00704   pTsqr_ = Teuchos::rcp (new tsqr_type (params_, getLabel()));
00705     }
00706     void 
00707     ensureDgksInit () const
00708     {
00709       // NOTE (mfh 11 Jan 2011) DGKS has a parameter that needs to be
00710       // set.  For now, we just use the default value of the
00711       // parameter.
00712       if (pDgks_.is_null())
00713   pDgks_ = Teuchos::rcp (new dgks_type (getLabel(), getOp()));
00714     }
00715 
00717     Teuchos::RCP<const Teuchos::ParameterList> params_;
00719     std::string label_;
00722     mutable Teuchos::RCP<tsqr_type> pTsqr_;
00725     mutable Teuchos::RCP<dgks_type> pDgks_;
00726   };
00727 
00728 } // namespace Belos
00729 
00730 #endif // __BelosTsqrOrthoManager_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines