Anasazi Version of the Day
AnasaziTsqrOrthoManager.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2010) 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 
00031 
00032 #ifndef __AnasaziTsqrOrthoManager_hpp
00033 #define __AnasaziTsqrOrthoManager_hpp
00034 
00035 #include "AnasaziTsqrOrthoManagerImpl.hpp"
00036 #include "AnasaziSVQBOrthoManager.hpp"
00037 
00038 
00039 namespace Anasazi {
00040 
00062   template<class Scalar, class MV>
00063   class OutOfPlaceNormalizerMixin {
00064   public:
00065     typedef Scalar scalar_type;
00066     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00069     typedef MV multivector_type;
00070 
00071     typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00072     typedef Teuchos::RCP<mat_type>                  mat_ptr;
00073 
00082     virtual int 
00083     normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const = 0;
00084 
00103     virtual int 
00104     projectAndNormalizeOutOfPlace (MV& X_in, 
00105            MV& X_out, 
00106            Teuchos::Array<mat_ptr> C,
00107            mat_ptr B,
00108            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const = 0;
00109 
00111     virtual ~OutOfPlaceNormalizerMixin () {}
00112   };
00113 
00120   template<class Scalar, class MV>
00121   class TsqrOrthoManager : 
00122     public OrthoManager<Scalar, MV>, 
00123     public OutOfPlaceNormalizerMixin<Scalar, MV>,
00124     public Teuchos::ParameterListAcceptor
00125   {
00126   public:
00127     typedef Scalar scalar_type;
00128     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00130     typedef MV multivector_type;
00131 
00132     typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00133     typedef Teuchos::RCP<mat_type>                  mat_ptr;
00134 
00135     void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
00136       impl_.setParameterList (params);
00137     }
00138 
00139     Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList () {
00140       return impl_.getNonconstParameterList ();
00141     }
00142 
00143     Teuchos::RCP<Teuchos::ParameterList> unsetParameterList () {
00144       return impl_.unsetParameterList ();
00145     }
00146 
00154     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
00155       return impl_.getValidParameters();
00156     }
00157 
00167     Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() const {
00168       return impl_.getFastParameters();
00169     }
00170 
00187     TsqrOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params, 
00188           const std::string& label = "Anasazi") :
00189       impl_ (params, label)
00190     {}
00191 
00196     TsqrOrthoManager (const std::string& label) :
00197       impl_ (label)
00198     {}
00199 
00201     virtual ~TsqrOrthoManager() {}
00202 
00203     void innerProd (const MV &X, const MV& Y, mat_type& Z) const {
00204       return impl_.innerProd (X, Y, Z);
00205     }
00206 
00207     void norm (const MV& X, std::vector<magnitude_type>& normVec) const {
00208       return impl_.norm (X, normVec);
00209     }
00210 
00211     void 
00212     project (MV &X, 
00213        Teuchos::Array<Teuchos::RCP<const MV> > Q,
00214        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > > C 
00215        = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,Scalar> >(Teuchos::null))) const
00216     {
00217       return impl_.project (X, C, Q);
00218     }
00219 
00220     int 
00221     normalize (MV &X, mat_ptr B = Teuchos::null) const
00222     {
00223       return impl_.normalize (X, B);
00224     }
00225 
00226     int 
00227     projectAndNormalize (MV &X, 
00228        Teuchos::Array<Teuchos::RCP<const MV> > Q,
00229        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > > C 
00230        = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,Scalar> >(Teuchos::null)),
00231        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > B = Teuchos::null) const
00232     {
00233       return impl_.projectAndNormalize (X, C, B, Q);
00234     }
00235 
00252     int 
00253     normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
00254     {
00255       return impl_.normalizeOutOfPlace (X, Q, B);
00256     }
00257 
00278     int 
00279     projectAndNormalizeOutOfPlace (MV& X_in, 
00280            MV& X_out,
00281            Teuchos::Array<mat_ptr> C,
00282            mat_ptr B,
00283            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00284     {
00285       return impl_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
00286     }
00287 
00288     magnitude_type orthonormError (const MV &X) const {
00289       return impl_.orthonormError (X);
00290     }
00291 
00292     magnitude_type orthogError (const MV &X1, const MV &X2) const {
00293       return impl_.orthogError (X1, X2);
00294     }
00295 
00296   private:
00302     mutable TsqrOrthoManagerImpl<Scalar, MV> impl_;
00303   };
00304 
00305 
00318   template<class Scalar, class MV, class OP>
00319   class TsqrMatOrthoManager : 
00320     public MatOrthoManager<Scalar, MV, OP>,
00321     public OutOfPlaceNormalizerMixin<Scalar, MV>,
00322     public Teuchos::ParameterListAcceptorDefaultBase
00323   {
00324   public:
00325     typedef Scalar scalar_type;
00326     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00328     typedef MV multivector_type;
00330     typedef OP operator_type;
00331 
00332     typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00333     typedef Teuchos::RCP<mat_type>                  mat_ptr;
00334 
00335   private:
00344     typedef MatOrthoManager<Scalar, MV, OP> base_type;
00345 
00348     typedef TsqrOrthoManagerImpl<Scalar, MV> tsqr_type;
00349 
00352     typedef SVQBOrthoManager<Scalar, MV, OP> svqb_type;
00353 
00356     typedef MultiVecTraits<Scalar, MV> MVT;
00357 
00358   public:
00379     TsqrMatOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params,
00380        const std::string& label = "Belos",
00381        Teuchos::RCP<const OP> Op = Teuchos::null) :
00382       MatOrthoManager<Scalar, MV, OP> (Op),
00383       tsqr_ (params, label),
00384       pSvqb_ (Teuchos::null) // Initialized lazily
00385     {}
00386 
00395     TsqrMatOrthoManager (const std::string& label = "Belos",
00396        Teuchos::RCP<const OP> Op = Teuchos::null) :
00397       MatOrthoManager<Scalar, MV, OP> (Op),
00398       tsqr_ (label),
00399       pSvqb_ (Teuchos::null) // Initialized lazily
00400     {}
00401 
00403     virtual ~TsqrMatOrthoManager() {}
00404 
00412     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
00413       return tsqr_.getValidParameters ();
00414     }
00415 
00425     Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() {
00426       return tsqr_.getFastParameters ();
00427     }
00428 
00429     void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
00430       tsqr_.setParameterList (params);
00431     }
00432 
00433     virtual void 
00434     setOp (Teuchos::RCP< const OP > Op) 
00435     {
00436       // We override the base class' setOp() so that the
00437       // SVQBOrthoManager instance gets the new op.
00438       //
00439       // The base_type notation helps C++ resolve the name for a
00440       // member function of a templated base class.
00441       base_type::setOp (Op); // base class gets a copy of the Op too
00442 
00443       if (! Op.is_null()) {
00444   ensureSvqbInit (); // Make sure the SVQB object has been initialized
00445   pSvqb_->setOp (Op);
00446       }
00447     }
00448 
00449     Teuchos::RCP<const OP> getOp () const { 
00450       // The base_type notation helps C++ resolve the name for a
00451       // member function of a templated base class.
00452       return base_type::getOp(); 
00453     }
00454 
00455     void 
00456     projectMat (MV &X, 
00457     Teuchos::Array<Teuchos::RCP<const MV> > Q,
00458     Teuchos::Array<Teuchos::RCP<mat_type> > C = 
00459       Teuchos::tuple (Teuchos::RCP<mat_type> (Teuchos::null)),
00460     Teuchos::RCP<MV> MX = Teuchos::null,
00461     Teuchos::Array<Teuchos::RCP<const MV> > MQ = 
00462       Teuchos::tuple (Teuchos::null)) const
00463     {
00464       if (getOp().is_null()) {
00465   // FIXME (mfh 12 Jan 2011): Do we need to check if C is null
00466   // and allocate space if so?
00467   tsqr_.project (X, C, Q);
00468   // FIXME (mfh 20 Jul 2010) What about MX and MQ?
00469   // 
00470   // FIXME (mfh 12 Jan 2011) Need to port MVT::Assign from Belos to Anasazi
00471 #if 0
00472   if (! MX.is_null()) {
00473     // MX gets a copy of X; M is the identity operator.
00474     MVT::Assign (X, *MX);
00475   }
00476 #endif // 0
00477       } else {
00478   ensureSvqbInit ();
00479   pSvqb_->projectMat (X, Q, C, MX, MQ);
00480       }
00481     }
00482 
00483     int 
00484     normalizeMat (MV &X, 
00485       mat_ptr B = Teuchos::null,
00486       Teuchos::RCP<MV> MX = Teuchos::null) const
00487     {
00488       if (getOp().is_null()) {
00489   // FIXME (mfh 12 Jan 2011): Do we need to check if B is null
00490   // and allocate space if so?
00491   const int rank = tsqr_.normalize (X, B);
00492   // FIXME (mfh 20 Jul 2010) What about MX?
00493   // 
00494   // FIXME (mfh 12 Jan 2011) Need to port MVT::Assign from Belos to Anasazi
00495 #if 0
00496   if (! MX.is_null()) {
00497     // MX gets a copy of X; M is the identity operator.
00498     MVT::Assign (X, *MX);
00499   }
00500 #endif // 0
00501   return rank;
00502       } else {
00503   ensureSvqbInit ();
00504   return pSvqb_->normalizeMat (X, B, MX);
00505       }
00506     }
00507 
00508     int 
00509     projectAndNormalizeMat (MV &X, 
00510           Teuchos::Array<Teuchos::RCP<const MV> > Q,
00511           Teuchos::Array<Teuchos::RCP<mat_type> > C = 
00512             Teuchos::tuple (Teuchos::RCP<mat_type> (Teuchos::null)),
00513           Teuchos::RCP<mat_type> B = Teuchos::null,
00514           Teuchos::RCP<MV> MX = Teuchos::null,
00515           Teuchos::Array<Teuchos::RCP<const MV> > MQ = 
00516             Teuchos::tuple (Teuchos::RCP<const MV> (Teuchos::null))) const 
00517     {
00518       if (getOp().is_null()) {
00519   // FIXME (mfh 12 Jan 2011): Do we need to check if C and B
00520   // are null and allocate space if so?
00521   const int rank = tsqr_.projectAndNormalize (X, C, B, Q); 
00522   // FIXME (mfh 20 Jul 2010) What about MX and MQ?
00523   // mfh 12 Jan 2011: Ignore MQ (assume Q == MQ).
00524   // 
00525   // FIXME (mfh 12 Jan 2011) Need to port MVT::Assign from Belos to Anasazi
00526 #if 0
00527   if (! MX.is_null()) {
00528     // MX gets a copy of X; M is the identity operator.
00529     MVT::Assign (X, *MX);
00530   }
00531 #endif // 0
00532   return rank;
00533       } else {
00534   ensureSvqbInit ();
00535   return pSvqb_->projectAndNormalizeMat (X, Q, C, B, MX, MQ);
00536       }
00537     }
00538 
00539     int 
00540     normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
00541     {
00542       if (getOp().is_null()) {
00543   return tsqr_.normalizeOutOfPlace (X, Q, B);
00544       } else {
00545   // SVQB normalizes in place, so we have to copy.
00546   ensureSvqbInit ();
00547   const int rank = pSvqb_->normalize (X, B);
00548   MVT::Assign (X, Q);
00549   return rank;
00550       }
00551     }
00552 
00553     int 
00554     projectAndNormalizeOutOfPlace (MV& X_in, 
00555            MV& X_out,
00556            Teuchos::Array<mat_ptr> C,
00557            mat_ptr B,
00558            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00559     {
00560       using Teuchos::null;
00561 
00562       if (getOp().is_null()) {
00563   return tsqr_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
00564       } else {
00565   ensureSvqbInit ();
00566   // SVQB's projectAndNormalize wants an Array, not an ArrayView.
00567   // Copy into a temporary Array and copy back afterwards.
00568   Teuchos::Array<Teuchos::RCP<const MV> > Q_array (Q);
00569   const int rank = pSvqb_->projectAndNormalize (X_in, Q_array, C, B);
00570   Q.assign (Q_array);
00571   // SVQB normalizes in place, so we have to copy X_in to X_out.
00572   MVT::Assign (X_in, X_out);
00573   return rank;
00574       }
00575     }
00576 
00577     magnitude_type
00578     orthonormErrorMat (const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const
00579     {
00580       if (getOp().is_null()) {
00581   return tsqr_.orthonormError (X);
00582   // FIXME (mfh 20 Jul 2010) What about MX?
00583       } else {
00584   ensureSvqbInit ();
00585   return pSvqb_->orthonormErrorMat (X, MX);
00586       }
00587     }
00588 
00589     magnitude_type
00590     orthogErrorMat (const MV &X, 
00591         const MV &Y,
00592         Teuchos::RCP<const MV> MX = Teuchos::null, 
00593         Teuchos::RCP<const MV> MY = Teuchos::null) const
00594     {
00595       if (getOp().is_null()) {
00596   return tsqr_.orthogError (X, Y);
00597   // FIXME (mfh 20 Jul 2010) What about MX and MY?
00598       } else {
00599   ensureSvqbInit ();
00600   return pSvqb_->orthogErrorMat (X, Y, MX, MY);
00601       }
00602     }
00603 
00604   private:
00606     void 
00607     ensureSvqbInit () const
00608     {
00609       // NOTE (mfh 12 Oct 2011) For now, we rely on the default values
00610       // of SVQB parameters.
00611       if (pSvqb_.is_null()) {
00612   pSvqb_ = Teuchos::rcp (new svqb_type (getOp()));
00613       }
00614     }
00615 
00623     mutable tsqr_type tsqr_;
00624 
00629     mutable Teuchos::RCP<svqb_type> pSvqb_;
00630   };
00631 
00632 } // namespace Anasazi
00633 
00634 #endif // __AnasaziTsqrOrthoManager_hpp
00635 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends