Belos Version of the Day
BelosDGKSOrthoManager.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_DGKS_ORTHOMANAGER_HPP
00048 #define BELOS_DGKS_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   getDefaultDgksParameters()
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 DGKS orthogonalization.
00085     // Documentation will be embedded in the parameter list.
00086     const int defaultMaxNumOrthogPasses = 2;
00087     const magnitude_type one = STM::one();
00088     const magnitude_type eps = STM::eps();
00089     const magnitude_type squareRootTwo = 
00090       STM::squareroot (as<magnitude_type> (2));
00091     const magnitude_type defaultBlkTol = 
00092       as<magnitude_type> (10) * STM::squareroot (eps);
00093     const magnitude_type defaultDepTol = one / squareRootTwo;
00094     const magnitude_type defaultSingTol = as<magnitude_type> (10) * eps;
00095 
00096     Teuchos::RCP<Teuchos::ParameterList> params = 
00097       Teuchos::parameterList ("DGKS");
00098     params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
00099      "Maximum number of orthogonalization passes "
00100      "(includes the first).  Default is 2, since "
00101      "\"twice is enough\" for Krylov methods.");
00102     params->set ("blkTol", defaultBlkTol, 
00103      "Block reorthogonalization threshhold.");
00104     params->set ("depTol", defaultDepTol, 
00105      "(Non-block) reorthogonalization threshold.");
00106     params->set ("singTol", defaultSingTol, 
00107      "Singular block detection threshold.");
00108     return params;
00109   }
00110 
00116   template<class ScalarType>
00117   Teuchos::RCP<const Teuchos::ParameterList> TEUCHOS_DEPRECATED
00118   getFastDgksParameters()
00119   {
00120     using Teuchos::ParameterList;
00121     using Teuchos::RCP;
00122     using Teuchos::rcp;
00123     using Teuchos::ScalarTraits;
00124     typedef typename ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00125     typedef ScalarTraits<magnitude_type> STM;
00126 
00127     RCP<const ParameterList> defaultParams = getDefaultDgksParameters<ScalarType> ();
00128     // Start with a clone of the default parameters.
00129     RCP<ParameterList> params = rcp (new ParameterList (*defaultParams));
00130 
00131     const int maxBlkOrtho = 1;
00132     const magnitude_type blkTol = STM::zero();
00133     const magnitude_type depTol = STM::zero();
00134     const magnitude_type singTol = STM::zero();
00135 
00136     params->set ("maxNumOrthogPasses", maxBlkOrtho);
00137     params->set ("blkTol", blkTol);
00138     params->set ("depTol", depTol);
00139     params->set ("singTol", singTol);
00140 
00141     return params;
00142   }
00143 
00153   template<class ScalarType>
00154   void TEUCHOS_DEPRECATED
00155   readDgksParameters (const Teuchos::RCP<const Teuchos::ParameterList>& params,
00156           int& maxNumOrthogPasses,
00157           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& blkTol,
00158           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& depTol,
00159           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& singTol)
00160   {
00161     using Teuchos::ParameterList;
00162     using Teuchos::RCP;
00163     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00164     const magnitude_type zero = Teuchos::ScalarTraits<magnitude_type>::zero();
00165     RCP<const ParameterList> defaultParams = getDefaultDgksParameters<ScalarType>();
00166 
00167     // Using temporary variables and fetching all values before
00168     // setting the output arguments ensures the strong exception
00169     // guarantee for this function: if an exception is thrown, no
00170     // externally visible side effects (in this case, setting the
00171     // output arguments) have taken place.
00172     int _maxNumOrthogPasses;
00173     magnitude_type _blkTol, _depTol, _singTol;
00174     if (params.is_null())
00175       {
00176   _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00177   _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00178   _depTol = defaultParams->get<magnitude_type> ("depTol");  
00179   _singTol = defaultParams->get<magnitude_type> ("singTol");
00180       }
00181     else
00182       {
00183   try {
00184     _maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
00185     if (_maxNumOrthogPasses < 1)
00186       _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00187   } catch (Teuchos::Exceptions::InvalidParameter&) {
00188     _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00189   }
00190 
00191   try {
00192     _blkTol = params->get<magnitude_type> ("blkTol");
00193     if (_blkTol < zero)
00194       _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00195   } catch (Teuchos::Exceptions::InvalidParameter&) {
00196     _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00197   }
00198 
00199   try {
00200     _depTol = params->get<magnitude_type> ("depTol");
00201     if (_depTol < zero)
00202       _depTol = defaultParams->get<magnitude_type> ("depTol");
00203   } catch (Teuchos::Exceptions::InvalidParameter&) {
00204     _depTol = defaultParams->get<magnitude_type> ("depTol");
00205   }
00206 
00207   try {
00208     _singTol = params->get<magnitude_type> ("singTol");
00209     if (_singTol < zero)
00210       _singTol = defaultParams->get<magnitude_type> ("singTol");
00211   } catch (Teuchos::Exceptions::InvalidParameter&) {
00212     _singTol = defaultParams->get<magnitude_type> ("singTol");
00213   }
00214       }
00215     maxNumOrthogPasses = _maxNumOrthogPasses;
00216     blkTol = _blkTol;
00217     depTol = _depTol;
00218     singTol = _singTol;
00219   }
00220 
00221   template<class ScalarType, class MV, class OP>
00222   class DGKSOrthoManager : 
00223     public MatOrthoManager<ScalarType,MV,OP>,
00224     public Teuchos::ParameterListAcceptorDefaultBase
00225   {
00226   private:
00227     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00228     typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
00229     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
00230     typedef MultiVecTraits<ScalarType,MV>      MVT;
00231     typedef MultiVecTraitsExt<ScalarType,MV>   MVText;
00232     typedef OperatorTraits<ScalarType,MV,OP>   OPT;
00233 
00234   public:
00236 
00237 
00239     DGKSOrthoManager( const std::string& label = "Belos",
00240                       Teuchos::RCP<const OP> Op = Teuchos::null,
00241           const int max_blk_ortho = 2,
00242           const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
00243           const MagnitudeType dep_tol = MGT::one()/MGT::squareroot( 2*MGT::one() ),
00244           const MagnitudeType sing_tol = 10*MGT::eps() )
00245       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00246   max_blk_ortho_( max_blk_ortho ),
00247   blk_tol_( blk_tol ),
00248   dep_tol_( dep_tol ),
00249   sing_tol_( sing_tol ),
00250   label_( label )
00251     {
00252 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00253       std::string orthoLabel = label_ + ": Orthogonalization";
00254       timerOrtho_ = Teuchos::TimeMonitor::getNewCounter( orthoLabel );
00255 #endif
00256     }    
00257 
00259     DGKSOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& plist,
00260           const std::string& label = "Belos",
00261                       Teuchos::RCP<const OP> Op = Teuchos::null)
00262       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00263   max_blk_ortho_ (2),
00264   blk_tol_ (Teuchos::as<MagnitudeType>(10) * MGT::squareroot(MGT::eps())),
00265   dep_tol_ (MGT::one() / MGT::squareroot (Teuchos::as<MagnitudeType>(2))),
00266   sing_tol_ (Teuchos::as<MagnitudeType>(10) * MGT::eps()),
00267   label_( label )
00268     {
00269       setParameterList (plist);
00270 
00271 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00272       std::string orthoLabel = label_ + ": Orthogonalization";
00273       timerOrtho_ = Teuchos::TimeMonitor::getNewCounter( orthoLabel );
00274 #endif
00275     }    
00276 
00278     ~DGKSOrthoManager() {}
00280 
00282 
00283 
00284     void 
00285     setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
00286     {
00287       using Teuchos::ParameterList;
00288       using Teuchos::parameterList;
00289       using Teuchos::RCP;
00290 
00291       RCP<const ParameterList> defaultParams = getValidParameters();
00292       RCP<ParameterList> params;
00293       if (plist.is_null()) {
00294   // No need to validate default parameters.
00295   params = parameterList (*defaultParams);
00296       } else {
00297   params = plist;
00298   params->validateParametersAndSetDefaults (*defaultParams);
00299       }
00300   
00301       // Using temporary variables and fetching all values before
00302       // setting the output arguments ensures the strong exception
00303       // guarantee for this function: if an exception is thrown, no
00304       // externally visible side effects (in this case, setting the
00305       // output arguments) have taken place.
00306       const int maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
00307       const MagnitudeType blkTol = params->get<MagnitudeType> ("blkTol");
00308       const MagnitudeType depTol = params->get<MagnitudeType> ("depTol");
00309       const MagnitudeType singTol = params->get<MagnitudeType> ("singTol");
00310 
00311       max_blk_ortho_ = maxNumOrthogPasses;
00312       blk_tol_ = blkTol;
00313       dep_tol_ = depTol;
00314       sing_tol_ = singTol;
00315 
00316       setMyParamList (params);
00317     }
00318 
00319     Teuchos::RCP<const Teuchos::ParameterList>
00320     getValidParameters () const
00321     {
00322       using Teuchos::as;
00323       using Teuchos::ParameterList;
00324       using Teuchos::parameterList;
00325       using Teuchos::RCP;
00326 
00327       if (defaultParams_.is_null()) {
00328   RCP<ParameterList> params = parameterList ("DGKS");
00329   const MagnitudeType eps = MGT::eps ();
00330 
00331   // Default parameter values for DGKS orthogonalization.
00332   // Documentation will be embedded in the parameter list.
00333   const int defaultMaxNumOrthogPasses = 2;
00334   const MagnitudeType defaultBlkTol = 
00335     as<MagnitudeType> (10) * MGT::squareroot (eps);
00336   const MagnitudeType defaultDepTol =
00337     MGT::one() / MGT::squareroot (as<MagnitudeType> (2));
00338 
00339   const MagnitudeType defaultSingTol = as<MagnitudeType> (10) * eps;
00340 
00341   params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
00342          "Maximum number of orthogonalization passes (includes the "
00343          "first).  Default is 2, since \"twice is enough\" for Krylov "
00344          "methods.");
00345   params->set ("blkTol", defaultBlkTol, "Block reorthogonalization "
00346          "threshhold.");
00347   params->set ("depTol", defaultDepTol, 
00348          "(Non-block) reorthogonalization threshold.");
00349   params->set ("singTol", defaultSingTol, "Singular block detection "
00350          "threshold.");
00351   defaultParams_ = params;
00352       }
00353       return defaultParams_;
00354     }
00355 
00357 
00358     Teuchos::RCP<const Teuchos::ParameterList>
00359     getFastParameters() const
00360     {
00361       using Teuchos::ParameterList;
00362       using Teuchos::RCP;
00363       using Teuchos::rcp;
00364 
00365       RCP<const ParameterList> defaultParams = getValidParameters ();
00366       // Start with a clone of the default parameters.
00367       RCP<ParameterList> params = rcp (new ParameterList (*defaultParams));
00368 
00369       const int maxBlkOrtho = 1;
00370       const MagnitudeType blkTol = MGT::zero();
00371       const MagnitudeType depTol = MGT::zero();
00372       const MagnitudeType singTol = MGT::zero();
00373 
00374       params->set ("maxNumOrthogPasses", maxBlkOrtho);
00375       params->set ("blkTol", blkTol);
00376       params->set ("depTol", depTol);
00377       params->set ("singTol", singTol);
00378 
00379       return params;
00380     }
00381 
00383 
00384 
00386     void setBlkTol( const MagnitudeType blk_tol ) { 
00387       // Update the parameter list as well.
00388       Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
00389       if (! params.is_null()) {
00390   // If it's null, then we haven't called setParameterList()
00391   // yet.  It's entirely possible to construct the parameter
00392   // list on demand, so we don't try to create the parameter
00393   // list here.
00394   params->set ("blkTol", blk_tol);
00395       }
00396       blk_tol_ = blk_tol; 
00397     }
00398 
00400     void setDepTol( const MagnitudeType dep_tol ) { 
00401       // Update the parameter list as well.
00402       Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
00403       if (! params.is_null()) {
00404   params->set ("depTol", dep_tol);
00405       }
00406       dep_tol_ = dep_tol; 
00407     }
00408 
00410     void setSingTol( const MagnitudeType sing_tol ) { 
00411       // Update the parameter list as well.
00412       Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
00413       if (! params.is_null()) {
00414   params->set ("singTol", sing_tol);
00415       }
00416       sing_tol_ = sing_tol; 
00417     }
00418 
00420     MagnitudeType getBlkTol() const { return blk_tol_; } 
00421 
00423     MagnitudeType getDepTol() const { return dep_tol_; } 
00424 
00426     MagnitudeType getSingTol() const { return sing_tol_; } 
00427 
00429 
00430 
00432 
00433 
00461     void project ( MV &X, Teuchos::RCP<MV> MX, 
00462                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00463                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00464 
00465 
00468     void project ( MV &X, 
00469                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00470                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00471       project(X,Teuchos::null,C,Q);
00472     }
00473 
00474 
00475  
00500     int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00501                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
00502 
00503 
00506     int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00507       return normalize(X,Teuchos::null,B);
00508     }
00509 
00510 
00566     int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00567                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00568                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00569                               Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00570 
00573     int projectAndNormalize ( MV &X, 
00574                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00575                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00576                               Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const {
00577       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00578     }
00579 
00581 
00583 
00584 
00590     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00591     orthonormError(const MV &X) const {
00592       return orthonormError(X,Teuchos::null);
00593     }
00594 
00601     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00602     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
00603 
00607     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00608     orthogError(const MV &X1, const MV &X2) const {
00609       return orthogError(X1,Teuchos::null,X2);
00610     }
00611 
00616     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00617     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
00618 
00620 
00622 
00623 
00626     void setLabel(const std::string& label);
00627 
00630     const std::string& getLabel() const { return label_; }
00631 
00633 
00634   private:
00635     
00637     int max_blk_ortho_;
00639     MagnitudeType blk_tol_;
00641     MagnitudeType dep_tol_;
00643     MagnitudeType sing_tol_;
00644 
00646     std::string label_;
00647 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00648     Teuchos::RCP<Teuchos::Time> timerOrtho_;
00649 #endif // BELOS_TEUCHOS_TIME_MONITOR
00650 
00652     mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
00653 
00655     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00656       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C, 
00657       bool completeBasis, int howMany = -1 ) const;
00658     
00660     bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00661         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00662         Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00663 
00677     int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00678            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00679            Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00680            Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;    
00681   };
00682 
00684   // Set the label for this orthogonalization manager and create new timers if it's changed
00685   template<class ScalarType, class MV, class OP>
00686   void DGKSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
00687   {
00688     if (label != label_) {
00689       label_ = label;
00690       std::string orthoLabel = label_ + ": Orthogonalization";
00691 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00692       timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
00693 #endif
00694     }
00695   }
00696 
00698   // Compute the distance from orthonormality
00699   template<class ScalarType, class MV, class OP>
00700   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00701   DGKSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
00702     const ScalarType ONE = SCT::one();
00703     int rank = MVT::GetNumberVecs(X);
00704     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00705     MatOrthoManager<ScalarType,MV,OP>::innerProd(X,X,MX,xTx);
00706     for (int i=0; i<rank; i++) {
00707       xTx(i,i) -= ONE;
00708     }
00709     return xTx.normFrobenius();
00710   }
00711 
00713   // Compute the distance from orthogonality
00714   template<class ScalarType, class MV, class OP>
00715   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00716   DGKSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
00717     int r1 = MVT::GetNumberVecs(X1);
00718     int r2  = MVT::GetNumberVecs(X2);
00719     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
00720     MatOrthoManager<ScalarType,MV,OP>::innerProd(X2,X1,MX1,xTx);
00721     return xTx.normFrobenius();
00722   }
00723 
00725   // Find an Op-orthonormal basis for span(X) - span(W)
00726   template<class ScalarType, class MV, class OP>
00727   int DGKSOrthoManager<ScalarType, MV, OP>::projectAndNormalize(
00728                                     MV &X, Teuchos::RCP<MV> MX, 
00729                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00730                                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00731                                     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const 
00732   {
00733     using Teuchos::Array;
00734     using Teuchos::null;
00735     using Teuchos::is_null;
00736     using Teuchos::RCP;
00737     using Teuchos::rcp;
00738     using Teuchos::SerialDenseMatrix;
00739     typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
00740     typedef typename Array< RCP< const MV > >::size_type size_type;
00741     
00742 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00743     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00744 #endif
00745 
00746     ScalarType    ONE  = SCT::one();
00747     ScalarType    ZERO  = SCT::zero();
00748 
00749     int nq = Q.size();
00750     int xc = MVT::GetNumberVecs( X );
00751     ptrdiff_t xr = MVText::GetGlobalLength( X );
00752     int rank = xc;
00753 
00754     // If the user doesn't want to store the normalization
00755     // coefficients, allocate some local memory for them.  This will
00756     // go away at the end of this method.
00757     if (is_null (B)) {
00758       B = rcp (new serial_dense_matrix_type (xc, xc));
00759     }
00760     // Likewise, if the user doesn't want to store the projection
00761     // coefficients, allocate some local memory for them.  Also make
00762     // sure that all the entries of C are the right size.  We're going
00763     // to overwrite them anyway, so we don't have to worry about the
00764     // contents (other than to resize them if they are the wrong
00765     // size).
00766     if (C.size() < nq)
00767       C.resize (nq);
00768     for (size_type k = 0; k < nq; ++k)
00769       {
00770   const int numRows = MVT::GetNumberVecs (*Q[k]);
00771   const int numCols = xc; // Number of vectors in X
00772   
00773   if (is_null (C[k]))
00774     C[k] = rcp (new serial_dense_matrix_type (numRows, numCols));
00775   else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
00776     {
00777       int err = C[k]->reshape (numRows, numCols);
00778       TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, 
00779              "DGKS orthogonalization: failed to reshape "
00780              "C[" << k << "] (the array of block "
00781              "coefficients resulting from projecting X "
00782              "against Q[1:" << nq << "]).");
00783     }
00784       }
00785 
00786     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00787     if (this->_hasOp) {
00788       if (MX == Teuchos::null) {
00789         // we need to allocate space for MX
00790         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00791         OPT::Apply(*(this->_Op),X,*MX);
00792       }
00793     }
00794     else {
00795       // Op == I  -->  MX = X (ignore it if the user passed it in)
00796       MX = Teuchos::rcp( &X, false );
00797     }
00798 
00799     int mxc = MVT::GetNumberVecs( *MX );
00800     ptrdiff_t mxr = MVText::GetGlobalLength( *MX );
00801 
00802     // short-circuit
00803     TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
00804 
00805     int numbas = 0;
00806     for (int i=0; i<nq; i++) {
00807       numbas += MVT::GetNumberVecs( *Q[i] );
00808     }
00809 
00810     // check size of B
00811     TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00812                         "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
00813     // check size of X and MX
00814     TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00815                         "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
00816     // check size of X w.r.t. MX 
00817     TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00818                         "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
00819     // check feasibility
00820     //TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00821     //                    "Belos::DGKSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
00822 
00823     // Some flags for checking dependency returns from the internal orthogonalization methods
00824     bool dep_flg = false;
00825 
00826     // Make a temporary copy of X and MX, just in case a block dependency is detected.
00827     Teuchos::RCP<MV> tmpX, tmpMX;
00828     tmpX = MVT::CloneCopy(X);
00829     if (this->_hasOp) {
00830       tmpMX = MVT::CloneCopy(*MX);
00831     }
00832 
00833     // Use the cheaper block orthogonalization.
00834     dep_flg = blkOrtho( X, MX, C, Q );
00835 
00836     // If a dependency has been detected in this block, then perform
00837     // the more expensive single-vector orthogonalization.
00838     if (dep_flg) {
00839       rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00840 
00841       // Copy tmpX back into X.
00842       MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00843       if (this->_hasOp) {
00844   MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00845       }
00846     } 
00847     else {
00848       // There is no dependency, so orthonormalize new block X
00849       rank = findBasis( X, MX, B, false );
00850       if (rank < xc) {
00851   // A dependency was found during orthonormalization of X,
00852   // rerun orthogonalization using more expensive single-vector orthogonalization.
00853   rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00854 
00855   // Copy tmpX back into X.
00856   MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00857   if (this->_hasOp) {
00858     MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00859   }
00860       }    
00861     }
00862 
00863     // this should not raise an std::exception; but our post-conditions oblige us to check
00864     TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00865                         "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
00866 
00867     // Return the rank of X.
00868     return rank;
00869   }
00870 
00871 
00872 
00874   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00875   template<class ScalarType, class MV, class OP>
00876   int DGKSOrthoManager<ScalarType, MV, OP>::normalize(
00877                                 MV &X, Teuchos::RCP<MV> MX, 
00878                                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00879 
00880 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00881     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00882 #endif
00883 
00884     // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
00885     return findBasis(X, MX, B, true);
00886   }
00887 
00888 
00889 
00891   template<class ScalarType, class MV, class OP>
00892   void DGKSOrthoManager<ScalarType, MV, OP>::project(
00893                           MV &X, Teuchos::RCP<MV> MX, 
00894                           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00895                           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00896     // For the inner product defined by the operator Op or the identity (Op == 0)
00897     //   -> Orthogonalize X against each Q[i]
00898     // Modify MX accordingly
00899     //
00900     // Note that when Op is 0, MX is not referenced
00901     //
00902     // Parameter variables
00903     //
00904     // X  : Vectors to be transformed
00905     //
00906     // MX : Image of the block vector X by the mass matrix
00907     //
00908     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
00909     //
00910 
00911 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00912     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00913 #endif
00914 
00915     int xc = MVT::GetNumberVecs( X );
00916     ptrdiff_t xr = MVText::GetGlobalLength( X );
00917     int nq = Q.size();
00918     std::vector<int> qcs(nq);
00919     // short-circuit
00920     if (nq == 0 || xc == 0 || xr == 0) {
00921       return;
00922     }
00923     ptrdiff_t qr = MVText::GetGlobalLength ( *Q[0] );
00924     // if we don't have enough C, expand it with null references
00925     // if we have too many, resize to throw away the latter ones
00926     // if we have exactly as many as we have Q, this call has no effect
00927     C.resize(nq);
00928 
00929 
00930     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00931     if (this->_hasOp) {
00932       if (MX == Teuchos::null) {
00933         // we need to allocate space for MX
00934         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00935         OPT::Apply(*(this->_Op),X,*MX);
00936       }
00937     }
00938     else {
00939       // Op == I  -->  MX = X (ignore it if the user passed it in)
00940       MX = Teuchos::rcp( &X, false );
00941     }
00942     int mxc = MVT::GetNumberVecs( *MX );
00943     ptrdiff_t mxr = MVText::GetGlobalLength( *MX );
00944 
00945     // check size of X and Q w.r.t. common sense
00946     TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00947                         "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
00948     // check size of X w.r.t. MX and Q
00949     TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
00950                         "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
00951 
00952     // tally up size of all Q and check/allocate C
00953     int baslen = 0;
00954     for (int i=0; i<nq; i++) {
00955       TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument, 
00956                           "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
00957       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00958       TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
00959                           "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
00960       baslen += qcs[i];
00961 
00962       // check size of C[i]
00963       if ( C[i] == Teuchos::null ) {
00964         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00965       }
00966       else {
00967         TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
00968                            "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
00969       }
00970     }
00971 
00972     // Use the cheaper block orthogonalization, don't check for rank deficiency.
00973     blkOrtho( X, MX, C, Q );
00974 
00975   }  
00976 
00978   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
00979   // the rank is numvectors(X)
00980   template<class ScalarType, class MV, class OP>
00981   int DGKSOrthoManager<ScalarType, MV, OP>::findBasis(
00982                   MV &X, Teuchos::RCP<MV> MX, 
00983                   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00984                   bool completeBasis, int howMany ) const {
00985     // For the inner product defined by the operator Op or the identity (Op == 0)
00986     //   -> Orthonormalize X 
00987     // Modify MX accordingly
00988     //
00989     // Note that when Op is 0, MX is not referenced
00990     //
00991     // Parameter variables
00992     //
00993     // X  : Vectors to be orthonormalized
00994     //
00995     // MX : Image of the multivector X under the operator Op
00996     //
00997     // Op  : Pointer to the operator for the inner product
00998     //
00999     //
01000 
01001     const ScalarType ONE  = SCT::one();
01002     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
01003 
01004     int xc = MVT::GetNumberVecs( X );
01005     ptrdiff_t xr = MVText::GetGlobalLength( X );
01006 
01007     if (howMany == -1) {
01008       howMany = xc;
01009     }
01010 
01011     /*******************************************************
01012      *  If _hasOp == false, we will not reference MX below *
01013      *******************************************************/
01014 
01015     // if Op==null, MX == X (via pointer)
01016     // Otherwise, either the user passed in MX or we will allocated and compute it
01017     if (this->_hasOp) {
01018       if (MX == Teuchos::null) {
01019         // we need to allocate space for MX
01020         MX = MVT::Clone(X,xc);
01021         OPT::Apply(*(this->_Op),X,*MX);
01022       }
01023     }
01024 
01025     /* if the user doesn't want to store the coefficienets, 
01026      * allocate some local memory for them 
01027      */
01028     if ( B == Teuchos::null ) {
01029       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
01030     }
01031 
01032     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
01033     ptrdiff_t mxr = (this->_hasOp) ? MVText::GetGlobalLength( *MX ) : xr;
01034 
01035     // check size of C, B
01036     TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
01037                         "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
01038     TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
01039                         "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
01040     TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
01041                         "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
01042     TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument, 
01043                         "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
01044     TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
01045                         "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
01046 
01047     /* xstart is which column we are starting the process with, based on howMany
01048      * columns before xstart are assumed to be Op-orthonormal already
01049      */
01050     int xstart = xc - howMany;
01051 
01052     for (int j = xstart; j < xc; j++) {
01053 
01054       // numX is 
01055       // * number of currently orthonormal columns of X
01056       // * the index of the current column of X
01057       int numX = j;
01058       bool addVec = false;
01059 
01060       // Get a view of the vector currently being worked on.
01061       std::vector<int> index(1);
01062       index[0] = numX;
01063       Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
01064       Teuchos::RCP<MV> MXj;
01065       if ((this->_hasOp)) {
01066         // MXj is a view of the current vector in MX
01067         MXj = MVT::CloneViewNonConst( *MX, index );
01068       }
01069       else {
01070         // MXj is a pointer to Xj, and MUST NOT be modified
01071         MXj = Xj;
01072       }
01073 
01074       // Get a view of the previous vectors.
01075       std::vector<int> prev_idx( numX );
01076       Teuchos::RCP<const MV> prevX, prevMX;
01077 
01078       if (numX > 0) {
01079         for (int i=0; i<numX; i++) {
01080           prev_idx[i] = i;
01081         }
01082         prevX = MVT::CloneView( X, prev_idx );
01083         if (this->_hasOp) {
01084           prevMX = MVT::CloneView( *MX, prev_idx );
01085         }
01086       } 
01087 
01088       // Make storage for these Gram-Schmidt iterations.
01089       Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
01090       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01091       //
01092       // Save old MXj vector and compute Op-norm
01093       //
01094       Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
01095       MVT::MvDot( *Xj, *MXj, oldDot );
01096       // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
01097       TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
01098         "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
01099 
01100       if (numX > 0) {
01101   // Apply the first step of Gram-Schmidt
01102   
01103   // product <- prevX^T MXj
01104   MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,product);
01105   
01106   // Xj <- Xj - prevX prevX^T MXj   
01107   //     = Xj - prevX product
01108   MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
01109   
01110   // Update MXj
01111   if (this->_hasOp) {
01112     // MXj <- Op*Xj_new
01113     //      = Op*(Xj_old - prevX prevX^T MXj)
01114     //      = MXj - prevMX product
01115     MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
01116   }
01117   
01118   // Compute new Op-norm
01119   MVT::MvDot( *Xj, *MXj, newDot );
01120   
01121   // Check if a correction is needed.
01122   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(dep_tol_*oldDot[0]) ) {
01123     // Apply the second step of Gram-Schmidt
01124     // This is the same as above
01125     Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
01126     
01127     MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,P2);
01128     product += P2;
01129     MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
01130     if ((this->_hasOp)) {
01131       MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
01132     }
01133   } // if (newDot[0] < dep_tol_*oldDot[0])
01134   
01135       } // if (numX > 0)
01136 
01137         // Compute Op-norm with old MXj
01138       MVT::MvDot( *Xj, *oldMXj, newDot );
01139 
01140       // Check to see if the new vector is dependent.
01141       if (completeBasis) {
01142   //
01143   // We need a complete basis, so add random vectors if necessary
01144   //
01145   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
01146     
01147     // Add a random vector and orthogonalize it against previous vectors in block.
01148     addVec = true;
01149 #ifdef ORTHO_DEBUG
01150     std::cout << "Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
01151 #endif
01152     //
01153     Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01154     Teuchos::RCP<MV> tempMXj;
01155     MVT::MvRandom( *tempXj );
01156     if (this->_hasOp) {
01157       tempMXj = MVT::Clone( X, 1 );
01158       OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01159     } 
01160     else {
01161       tempMXj = tempXj;
01162     }
01163     MVT::MvDot( *tempXj, *tempMXj, oldDot );
01164     //
01165     for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
01166       MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,product);
01167       MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
01168       if (this->_hasOp) {
01169         MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
01170       }
01171     }
01172     // Compute new Op-norm
01173     MVT::MvDot( *tempXj, *tempMXj, newDot );
01174     //
01175     if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01176       // Copy vector into current column of _basisvecs
01177       MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
01178       if (this->_hasOp) {
01179         MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
01180       }
01181     }
01182     else {
01183       return numX;
01184     } 
01185   }
01186       }
01187       else {
01188   //
01189   // We only need to detect dependencies.
01190   //
01191   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
01192     return numX;
01193   }
01194       }
01195       
01196       // If we haven't left this method yet, then we can normalize the new vector Xj.
01197       // Normalize Xj.
01198       // Xj <- Xj / std::sqrt(newDot)
01199       ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01200 
01201       if (SCT::magnitude(diag) > ZERO) {      
01202         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01203         if (this->_hasOp) {
01204     // Update MXj.
01205     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01206         }
01207       }
01208 
01209       // If we've added a random vector, enter a zero in the j'th diagonal element.
01210       if (addVec) {
01211   (*B)(j,j) = ZERO;
01212       }
01213       else {
01214   (*B)(j,j) = diag;
01215       }
01216 
01217       // Save the coefficients, if we are working on the original vector and not a randomly generated one
01218       if (!addVec) {
01219   for (int i=0; i<numX; i++) {
01220     (*B)(i,j) = product(i,0);
01221   }
01222       }
01223 
01224     } // for (j = 0; j < xc; ++j)
01225 
01226     return xc;
01227   }
01228 
01230   // Routine to compute the block orthogonalization
01231   template<class ScalarType, class MV, class OP>
01232   bool 
01233   DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
01234                Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01235                Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
01236   {
01237     int nq = Q.size();
01238     int xc = MVT::GetNumberVecs( X );
01239     bool dep_flg = false;
01240     const ScalarType ONE  = SCT::one();
01241 
01242     std::vector<int> qcs( nq );
01243     for (int i=0; i<nq; i++) {
01244       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01245     }
01246 
01247     // Perform the Gram-Schmidt transformation for a block of vectors
01248     
01249     // Compute the initial Op-norms
01250     std::vector<ScalarType> oldDot( xc );
01251     MVT::MvDot( X, *MX, oldDot );
01252 
01253     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
01254     // Define the product Q^T * (Op*X)
01255     for (int i=0; i<nq; i++) {
01256       // Multiply Q' with MX
01257       MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,*C[i]);
01258       // Multiply by Q and subtract the result in X
01259       MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
01260 
01261       // Update MX, with the least number of applications of Op as possible
01262       if (this->_hasOp) {
01263         if (xc <= qcs[i]) {
01264           OPT::Apply( *(this->_Op), X, *MX);
01265         }
01266         else {
01267           // this will possibly be used again below; don't delete it
01268           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01269           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01270           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01271         }
01272       }
01273     }
01274 
01275     // Do as many steps of classical Gram-Schmidt as required by max_blk_ortho_
01276     for (int j = 1; j < max_blk_ortho_; ++j) {
01277       
01278       for (int i=0; i<nq; i++) {
01279   Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
01280         
01281   // Apply another step of classical Gram-Schmidt
01282   MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,C2);
01283   *C[i] += C2;
01284   MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
01285         
01286   // Update MX, with the least number of applications of Op as possible
01287   if (this->_hasOp) {
01288     if (MQ[i].get()) {
01289       // MQ was allocated and computed above; use it
01290       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01291     }
01292     else if (xc <= qcs[i]) {
01293       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01294       OPT::Apply( *(this->_Op), X, *MX);
01295     }
01296   }
01297       } // for (int i=0; i<nq; i++)
01298     } // for (int j = 0; j < max_blk_ortho; ++j)
01299   
01300     // Compute new Op-norms
01301     std::vector<ScalarType> newDot(xc);
01302     MVT::MvDot( X, *MX, newDot );
01303  
01304     // Check to make sure the new block of vectors are not dependent on previous vectors
01305     for (int i=0; i<xc; i++){
01306       if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
01307   dep_flg = true;
01308   break;
01309       }
01310     } // end for (i=0;...)
01311 
01312     return dep_flg;
01313   }
01314   
01315 
01316   template<class ScalarType, class MV, class OP>
01317   int
01318   DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
01319                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01320                    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
01321                    Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
01322   {
01323     Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
01324 
01325     const ScalarType ONE  = SCT::one();
01326     const ScalarType ZERO  = SCT::zero();
01327     
01328     int nq = Q.size();
01329     int xc = MVT::GetNumberVecs( X );
01330     std::vector<int> indX( 1 );
01331     std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01332 
01333     std::vector<int> qcs( nq );
01334     for (int i=0; i<nq; i++) {
01335       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01336     }
01337 
01338     // Create pointers for the previous vectors of X that have already been orthonormalized.
01339     Teuchos::RCP<const MV> lastQ;
01340     Teuchos::RCP<MV> Xj, MXj;
01341     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
01342 
01343     // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
01344     for (int j=0; j<xc; j++) {
01345       
01346       bool dep_flg = false;
01347       
01348       // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
01349       if (j > 0) {
01350   std::vector<int> index( j );
01351   for (int ind=0; ind<j; ind++) {
01352     index[ind] = ind;
01353   }
01354   lastQ = MVT::CloneView( X, index );
01355 
01356   // Add these views to the Q and C arrays.
01357   Q.push_back( lastQ );
01358   C.push_back( B );
01359   qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
01360       }
01361       
01362       // Get a view of the current vector in X to orthogonalize.
01363       indX[0] = j;
01364       Xj = MVT::CloneViewNonConst( X, indX );
01365       if (this->_hasOp) {
01366   MXj = MVT::CloneViewNonConst( *MX, indX );
01367       }
01368       else {
01369   MXj = Xj;
01370       }
01371       
01372       // Compute the initial Op-norms
01373       MVT::MvDot( *Xj, *MXj, oldDot );
01374       
01375       Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
01376       // Define the product Q^T * (Op*X)
01377       for (int i=0; i<Q.size(); i++) {
01378 
01379   // Get a view of the current serial dense matrix
01380   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01381 
01382   // Multiply Q' with MX
01383   MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,tempC);
01384   // Multiply by Q and subtract the result in Xj
01385   MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
01386   
01387   // Update MXj, with the least number of applications of Op as possible
01388   if (this->_hasOp) {
01389     if (xc <= qcs[i]) {
01390       OPT::Apply( *(this->_Op), *Xj, *MXj);
01391     }
01392     else {
01393       // this will possibly be used again below; don't delete it
01394       MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01395       OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01396       MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
01397     }
01398   }
01399       }
01400       
01401       // Compute the Op-norms
01402       MVT::MvDot( *Xj, *MXj, newDot );
01403       
01404       // Do one step of classical Gram-Schmidt orthogonalization 
01405       // with a second correction step if needed.
01406       
01407       if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
01408   
01409   for (int i=0; i<Q.size(); i++) {
01410     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01411     Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
01412     
01413     // Apply another step of classical Gram-Schmidt
01414     MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,C2);
01415     tempC += C2;
01416     MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
01417     
01418     // Update MXj, with the least number of applications of Op as possible
01419     if (this->_hasOp) {
01420       if (MQ[i].get()) {
01421         // MQ was allocated and computed above; use it
01422         MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
01423       }
01424       else if (xc <= qcs[i]) {
01425         // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01426         OPT::Apply( *(this->_Op), *Xj, *MXj);
01427       }
01428     }
01429   } // for (int i=0; i<Q.size(); i++)
01430   
01431   // Compute the Op-norms after the correction step.
01432   MVT::MvDot( *Xj, *MXj, newDot );
01433   
01434       } // if ()
01435       
01436       // Check for linear dependence.
01437       if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
01438   dep_flg = true;
01439       }
01440       
01441       // Normalize the new vector if it's not dependent
01442       if (!dep_flg) {
01443   ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01444   
01445   MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01446   if (this->_hasOp) {
01447     // Update MXj.
01448     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01449   }
01450   
01451   // Enter value on diagonal of B.
01452   (*B)(j,j) = diag;
01453       }
01454       else {
01455   // Create a random vector and orthogonalize it against all previous columns of Q.
01456   Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01457   Teuchos::RCP<MV> tempMXj;
01458   MVT::MvRandom( *tempXj );
01459   if (this->_hasOp) {
01460     tempMXj = MVT::Clone( X, 1 );
01461     OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01462   } 
01463   else {
01464     tempMXj = tempXj;
01465   }
01466   MVT::MvDot( *tempXj, *tempMXj, oldDot );
01467   //
01468   for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
01469     
01470     for (int i=0; i<Q.size(); i++) {
01471       Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
01472       
01473       // Apply another step of classical Gram-Schmidt
01474       MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*tempXj,tempMXj,product);
01475       MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
01476       
01477       // Update MXj, with the least number of applications of Op as possible
01478       if (this->_hasOp) {
01479         if (MQ[i].get()) {
01480     // MQ was allocated and computed above; use it
01481     MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
01482         }
01483         else if (xc <= qcs[i]) {
01484     // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01485     OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
01486         }
01487       }
01488     } // for (int i=0; i<nq; i++)
01489     
01490   }
01491   
01492   // Compute the Op-norms after the correction step.
01493   MVT::MvDot( *tempXj, *tempMXj, newDot );
01494   
01495   // Copy vector into current column of Xj
01496   if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01497     ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01498     
01499     // Enter value on diagonal of B.
01500     (*B)(j,j) = ZERO;
01501 
01502     // Copy vector into current column of _basisvecs
01503     MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
01504     if (this->_hasOp) {
01505       MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
01506     }
01507   }
01508   else {
01509     return j;
01510   } 
01511       } // if (!dep_flg)
01512 
01513       // Remove the vectors from array
01514       if (j > 0) {
01515   Q.resize( nq );
01516   C.resize( nq );
01517   qcs.resize( nq );
01518       }
01519 
01520     } // for (int j=0; j<xc; j++)
01521 
01522     return xc;
01523   }
01524 
01525 } // namespace Belos
01526 
01527 #endif // BELOS_DGKS_ORTHOMANAGER_HPP
01528 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines