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