Belos Package Browser (Single Doxygen Collection) Development
BelosICGSOrthoManager.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 
00042 
00047 #ifndef BELOS_ICGS_ORTHOMANAGER_HPP
00048 #define BELOS_ICGS_ORTHOMANAGER_HPP
00049 
00057 // #define ORTHO_DEBUG
00058 
00059 #include "BelosConfigDefs.hpp"
00060 #include "BelosMultiVecTraits.hpp"
00061 #include "BelosOperatorTraits.hpp"
00062 #include "BelosMatOrthoManager.hpp"
00063 
00064 #include "Teuchos_as.hpp"
00065 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
00066 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00067 #include "Teuchos_TimeMonitor.hpp"
00068 #endif // BELOS_TEUCHOS_TIME_MONITOR
00069 
00070 namespace Belos {
00071 
00076   template<class ScalarType>
00077   Teuchos::RCP<const Teuchos::ParameterList> TEUCHOS_DEPRECATED
00078   getDefaultIcgsParameters()
00079   {
00080     using Teuchos::as;
00081     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00082     typedef Teuchos::ScalarTraits<magnitude_type> STM;
00083 
00084     // Default parameter values for ICGS orthogonalization.
00085     // Documentation will be embedded in the parameter list.
00086     const int defaultMaxNumOrthogPasses = 2;
00087     const magnitude_type eps = STM::eps();
00088     const magnitude_type defaultBlkTol = 
00089       as<magnitude_type> (10) * STM::squareroot (eps);
00090     const magnitude_type defaultSingTol = as<magnitude_type> (10) * eps;
00091 
00092     Teuchos::RCP<Teuchos::ParameterList> params = 
00093       Teuchos::parameterList ("ICGS");
00094     params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
00095                  "Maximum number of orthogonalization passes "
00096                  "(includes the first).  Default is 2, since "
00097                  "\"twice is enough\" for Krylov methods.");
00098     params->set ("blkTol", defaultBlkTol, 
00099                  "Block reorthogonalization threshhold.");
00100     params->set ("singTol", defaultSingTol, 
00101                  "Singular block detection threshold.");
00102     return params;
00103   }
00104 
00109   template<class ScalarType>
00110   Teuchos::RCP<const Teuchos::ParameterList> TEUCHOS_DEPRECATED
00111   getFastIcgsParameters()
00112   {
00113     using Teuchos::ParameterList;
00114     using Teuchos::RCP;
00115     using Teuchos::rcp;
00116     using Teuchos::ScalarTraits;
00117     typedef typename ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00118     typedef ScalarTraits<magnitude_type> STM;
00119 
00120     RCP<const ParameterList> defaultParams = getDefaultIcgsParameters<ScalarType>();
00121     // Start with a clone of the default parameters.
00122     RCP<ParameterList> params = rcp (new ParameterList (*defaultParams));
00123 
00124     const int maxBlkOrtho = 1;
00125     const magnitude_type blkTol = STM::zero();
00126     const magnitude_type singTol = STM::zero();
00127 
00128     params->set ("maxNumOrthogPasses", maxBlkOrtho);
00129     params->set ("blkTol", blkTol);
00130     params->set ("singTol", singTol);
00131 
00132     return params;
00133   }
00134 
00144   template<class ScalarType>
00145   void TEUCHOS_DEPRECATED
00146   readIcgsParameters (const Teuchos::RCP<const Teuchos::ParameterList>& params,
00147                       int& maxNumOrthogPasses,
00148                       typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& blkTol,
00149                       typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& singTol)
00150   {
00151     using Teuchos::ParameterList;
00152     using Teuchos::RCP;
00153     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00154     const magnitude_type zero = Teuchos::ScalarTraits<magnitude_type>::zero();
00155     RCP<const ParameterList> defaultParams = getDefaultIcgsParameters<ScalarType>();
00156 
00157     // Using temporary variables and fetching all values before
00158     // setting the output arguments ensures the strong exception
00159     // guarantee for this function: if an exception is thrown, no
00160     // externally visible side effects (in this case, setting the
00161     // output arguments) have taken place.
00162     int _maxNumOrthogPasses;
00163     magnitude_type _blkTol, _singTol;
00164     if (params.is_null()) {
00165       _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00166       _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00167       _singTol = defaultParams->get<magnitude_type> ("singTol");
00168     } else {
00169       try {
00170         _maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
00171         if (_maxNumOrthogPasses < 1)
00172           _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00173       } catch (Teuchos::Exceptions::InvalidParameter&) {
00174         _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00175       }
00176 
00177       try {
00178   _blkTol = params->get<magnitude_type> ("blkTol");
00179   if (_blkTol < zero)
00180     _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00181       } catch (Teuchos::Exceptions::InvalidParameter&) {
00182   try {
00183     // People may have used depTol instead of blkTol for this
00184     // parameter's name, by analogy with DGKS.
00185     _blkTol = params->get<magnitude_type> ("depTol");
00186   } catch (Teuchos::Exceptions::InvalidParameter&) {
00187     _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00188   }
00189       }
00190 
00191       try {
00192         _singTol = params->get<magnitude_type> ("singTol");
00193         if (_singTol < zero)
00194           _singTol = defaultParams->get<magnitude_type> ("singTol");
00195       } catch (Teuchos::Exceptions::InvalidParameter&) {
00196         _singTol = defaultParams->get<magnitude_type> ("singTol");
00197       }
00198     }
00199     maxNumOrthogPasses = _maxNumOrthogPasses;
00200     blkTol = _blkTol;
00201     singTol = _singTol;
00202   }
00203 
00204 
00205   template<class ScalarType, class MV, class OP>
00206   class ICGSOrthoManager : 
00207     public MatOrthoManager<ScalarType,MV,OP>, 
00208     public Teuchos::ParameterListAcceptorDefaultBase
00209   {
00210   private:
00211     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00212     typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
00213     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
00214     typedef MultiVecTraits<ScalarType,MV>      MVT;
00215     typedef OperatorTraits<ScalarType,MV,OP>   OPT;
00216 
00217   public:
00219 
00220 
00222     ICGSOrthoManager( const std::string& label = "Belos",
00223                       Teuchos::RCP<const OP> Op = Teuchos::null,
00224                       const int max_ortho_steps = 2,
00225                       const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
00226                       const MagnitudeType sing_tol = 10*MGT::eps() )
00227       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00228       max_ortho_steps_( max_ortho_steps ),
00229       blk_tol_( blk_tol ),
00230       sing_tol_( sing_tol ),
00231       label_( label )
00232     {
00233 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00234       std::string orthoLabel = label_ + ": Orthogonalization";
00235       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00236 
00237       std::string updateLabel = label_ + ": Ortho (Update)";
00238       timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00239 
00240       std::string normLabel = label_ + ": Ortho (Norm)";
00241       timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00242 
00243       std::string ipLabel = label_ + ": Ortho (Inner Product)";
00244       timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel); 
00245 #endif
00246     }
00247 
00249     ICGSOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& plist,
00250           const std::string& label = "Belos",
00251                       Teuchos::RCP<const OP> Op = Teuchos::null) :
00252       MatOrthoManager<ScalarType,MV,OP>(Op), 
00253       max_ortho_steps_ (2),
00254       blk_tol_ (10 * MGT::squareroot (MGT::eps())),
00255       sing_tol_ (10 * MGT::eps()),
00256       label_ (label)
00257     {
00258       setParameterList (plist);
00259 
00260 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00261       std::string orthoLabel = label_ + ": Orthogonalization";
00262       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00263 
00264       std::string updateLabel = label_ + ": Ortho (Update)";
00265       timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00266 
00267       std::string normLabel = label_ + ": Ortho (Norm)";
00268       timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00269 
00270       std::string ipLabel = label_ + ": Ortho (Inner Product)";
00271       timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel); 
00272 #endif
00273     }
00274 
00276     ~ICGSOrthoManager() {}
00278 
00280 
00281 
00282     void 
00283     setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
00284     {
00285       using Teuchos::Exceptions::InvalidParameterName;
00286       using Teuchos::ParameterList;
00287       using Teuchos::parameterList;
00288       using Teuchos::RCP;
00289 
00290       RCP<const ParameterList> defaultParams = getValidParameters();
00291       RCP<ParameterList> params;
00292       if (plist.is_null()) {
00293   params = parameterList (*defaultParams);
00294       } else {
00295   params = plist;
00296   // Some users might want to specify "blkTol" as "depTol".  Due
00297   // to this case, we don't invoke
00298   // validateParametersAndSetDefaults on params.  Instead, we go
00299   // through the parameter list one parameter at a time and look
00300   // for alternatives.
00301       }
00302   
00303       // Using temporary variables and fetching all values before
00304       // setting the output arguments ensures the strong exception
00305       // guarantee for this function: if an exception is thrown, no
00306       // externally visible side effects (in this case, setting the
00307       // output arguments) have taken place.
00308       int maxNumOrthogPasses;
00309       MagnitudeType blkTol;
00310       MagnitudeType singTol;
00311 
00312       try {
00313   maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
00314       } catch (InvalidParameterName&) {
00315   maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00316   params->set ("maxNumOrthogPasses", maxNumOrthogPasses);
00317       }
00318 
00319       // Handling of the "blkTol" parameter is a special case.  This
00320       // is because some users may prefer to call this parameter
00321       // "depTol" for consistency with DGKS.  However, our default
00322       // parameter list calls this "blkTol", and we don't want the
00323       // default list's value to override the user's value.  Thus, we
00324       // first check the user's parameter list for both names, and
00325       // only then access the default parameter list.
00326       try {
00327   blkTol = params->get<MagnitudeType> ("blkTol");
00328       } catch (InvalidParameterName&) {
00329   try {
00330     blkTol = params->get<MagnitudeType> ("depTol");
00331     // "depTol" is the wrong name, so remove it and replace with
00332     // "blkTol".  We'll set "blkTol" below.
00333     params->remove ("depTol");
00334   } catch (InvalidParameterName&) {
00335     blkTol = defaultParams->get<MagnitudeType> ("blkTol");
00336   }
00337   params->set ("blkTol", blkTol);
00338       }
00339 
00340       try {
00341   singTol = params->get<MagnitudeType> ("singTol");
00342       } catch (InvalidParameterName&) {
00343   singTol = defaultParams->get<MagnitudeType> ("singTol");
00344   params->set ("singTol", singTol);
00345       }
00346 
00347       max_ortho_steps_ = maxNumOrthogPasses;
00348       blk_tol_ = blkTol;
00349       sing_tol_ = singTol;
00350 
00351       setMyParamList (params);
00352     }
00353 
00354     Teuchos::RCP<const Teuchos::ParameterList> 
00355     getValidParameters () const
00356     {
00357       using Teuchos::as;
00358       using Teuchos::ParameterList;
00359       using Teuchos::parameterList;
00360       using Teuchos::RCP;
00361 
00362       if (defaultParams_.is_null()) {
00363   RCP<ParameterList> params = parameterList ("ICGS");
00364 
00365   // Default parameter values for ICGS orthogonalization.
00366   // Documentation will be embedded in the parameter list.
00367   const int defaultMaxNumOrthogPasses = 2;
00368   const MagnitudeType eps = MGT::eps();
00369   const MagnitudeType defaultBlkTol = 
00370     as<MagnitudeType> (10) * MGT::squareroot (eps);
00371   const MagnitudeType defaultSingTol = as<MagnitudeType> (10) * eps;
00372 
00373   params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
00374          "Maximum number of orthogonalization passes (includes the "
00375          "first).  Default is 2, since \"twice is enough\" for Krylov "
00376          "methods.");
00377   params->set ("blkTol", defaultBlkTol, "Block reorthogonalization "
00378          "threshhold.");
00379   params->set ("singTol", defaultSingTol, "Singular block detection "
00380          "threshold.");
00381   defaultParams_ = params;
00382       }
00383       return defaultParams_;
00384     }
00386 
00391     Teuchos::RCP<const Teuchos::ParameterList> 
00392     getFastParameters () const 
00393     {
00394       using Teuchos::as;
00395       using Teuchos::ParameterList;
00396       using Teuchos::parameterList;
00397       using Teuchos::RCP;
00398 
00399       RCP<const ParameterList> defaultParams = getValidParameters ();
00400       // Start with a clone of the default parameters.
00401       RCP<ParameterList> params = parameterList (*defaultParams);
00402 
00403       const int maxBlkOrtho = 1; // No block reorthogonalization
00404       const MagnitudeType blkTol = MGT::zero();
00405       const MagnitudeType singTol = MGT::zero();
00406 
00407       params->set ("maxNumOrthogPasses", maxBlkOrtho);
00408       params->set ("blkTol", blkTol);
00409       params->set ("singTol", singTol);
00410 
00411       return params;
00412     }
00413 
00415 
00416 
00418     void setBlkTol( const MagnitudeType blk_tol ) { 
00419       // Update the parameter list as well.
00420       Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
00421       if (! params.is_null()) {
00422   // If it's null, then we haven't called setParameterList()
00423   // yet.  It's entirely possible to construct the parameter
00424   // list on demand, so we don't try to create the parameter
00425   // list here.
00426   params->set ("blkTol", blk_tol);
00427       }
00428       blk_tol_ = blk_tol; 
00429     }
00430 
00432     void setSingTol( const MagnitudeType sing_tol ) { 
00433       // Update the parameter list as well.
00434       Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
00435       if (! params.is_null()) {
00436   // If it's null, then we haven't called setParameterList()
00437   // yet.  It's entirely possible to construct the parameter
00438   // list on demand, so we don't try to create the parameter
00439   // list here.
00440   params->set ("singTol", sing_tol);
00441       }
00442       sing_tol_ = sing_tol; 
00443     }
00444 
00446     MagnitudeType getBlkTol() const { return blk_tol_; } 
00447 
00449     MagnitudeType getSingTol() const { return sing_tol_; } 
00450 
00452 
00453 
00455 
00456 
00484     void project ( MV &X, Teuchos::RCP<MV> MX, 
00485                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00486                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00487 
00488 
00491     void project ( MV &X, 
00492                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00493                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00494       project(X,Teuchos::null,C,Q);
00495     }
00496 
00497 
00498  
00523     int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00524                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
00525 
00526 
00529     int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00530       return normalize(X,Teuchos::null,B);
00531     }
00532 
00533 
00575     int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00576                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00577                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00578                               Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00579 
00582     int projectAndNormalize ( MV &X, 
00583                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00584                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00585                               Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const {
00586       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00587     }
00588 
00590 
00592 
00593 
00597     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00598     orthonormError(const MV &X) const {
00599       return orthonormError(X,Teuchos::null);
00600     }
00601 
00606     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00607     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
00608 
00612     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00613     orthogError(const MV &X1, const MV &X2) const {
00614       return orthogError(X1,Teuchos::null,X2);
00615     }
00616 
00621     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00622     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
00623 
00625 
00627 
00628 
00631     void setLabel(const std::string& label);
00632 
00635     const std::string& getLabel() const { return label_; }
00636 
00638 
00639   private:
00641     int max_ortho_steps_;
00643     MagnitudeType blk_tol_;
00645     MagnitudeType sing_tol_;
00646 
00648     std::string label_;
00649 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00650     Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, 
00651       timerNorm_, timerScale_, timerInnerProd_;
00652 #endif // BELOS_TEUCHOS_TIME_MONITOR
00653 
00655     mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
00656   
00658     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00659                   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C, 
00660                   bool completeBasis, int howMany = -1 ) const;
00661     
00663     bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX, 
00664                      Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00665                      Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00666 
00668     bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00669                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00670                     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00671 
00685     int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00686                        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00687                        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00688                        Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;    
00689   };
00690 
00692   // Set the label for this orthogonalization manager and create new timers if it's changed
00693   template<class ScalarType, class MV, class OP>
00694   void ICGSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
00695   { 
00696     if (label != label_) {
00697       label_ = label;
00698       std::string orthoLabel = label_ + ": Orthogonalization";
00699 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00700       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00701 #endif
00702 
00703       std::string updateLabel = label_ + ": Ortho (Update)";
00704 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00705       timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00706 #endif
00707 
00708       std::string normLabel = label_ + ": Ortho (Norm)";
00709 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00710       timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00711 #endif
00712 
00713       std::string ipLabel = label_ + ": Ortho (Inner Product)";
00714 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00715       timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel);
00716 #endif
00717     }
00718   } 
00719 
00721   // Compute the distance from orthonormality
00722   template<class ScalarType, class MV, class OP>
00723   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00724   ICGSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
00725     const ScalarType ONE = SCT::one();
00726     int rank = MVT::GetNumberVecs(X);
00727     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00728     MatOrthoManager<ScalarType,MV,OP>::innerProd(X,X,MX,xTx);
00729     for (int i=0; i<rank; i++) {
00730       xTx(i,i) -= ONE;
00731     }
00732     return xTx.normFrobenius();
00733   }
00734 
00736   // Compute the distance from orthogonality
00737   template<class ScalarType, class MV, class OP>
00738   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00739   ICGSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
00740     int r1 = MVT::GetNumberVecs(X1);
00741     int r2  = MVT::GetNumberVecs(X2);
00742     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
00743     MatOrthoManager<ScalarType,MV,OP>::innerProd(X2,X1,MX1,xTx);
00744     return xTx.normFrobenius();
00745   }
00746 
00748   // Find an Op-orthonormal basis for span(X) - span(W)
00749   template<class ScalarType, class MV, class OP>
00750   int ICGSOrthoManager<ScalarType, MV, OP>::projectAndNormalize(
00751                                     MV &X, Teuchos::RCP<MV> MX, 
00752                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00753                                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00754                                     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const 
00755   {
00756     using Teuchos::Array;
00757     using Teuchos::null;
00758     using Teuchos::is_null;
00759     using Teuchos::RCP;
00760     using Teuchos::rcp;
00761     using Teuchos::SerialDenseMatrix;
00762     typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
00763     typedef typename Array< RCP< const MV > >::size_type size_type;
00764 
00765 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00766     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00767 #endif
00768 
00769     ScalarType    ONE  = SCT::one();
00770     ScalarType    ZERO  = SCT::zero();
00771 
00772     int nq = Q.size();
00773     int xc = MVT::GetNumberVecs( X );
00774     int xr = MVT::GetVecLength( X );
00775     int rank = xc;
00776 
00777     // If the user doesn't want to store the normalization
00778     // coefficients, allocate some local memory for them.  This will
00779     // go away at the end of this method.
00780     if (is_null (B)) {
00781       B = rcp (new serial_dense_matrix_type (xc, xc));
00782     }
00783     // Likewise, if the user doesn't want to store the projection
00784     // coefficients, allocate some local memory for them.  Also make
00785     // sure that all the entries of C are the right size.  We're going
00786     // to overwrite them anyway, so we don't have to worry about the
00787     // contents (other than to resize them if they are the wrong
00788     // size).
00789     if (C.size() < nq)
00790       C.resize (nq);
00791     for (size_type k = 0; k < nq; ++k)
00792       {
00793         const int numRows = MVT::GetNumberVecs (*Q[k]);
00794         const int numCols = xc; // Number of vectors in X
00795 
00796         if (is_null (C[k]))
00797           C[k] = rcp (new serial_dense_matrix_type (numRows, numCols));
00798         else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
00799         {
00800           int err = C[k]->reshape (numRows, numCols);
00801           TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, 
00802               "IMGS orthogonalization: failed to reshape "
00803               "C[" << k << "] (the array of block "
00804               "coefficients resulting from projecting X "
00805               "against Q[1:" << nq << "]).");
00806         }
00807       }
00808 
00809     /******   DO NOT MODIFY *MX IF _hasOp == false   ******/
00810     if (this->_hasOp) {
00811       if (MX == Teuchos::null) {
00812         // we need to allocate space for MX
00813         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00814         OPT::Apply(*(this->_Op),X,*MX);
00815       }
00816     }
00817     else {
00818       // Op == I  -->  MX = X (ignore it if the user passed it in)
00819       MX = Teuchos::rcp( &X, false );
00820     }
00821 
00822     int mxc = MVT::GetNumberVecs( *MX );
00823     int mxr = MVT::GetVecLength( *MX );
00824 
00825     // short-circuit
00826     TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::ICGSOrthoManager::projectAndNormalize(): X must be non-empty" );
00827 
00828     int numbas = 0;
00829     for (int i=0; i<nq; i++) {
00830       numbas += MVT::GetNumberVecs( *Q[i] );
00831     }
00832 
00833     // check size of B
00834     TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00835                         "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
00836     // check size of X and MX
00837     TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00838                         "Belos::ICGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
00839     // check size of X w.r.t. MX 
00840     TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00841                         "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
00842     // check feasibility
00843     //TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00844     //                    "Belos::ICGSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
00845 
00846     // Some flags for checking dependency returns from the internal orthogonalization methods
00847     bool dep_flg = false;
00848 
00849     // Make a temporary copy of X and MX, just in case a block dependency is detected.
00850     Teuchos::RCP<MV> tmpX, tmpMX;
00851     tmpX = MVT::CloneCopy(X);
00852     if (this->_hasOp) {
00853       tmpMX = MVT::CloneCopy(*MX);
00854     }
00855 
00856     if (xc == 1) {
00857 
00858       // Use the cheaper block orthogonalization.
00859       // NOTE: Don't check for dependencies because the update has one vector.
00860       dep_flg = blkOrtho1( X, MX, C, Q );
00861 
00862       // Normalize the new block X
00863       {
00864 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00865         Teuchos::TimeMonitor normTimer( *timerNorm_ );
00866 #endif
00867         if ( B == Teuchos::null ) {
00868           B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00869         }
00870         std::vector<ScalarType> diag(xc);
00871         MVT::MvDot( X, *MX, diag );
00872         (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
00873         rank = 1; 
00874         MVT::MvAddMv( ONE/(*B)(0,0), X, ZERO, X, X );
00875         if (this->_hasOp) {
00876           // Update MXj.
00877           MVT::MvAddMv( ONE/(*B)(0,0), *MX, ZERO, *MX, *MX );
00878         }
00879       }
00880 
00881     } 
00882     else {
00883 
00884       // Use the cheaper block orthogonalization.
00885       dep_flg = blkOrtho( X, MX, C, Q );
00886 
00887       // If a dependency has been detected in this block, then perform
00888       // the more expensive nonblock (single vector at a time)
00889       // orthogonalization.
00890       if (dep_flg) {
00891         rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00892 
00893         // Copy tmpX back into X.
00894         MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00895         if (this->_hasOp) {
00896           MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00897         }
00898       } 
00899       else {
00900         // There is no dependency, so orthonormalize new block X
00901         rank = findBasis( X, MX, B, false );
00902         if (rank < xc) {
00903           // A dependency was found during orthonormalization of X,
00904           // rerun orthogonalization using more expensive nonblock
00905           // orthogonalization.
00906           rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00907 
00908           // Copy tmpX back into X.
00909           MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00910           if (this->_hasOp) {
00911             MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00912           }
00913         }    
00914       }
00915     } // if (xc == 1) {
00916 
00917     // this should not raise an std::exception; but our post-conditions oblige us to check
00918     TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00919                         "Belos::ICGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
00920 
00921     // Return the rank of X.
00922     return rank;
00923   }
00924 
00925 
00926 
00928   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00929   template<class ScalarType, class MV, class OP>
00930   int ICGSOrthoManager<ScalarType, MV, OP>::normalize(
00931                                 MV &X, Teuchos::RCP<MV> MX, 
00932                                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00933 
00934 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00935     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00936 #endif
00937 
00938     // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
00939     return findBasis(X, MX, B, true);
00940   }
00941 
00942 
00943 
00945   template<class ScalarType, class MV, class OP>
00946   void ICGSOrthoManager<ScalarType, MV, OP>::project(
00947                           MV &X, Teuchos::RCP<MV> MX, 
00948                           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00949                           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00950     // For the inner product defined by the operator Op or the identity (Op == 0)
00951     //   -> Orthogonalize X against each Q[i]
00952     // Modify MX accordingly
00953     //
00954     // Note that when Op is 0, MX is not referenced
00955     //
00956     // Parameter variables
00957     //
00958     // X  : Vectors to be transformed
00959     //
00960     // MX : Image of the block of vectors X by the mass matrix
00961     //
00962     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
00963     //
00964     
00965 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00966     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00967 #endif
00968     
00969     int xc = MVT::GetNumberVecs( X );
00970     int xr = MVT::GetVecLength( X );
00971     int nq = Q.size();
00972     std::vector<int> qcs(nq);
00973     // short-circuit
00974     if (nq == 0 || xc == 0 || xr == 0) {
00975       return;
00976     }
00977     int qr = MVT::GetVecLength ( *Q[0] );
00978     // if we don't have enough C, expand it with null references
00979     // if we have too many, resize to throw away the latter ones
00980     // if we have exactly as many as we have Q, this call has no effect
00981     C.resize(nq);
00982 
00983 
00984     /******   DO NOT MODIFY *MX IF _hasOp == false   ******/
00985     if (this->_hasOp) {
00986       if (MX == Teuchos::null) {
00987         // we need to allocate space for MX
00988         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00989         OPT::Apply(*(this->_Op),X,*MX);
00990       }
00991     }
00992     else {
00993       // Op == I  -->  MX = X (ignore it if the user passed it in)
00994       MX = Teuchos::rcp( &X, false );
00995     }
00996     int mxc = MVT::GetNumberVecs( *MX );
00997     int mxr = MVT::GetVecLength( *MX );
00998 
00999     // check size of X and Q w.r.t. common sense
01000     TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
01001                         "Belos::ICGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
01002     // check size of X w.r.t. MX and Q
01003     TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
01004                         "Belos::ICGSOrthoManager::project(): Size of X not consistant with MX,Q" );
01005 
01006     // tally up size of all Q and check/allocate C
01007     int baslen = 0;
01008     for (int i=0; i<nq; i++) {
01009       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 
01010                           "Belos::ICGSOrthoManager::project(): Q lengths not mutually consistant" );
01011       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01012       TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
01013                           "Belos::ICGSOrthoManager::project(): Q has less rows than columns" );
01014       baslen += qcs[i];
01015 
01016       // check size of C[i]
01017       if ( C[i] == Teuchos::null ) {
01018         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
01019       }
01020       else {
01021         TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
01022                            "Belos::ICGSOrthoManager::project(): Size of Q not consistant with size of C" );
01023       }
01024     }
01025 
01026     // Use the cheaper block orthogonalization, don't check for rank deficiency.
01027     blkOrtho( X, MX, C, Q );
01028 
01029   }  
01030 
01032   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
01033   // the rank is numvectors(X)
01034   template<class ScalarType, class MV, class OP>
01035   int ICGSOrthoManager<ScalarType, MV, OP>::findBasis(
01036                 MV &X, Teuchos::RCP<MV> MX, 
01037                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
01038                 bool completeBasis, int howMany ) const {
01039     // For the inner product defined by the operator Op or the identity (Op == 0)
01040     //   -> Orthonormalize X 
01041     // Modify MX accordingly
01042     //
01043     // Note that when Op is 0, MX is not referenced
01044     //
01045     // Parameter variables
01046     //
01047     // X  : Vectors to be orthonormalized
01048     //
01049     // MX : Image of the multivector X under the operator Op
01050     //
01051     // Op  : Pointer to the operator for the inner product
01052     //
01053     //
01054 
01055 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01056     Teuchos::TimeMonitor normTimer( *timerNorm_ );
01057 #endif
01058 
01059     const ScalarType ONE  = SCT::one();
01060     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
01061 
01062     int xc = MVT::GetNumberVecs( X );
01063     int xr = MVT::GetVecLength( X );
01064 
01065     if (howMany == -1) {
01066       howMany = xc;
01067     }
01068 
01069     /*******************************************************
01070      *  If _hasOp == false, we will not reference MX below *
01071      *******************************************************/
01072 
01073     // if Op==null, MX == X (via pointer)
01074     // Otherwise, either the user passed in MX or we will allocated and compute it
01075     if (this->_hasOp) {
01076       if (MX == Teuchos::null) {
01077         // we need to allocate space for MX
01078         MX = MVT::Clone(X,xc);
01079         OPT::Apply(*(this->_Op),X,*MX);
01080       }
01081     }
01082 
01083     /* if the user doesn't want to store the coefficienets, 
01084      * allocate some local memory for them 
01085      */
01086     if ( B == Teuchos::null ) {
01087       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
01088     }
01089 
01090     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
01091     int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX )  : xr;
01092 
01093     // check size of C, B
01094     TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
01095                         "Belos::ICGSOrthoManager::findBasis(): X must be non-empty" );
01096     TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
01097                         "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
01098     TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
01099                         "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
01100     TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 
01101                         "Belos::ICGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
01102     TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
01103                         "Belos::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
01104 
01105     /* xstart is which column we are starting the process with, based on howMany
01106      * columns before xstart are assumed to be Op-orthonormal already
01107      */
01108     int xstart = xc - howMany;
01109 
01110     for (int j = xstart; j < xc; j++) {
01111 
01112       // numX is 
01113       // * number of currently orthonormal columns of X
01114       // * the index of the current column of X
01115       int numX = j;
01116       bool addVec = false;
01117 
01118       // Get a view of the vector currently being worked on.
01119       std::vector<int> index(1);
01120       index[0] = numX;
01121       Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
01122       Teuchos::RCP<MV> MXj;
01123       if ((this->_hasOp)) {
01124         // MXj is a view of the current vector in MX
01125         MXj = MVT::CloneViewNonConst( *MX, index );
01126       }
01127       else {
01128         // MXj is a pointer to Xj, and MUST NOT be modified
01129         MXj = Xj;
01130       }
01131 
01132       // Get a view of the previous vectors.
01133       std::vector<int> prev_idx( numX );
01134       Teuchos::RCP<const MV> prevX, prevMX;
01135 
01136       if (numX > 0) {
01137         for (int i=0; i<numX; i++) {
01138           prev_idx[i] = i;
01139         }
01140         prevX = MVT::CloneView( X, prev_idx );
01141         if (this->_hasOp) {
01142           prevMX = MVT::CloneView( *MX, prev_idx );
01143         }
01144       } 
01145 
01146       // Make storage for these Gram-Schmidt iterations.
01147       Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
01148       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01149       //
01150       // Save old MXj vector and compute Op-norm
01151       //
01152       Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
01153       MVT::MvDot( *Xj, *MXj, oldDot );
01154       // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
01155       TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
01156           "Belos::ICGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
01157 
01158       if (numX > 0) {
01159 
01160         Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
01161 
01162         for (int i=0; i<max_ortho_steps_; ++i) {
01163 
01164           // product <- prevX^T MXj
01165           {
01166 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01167             Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01168 #endif
01169             MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,P2);
01170           }
01171 
01172           // Xj <- Xj - prevX prevX^T MXj   
01173           //     = Xj - prevX product
01174           {
01175 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01176             Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01177 #endif
01178             MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
01179           }
01180 
01181           // Update MXj
01182           if (this->_hasOp) {
01183             // MXj <- Op*Xj_new
01184             //      = Op*(Xj_old - prevX prevX^T MXj)
01185             //      = MXj - prevMX product
01186 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01187             Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01188 #endif
01189             MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
01190           }
01191 
01192           // Set coefficients
01193           if ( i==0 )
01194             product = P2;
01195           else
01196             product += P2;
01197         }
01198 
01199       } // if (numX > 0)
01200 
01201       // Compute Op-norm with old MXj
01202       MVT::MvDot( *Xj, *oldMXj, newDot );
01203 
01204       // Check to see if the new vector is dependent.
01205       if (completeBasis) {
01206         //
01207         // We need a complete basis, so add random vectors if necessary
01208         //
01209         if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
01210 
01211           // Add a random vector and orthogonalize it against previous vectors in block.
01212           addVec = true;
01213 #ifdef ORTHO_DEBUG
01214           std::cout << "Belos::ICGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
01215 #endif
01216           //
01217           Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01218           Teuchos::RCP<MV> tempMXj;
01219           MVT::MvRandom( *tempXj );
01220           if (this->_hasOp) {
01221             tempMXj = MVT::Clone( X, 1 );
01222             OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01223           } 
01224           else {
01225             tempMXj = tempXj;
01226           }
01227           MVT::MvDot( *tempXj, *tempMXj, oldDot );
01228           //
01229           for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
01230             {
01231 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01232               Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01233 #endif
01234               MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,product);
01235             }
01236             {
01237 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01238               Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01239 #endif
01240               MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
01241             }
01242             if (this->_hasOp) {
01243 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01244               Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01245 #endif
01246               MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
01247             }
01248           }
01249           // Compute new Op-norm
01250           MVT::MvDot( *tempXj, *tempMXj, newDot );
01251           //
01252           if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01253             // Copy vector into current column of _basisvecs
01254             MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
01255             if (this->_hasOp) {
01256               MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
01257             }
01258           }
01259           else {
01260             return numX;
01261           } 
01262         }
01263       }
01264       else {
01265         //
01266         // We only need to detect dependencies.
01267         //
01268         if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
01269           return numX;
01270         }
01271       }
01272 
01273       // If we haven't left this method yet, then we can normalize the new vector Xj.
01274       // Normalize Xj.
01275       // Xj <- Xj / std::sqrt(newDot)
01276       ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01277       {
01278         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01279         if (this->_hasOp) {
01280           // Update MXj.
01281           MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01282         }
01283       }
01284 
01285       // If we've added a random vector, enter a zero in the j'th diagonal element.
01286       if (addVec) {
01287         (*B)(j,j) = ZERO;
01288       }
01289       else {
01290         (*B)(j,j) = diag;
01291       }
01292 
01293       // Save the coefficients, if we are working on the original vector and not a randomly generated one
01294       if (!addVec) {
01295         for (int i=0; i<numX; i++) {
01296           (*B)(i,j) = product(i,0);
01297         }
01298       }
01299 
01300     } // for (j = 0; j < xc; ++j)
01301 
01302     return xc;
01303   }
01304 
01306   // Routine to compute the block orthogonalization
01307   template<class ScalarType, class MV, class OP>
01308   bool 
01309   ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX, 
01310                                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01311                                                     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
01312   {
01313     int nq = Q.size();
01314     int xc = MVT::GetNumberVecs( X );
01315     const ScalarType ONE  = SCT::one();
01316 
01317     std::vector<int> qcs( nq );
01318     for (int i=0; i<nq; i++) {
01319       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01320     }
01321 
01322     // Perform the Gram-Schmidt transformation for a block of vectors
01323 
01324     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
01325     // Define the product Q^T * (Op*X)
01326     for (int i=0; i<nq; i++) {
01327       // Multiply Q' with MX
01328       {
01329 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01330         Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01331 #endif
01332         MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,*C[i]);
01333       }
01334       // Multiply by Q and subtract the result in X
01335       {
01336 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01337         Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01338 #endif
01339         MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
01340       }
01341 
01342       // Update MX, with the least number of applications of Op as possible
01343       if (this->_hasOp) {
01344         if (xc <= qcs[i]) {
01345           OPT::Apply( *(this->_Op), X, *MX);
01346         }
01347         else {
01348           // this will possibly be used again below; don't delete it
01349           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01350           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01351           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01352         }
01353       }
01354     }
01355 
01356     // Do as many steps of classical Gram-Schmidt as required by max_ortho_steps_
01357     for (int j = 1; j < max_ortho_steps_; ++j) {
01358       
01359       for (int i=0; i<nq; i++) {
01360         Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
01361 
01362         // Apply another step of classical Gram-Schmidt
01363         {
01364 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01365           Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01366 #endif
01367           MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,C2);
01368         }
01369         *C[i] += C2;
01370         {
01371 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01372           Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01373 #endif
01374           MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
01375         }
01376 
01377         // Update MX, with the least number of applications of Op as possible
01378         if (this->_hasOp) {
01379           if (MQ[i].get()) {
01380             // MQ was allocated and computed above; use it
01381             MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01382           }
01383           else if (xc <= qcs[i]) {
01384             // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01385             OPT::Apply( *(this->_Op), X, *MX);
01386           }
01387         }
01388       } // for (int i=0; i<nq; i++)
01389     } // for (int j = 0; j < max_ortho_steps; ++j)
01390 
01391     return false;
01392   }
01393 
01395   // Routine to compute the block orthogonalization
01396   template<class ScalarType, class MV, class OP>
01397   bool 
01398   ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
01399                                                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01400                                                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
01401   {
01402     int nq = Q.size();
01403     int xc = MVT::GetNumberVecs( X );
01404     bool dep_flg = false;
01405     const ScalarType ONE  = SCT::one();
01406 
01407     std::vector<int> qcs( nq );
01408     for (int i=0; i<nq; i++) {
01409       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01410     }
01411 
01412     // Perform the Gram-Schmidt transformation for a block of vectors
01413     
01414     // Compute the initial Op-norms
01415     std::vector<ScalarType> oldDot( xc );
01416     MVT::MvDot( X, *MX, oldDot );
01417 
01418     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
01419     // Define the product Q^T * (Op*X)
01420     for (int i=0; i<nq; i++) {
01421       // Multiply Q' with MX
01422       {
01423 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01424         Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01425 #endif
01426         MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,*C[i]);
01427       }
01428       // Multiply by Q and subtract the result in X
01429       {
01430 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01431         Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01432 #endif
01433         MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
01434       }
01435       // Update MX, with the least number of applications of Op as possible
01436       if (this->_hasOp) {
01437         if (xc <= qcs[i]) {
01438           OPT::Apply( *(this->_Op), X, *MX);
01439         }
01440         else {
01441           // this will possibly be used again below; don't delete it
01442           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01443           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01444           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01445         }
01446       }
01447     }
01448 
01449     // Do as many steps of classical Gram-Schmidt as required by max_ortho_steps_
01450     for (int j = 1; j < max_ortho_steps_; ++j) {
01451       
01452       for (int i=0; i<nq; i++) {
01453         Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
01454 
01455         // Apply another step of classical Gram-Schmidt
01456         {
01457 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01458           Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01459 #endif
01460           MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,C2);
01461         }
01462         *C[i] += C2;
01463         {
01464 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01465           Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01466 #endif
01467           MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
01468         }
01469 
01470         // Update MX, with the least number of applications of Op as possible
01471         if (this->_hasOp) {
01472           if (MQ[i].get()) {
01473             // MQ was allocated and computed above; use it
01474             MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01475           }
01476           else if (xc <= qcs[i]) {
01477             // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01478             OPT::Apply( *(this->_Op), X, *MX);
01479           }
01480         }
01481       } // for (int i=0; i<nq; i++)
01482     } // for (int j = 0; j < max_ortho_steps; ++j)
01483 
01484     // Compute new Op-norms
01485     std::vector<ScalarType> newDot(xc);
01486     MVT::MvDot( X, *MX, newDot );
01487 
01488     // Check to make sure the new block of vectors are not dependent on previous vectors
01489     for (int i=0; i<xc; i++){
01490       if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
01491         dep_flg = true;
01492         break;
01493       }
01494     } // end for (i=0;...)
01495 
01496     return dep_flg;
01497   }
01498   
01499   template<class ScalarType, class MV, class OP>
01500   int
01501   ICGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
01502                                                        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01503                                                        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
01504                                                        Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
01505   {
01506     Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
01507 
01508     const ScalarType ONE  = SCT::one();
01509     const ScalarType ZERO  = SCT::zero();
01510     
01511     int nq = Q.size();
01512     int xc = MVT::GetNumberVecs( X );
01513     std::vector<int> indX( 1 );
01514     std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01515 
01516     std::vector<int> qcs( nq );
01517     for (int i=0; i<nq; i++) {
01518       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01519     }
01520 
01521     // Create pointers for the previous vectors of X that have already been orthonormalized.
01522     Teuchos::RCP<const MV> lastQ;
01523     Teuchos::RCP<MV> Xj, MXj;
01524     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
01525 
01526     // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
01527     for (int j=0; j<xc; j++) {
01528       
01529       bool dep_flg = false;
01530       
01531       // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
01532       if (j > 0) {
01533         std::vector<int> index( j );
01534         for (int ind=0; ind<j; ind++) {
01535           index[ind] = ind;
01536         }
01537         lastQ = MVT::CloneView( X, index );
01538 
01539         // Add these views to the Q and C arrays.
01540         Q.push_back( lastQ );
01541         C.push_back( B );
01542         qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
01543       }
01544 
01545       // Get a view of the current vector in X to orthogonalize.
01546       indX[0] = j;
01547       Xj = MVT::CloneViewNonConst( X, indX );
01548       if (this->_hasOp) {
01549         MXj = MVT::CloneViewNonConst( *MX, indX );
01550       }
01551       else {
01552         MXj = Xj;
01553       }
01554 
01555       // Compute the initial Op-norms
01556       MVT::MvDot( *Xj, *MXj, oldDot );
01557 
01558       Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
01559       // Define the product Q^T * (Op*X)
01560       for (int i=0; i<Q.size(); i++) {
01561 
01562         // Get a view of the current serial dense matrix
01563         Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01564 
01565         // Multiply Q' with MX
01566         MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,tempC);
01567         // Multiply by Q and subtract the result in Xj
01568         MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
01569 
01570         // Update MXj, with the least number of applications of Op as possible
01571         if (this->_hasOp) {
01572           if (xc <= qcs[i]) {
01573             OPT::Apply( *(this->_Op), *Xj, *MXj);
01574           }
01575           else {
01576             // this will possibly be used again below; don't delete it
01577             MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01578             OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01579             MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
01580           }
01581         }
01582       }
01583 
01584       // Do any additional steps of classical Gram-Schmidt orthogonalization 
01585       for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
01586 
01587         for (int i=0; i<Q.size(); i++) {
01588           Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01589           Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
01590 
01591           // Apply another step of classical Gram-Schmidt
01592           MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,C2);
01593           tempC += C2;
01594           MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
01595 
01596           // Update MXj, with the least number of applications of Op as possible
01597           if (this->_hasOp) {
01598             if (MQ[i].get()) {
01599               // MQ was allocated and computed above; use it
01600               MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
01601             }
01602             else if (xc <= qcs[i]) {
01603               // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01604               OPT::Apply( *(this->_Op), *Xj, *MXj);
01605             }
01606           }
01607         } // for (int i=0; i<Q.size(); i++)
01608 
01609       } // for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps)
01610 
01611       // Check for linear dependence.
01612       if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
01613         dep_flg = true;
01614       }
01615 
01616       // Normalize the new vector if it's not dependent
01617       if (!dep_flg) {
01618         ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01619 
01620         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01621         if (this->_hasOp) {
01622           // Update MXj.
01623           MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01624         }
01625 
01626         // Enter value on diagonal of B.
01627         (*B)(j,j) = diag;
01628       }
01629       else {
01630         // Create a random vector and orthogonalize it against all previous columns of Q.
01631         Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01632         Teuchos::RCP<MV> tempMXj;
01633         MVT::MvRandom( *tempXj );
01634         if (this->_hasOp) {
01635           tempMXj = MVT::Clone( X, 1 );
01636           OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01637         } 
01638         else {
01639           tempMXj = tempXj;
01640         }
01641         MVT::MvDot( *tempXj, *tempMXj, oldDot );
01642         //
01643         for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
01644 
01645           for (int i=0; i<Q.size(); i++) {
01646             Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
01647 
01648             // Apply another step of classical Gram-Schmidt
01649             MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*tempXj,tempMXj,product);
01650             MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
01651 
01652             // Update MXj, with the least number of applications of Op as possible
01653             if (this->_hasOp) {
01654               if (MQ[i].get()) {
01655                 // MQ was allocated and computed above; use it
01656                 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
01657               }
01658               else if (xc <= qcs[i]) {
01659                 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01660                 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
01661               }
01662             }
01663           } // for (int i=0; i<nq; i++)
01664 
01665         }
01666 
01667         // Compute the Op-norms after the correction step.
01668         MVT::MvDot( *tempXj, *tempMXj, newDot );
01669 
01670         // Copy vector into current column of Xj
01671         if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01672           ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01673 
01674           // Enter value on diagonal of B.
01675           (*B)(j,j) = ZERO;
01676 
01677           // Copy vector into current column of _basisvecs
01678           MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
01679           if (this->_hasOp) {
01680             MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
01681           }
01682         }
01683         else {
01684           return j;
01685         } 
01686       } // if (!dep_flg)
01687 
01688       // Remove the vectors from array
01689       if (j > 0) {
01690         Q.resize( nq );
01691         C.resize( nq );
01692         qcs.resize( nq );
01693       }
01694 
01695     } // for (int j=0; j<xc; j++)
01696 
01697     return xc;
01698   }
01699 
01700 } // namespace Belos
01701 
01702 #endif // BELOS_ICGS_ORTHOMANAGER_HPP
01703 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines