Anasazi Version of the Day
AnasaziBasicOrthoManager.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 
00034 #ifndef ANASAZI_BASIC_ORTHOMANAGER_HPP
00035 #define ANASAZI_BASIC_ORTHOMANAGER_HPP
00036 
00044 #include "AnasaziConfigDefs.hpp"
00045 #include "AnasaziMultiVecTraits.hpp"
00046 #include "AnasaziOperatorTraits.hpp"
00047 #include "AnasaziMatOrthoManager.hpp"
00048 #include "Teuchos_TimeMonitor.hpp"
00049 #include "Teuchos_LAPACK.hpp"
00050 #include "Teuchos_BLAS.hpp"
00051 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00052 #  include <Teuchos_FancyOStream.hpp>
00053 #endif
00054 
00055 namespace Anasazi {
00056 
00057   template<class ScalarType, class MV, class OP>
00058   class BasicOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
00059 
00060   private:
00061     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00062     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
00063     typedef MultiVecTraits<ScalarType,MV>      MVT;
00064     typedef OperatorTraits<ScalarType,MV,OP>   OPT;
00065 
00066   public:
00067 
00069 
00070 
00071     BasicOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null, 
00072                        typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa = 1.41421356 /* sqrt(2) */,
00073                        typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0,
00074                        typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 );
00075 
00076 
00078     ~BasicOrthoManager() {}
00080 
00081 
00083 
00084 
00085 
00124     void projectMat ( 
00125           MV &X, 
00126           Teuchos::Array<Teuchos::RCP<const MV> > Q,
00127           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00128               = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00129           Teuchos::RCP<MV> MX                        = Teuchos::null,
00130           Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00131         ) const;
00132 
00133 
00172     int normalizeMat ( 
00173           MV &X, 
00174           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00175           Teuchos::RCP<MV> MX                                         = Teuchos::null
00176         ) const;
00177 
00178 
00238     int projectAndNormalizeMat ( 
00239           MV &X, 
00240           Teuchos::Array<Teuchos::RCP<const MV> >  Q,
00241           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00242               = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00243           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B                  = Teuchos::null, 
00244           Teuchos::RCP<MV> MX                        = Teuchos::null,
00245           Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00246         ) const;
00247 
00249 
00251 
00252 
00257     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00258     orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
00259 
00264     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00265     orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const;
00266 
00268 
00270 
00271 
00273     void setKappa( typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa ) { kappa_ = kappa; }
00274 
00276     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType getKappa() const { return kappa_; } 
00277 
00279 
00280   private:
00282     MagnitudeType kappa_;
00283     MagnitudeType eps_;
00284     MagnitudeType tol_;
00285   
00286     // ! Routine to find an orthonormal basis
00287     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00288                          Teuchos::SerialDenseMatrix<int,ScalarType> &B, 
00289                          bool completeBasis, int howMany = -1 ) const;
00290 
00291     //
00292     // Internal timers
00293     //
00294     Teuchos::RCP<Teuchos::Time> timerReortho_;
00295     
00296   };
00297 
00298 
00300   // Constructor
00301   template<class ScalarType, class MV, class OP>
00302   BasicOrthoManager<ScalarType,MV,OP>::BasicOrthoManager( Teuchos::RCP<const OP> Op,
00303                                                           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa,
00304                                                           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps,
00305                                                           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol ) : 
00306     MatOrthoManager<ScalarType,MV,OP>(Op), kappa_(kappa), eps_(eps), tol_(tol)
00307 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00308     , timerReortho_(Teuchos::TimeMonitor::getNewTimer("Anasazi::BasicOrthoManager::Re-orthogonalization"))
00309 #endif
00310   {
00311     TEUCHOS_TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument,
00312         "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"eps\" must be non-negative.");
00313     if (eps_ == 0) {
00314       Teuchos::LAPACK<int,MagnitudeType> lapack;
00315       eps_ = lapack.LAMCH('E');
00316       eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.75);
00317     }
00318     TEUCHOS_TEST_FOR_EXCEPTION(
00319         tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()),
00320         std::invalid_argument,
00321         "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"tol\" must be in [0,1].");
00322   }
00323 
00324 
00325 
00327   // Compute the distance from orthonormality
00328   template<class ScalarType, class MV, class OP>
00329   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00330   BasicOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const {
00331     const ScalarType ONE = SCT::one();
00332     int rank = MVT::GetNumberVecs(X);
00333     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00334     MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X,X,xTx,MX,MX);
00335     for (int i=0; i<rank; i++) {
00336       xTx(i,i) -= ONE;
00337     }
00338     return xTx.normFrobenius();
00339   }
00340 
00341 
00342 
00344   // Compute the distance from orthogonality
00345   template<class ScalarType, class MV, class OP>
00346   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00347   BasicOrthoManager<ScalarType,MV,OP>::orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const {
00348     int r1 = MVT::GetNumberVecs(X1);
00349     int r2  = MVT::GetNumberVecs(X2);
00350     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
00351     MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X1,X2,xTx,MX1,MX2);
00352     return xTx.normFrobenius();
00353   }
00354 
00355 
00356 
00358   template<class ScalarType, class MV, class OP>
00359   void BasicOrthoManager<ScalarType, MV, OP>::projectMat(
00360           MV &X, 
00361           Teuchos::Array<Teuchos::RCP<const MV> > Q,
00362           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00363           Teuchos::RCP<MV> MX,
00364           Teuchos::Array<Teuchos::RCP<const MV> > MQ
00365       ) const {
00366     // For the inner product defined by the operator Op or the identity (Op == 0)
00367     //   -> Orthogonalize X against each Q[i]
00368     // Modify MX accordingly
00369     //
00370     // Note that when Op is 0, MX is not referenced
00371     //
00372     // Parameter variables
00373     //
00374     // X  : Vectors to be transformed
00375     //
00376     // MX : Image of the block vector X by the mass matrix
00377     //
00378     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
00379     //
00380 
00381 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00382     // Get a FancyOStream from out_arg or create a new one ...
00383     Teuchos::RCP<Teuchos::FancyOStream>
00384       out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
00385     out->setShowAllFrontMatter(false).setShowProcRank(true);
00386     *out << "Entering Anasazi::BasicOrthoManager::projectMat(...)\n";
00387 #endif
00388 
00389     ScalarType ONE  = SCT::one();
00390 
00391     int xc = MVT::GetNumberVecs( X );
00392     int xr = MVT::GetVecLength( X );
00393     int nq = Q.length();
00394     std::vector<int> qcs(nq);
00395     // short-circuit
00396     if (nq == 0 || xc == 0 || xr == 0) {
00397 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00398       *out << "Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
00399 #endif
00400       return;
00401     }
00402     int qr = MVT::GetVecLength ( *Q[0] );
00403     // if we don't have enough C, expand it with null references
00404     // if we have too many, resize to throw away the latter ones
00405     // if we have exactly as many as we have Q, this call has no effect
00406     C.resize(nq);
00407 
00408 
00409     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00410     if (this->_hasOp) {
00411 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00412       *out << "Allocating MX...\n";
00413 #endif
00414       if (MX == Teuchos::null) {
00415         // we need to allocate space for MX
00416         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00417         OPT::Apply(*(this->_Op),X,*MX);
00418         this->_OpCounter += MVT::GetNumberVecs(X);
00419       }
00420     }
00421     else {
00422       // Op == I  -->  MX = X (ignore it if the user passed it in)
00423       MX = Teuchos::rcpFromRef(X);
00424     }
00425     int mxc = MVT::GetNumberVecs( *MX );
00426     int mxr = MVT::GetVecLength( *MX );
00427 
00428     // check size of X and Q w.r.t. common sense
00429     TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00430                         "Anasazi::BasicOrthoManager::projectMat(): MVT returned negative dimensions for X,MX" );
00431     // check size of X w.r.t. MX and Q
00432     TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
00433                         "Anasazi::BasicOrthoManager::projectMat(): Size of X not consistent with MX,Q" );
00434 
00435     // tally up size of all Q and check/allocate C
00436     int baslen = 0;
00437     for (int i=0; i<nq; i++) {
00438       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 
00439                           "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" );
00440       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00441       TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
00442                           "Anasazi::BasicOrthoManager::projectMat(): Q has less rows than columns" );
00443       baslen += qcs[i];
00444 
00445       // check size of C[i]
00446       if ( C[i] == Teuchos::null ) {
00447         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00448       }
00449       else {
00450         TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
00451                            "Anasazi::BasicOrthoManager::projectMat(): Size of Q not consistent with size of C" );
00452       }
00453     }
00454 
00455     // Perform the Gram-Schmidt transformation for a block of vectors
00456 
00457     // Compute the initial Op-norms
00458     std::vector<ScalarType> oldDot( xc );
00459     MVT::MvDot( X, *MX, oldDot );
00460 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00461     *out << "oldDot = { ";
00462     std::copy(oldDot.begin(), oldDot.end(), std::ostream_iterator<ScalarType>(*out, " "));
00463     *out << "}\n";
00464 #endif
00465 
00466     MQ.resize(nq);
00467     // Define the product Q^T * (Op*X)
00468     for (int i=0; i<nq; i++) {
00469       // Multiply Q' with MX
00470       MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Q[i],X,*C[i],MQ[i],MX);
00471       // Multiply by Q and subtract the result in X
00472 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00473       *out << "Applying projector P_Q[" << i << "]...\n";
00474 #endif
00475       MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
00476 
00477       // Update MX, with the least number of applications of Op as possible
00478       // Update MX. If we have MQ, use it. Otherwise, just multiply by Op
00479       if (this->_hasOp) {
00480         if (MQ[i] == Teuchos::null) {
00481 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00482           *out << "Updating MX via M*X...\n";
00483 #endif
00484           OPT::Apply( *(this->_Op), X, *MX);
00485           this->_OpCounter += MVT::GetNumberVecs(X);
00486         }
00487         else {
00488 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00489           *out << "Updating MX via M*Q...\n";
00490 #endif
00491           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
00492         }
00493       }
00494     }
00495 
00496     // Compute new Op-norms
00497     std::vector<ScalarType> newDot(xc);
00498     MVT::MvDot( X, *MX, newDot );
00499 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00500     *out << "newDot = { ";
00501     std::copy(newDot.begin(), newDot.end(), std::ostream_iterator<ScalarType>(*out, " "));
00502     *out << "}\n";
00503 #endif
00504 
00505     // determine (individually) whether to do another step of classical Gram-Schmidt
00506     for (int j = 0; j < xc; ++j) {
00507       
00508       if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
00509 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00510         *out << "kappa_*newDot[" <<j<< "] == " << kappa_*newDot[j] << "... another step of Gram-Schmidt.\n";
00511 #endif
00512 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00513         Teuchos::TimeMonitor lcltimer( *timerReortho_ );
00514 #endif
00515         for (int i=0; i<nq; i++) {
00516           Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
00517           
00518           // Apply another step of classical Gram-Schmidt
00519           MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Q[i],X,C2,MQ[i],MX);
00520           *C[i] += C2;
00521 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00522           *out << "Applying projector P_Q[" << i << "]...\n";
00523 #endif
00524           MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
00525           
00526           // Update MX as above
00527           if (this->_hasOp) {
00528             if (MQ[i] == Teuchos::null) {
00529 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00530               *out << "Updating MX via M*X...\n";
00531 #endif
00532               OPT::Apply( *(this->_Op), X, *MX);
00533               this->_OpCounter += MVT::GetNumberVecs(X);
00534             }
00535             else {
00536 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00537               *out << "Updating MX via M*Q...\n";
00538 #endif
00539               MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
00540             }
00541           }
00542         }
00543         break;
00544       } // if (kappa_*newDot[j] < oldDot[j])
00545     } // for (int j = 0; j < xc; ++j)
00546 
00547 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00548     *out << "Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
00549 #endif
00550   }
00551 
00552 
00553 
00555   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00556   template<class ScalarType, class MV, class OP>
00557   int BasicOrthoManager<ScalarType, MV, OP>::normalizeMat(
00558         MV &X, 
00559         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00560         Teuchos::RCP<MV> MX) const {
00561     // call findBasis(), with the instruction to try to generate a basis of rank numvecs(X)
00562     // findBasis() requires MX
00563 
00564     int xc = MVT::GetNumberVecs(X);
00565     int xr = MVT::GetVecLength(X);
00566 
00567     // if Op==null, MX == X (via pointer)
00568     // Otherwise, either the user passed in MX or we will allocated and compute it
00569     if (this->_hasOp) {
00570       if (MX == Teuchos::null) {
00571         // we need to allocate space for MX
00572         MX = MVT::Clone(X,xc);
00573         OPT::Apply(*(this->_Op),X,*MX);
00574         this->_OpCounter += MVT::GetNumberVecs(X);
00575       }
00576     }
00577 
00578     // if the user doesn't want to store the coefficients, 
00579     // allocate some local memory for them 
00580     if ( B == Teuchos::null ) {
00581       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00582     }
00583 
00584     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00585     int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX )  : xr;
00586 
00587     // check size of C, B
00588     TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
00589                         "Anasazi::BasicOrthoManager::normalizeMat(): X must be non-empty" );
00590     TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00591                         "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
00592     TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
00593                         "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
00594     TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 
00595                         "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
00596 
00597     return findBasis(X, MX, *B, true );
00598   }
00599 
00600 
00601 
00603   // Find an Op-orthonormal basis for span(X) - span(W)
00604   template<class ScalarType, class MV, class OP>
00605   int BasicOrthoManager<ScalarType, MV, OP>::projectAndNormalizeMat(
00606           MV &X, 
00607           Teuchos::Array<Teuchos::RCP<const MV> >  Q,
00608           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00609           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00610           Teuchos::RCP<MV> MX,
00611           Teuchos::Array<Teuchos::RCP<const MV> > MQ
00612       ) const {
00613 
00614 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00615     // Get a FancyOStream from out_arg or create a new one ...
00616     Teuchos::RCP<Teuchos::FancyOStream>
00617       out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
00618     out->setShowAllFrontMatter(false).setShowProcRank(true);
00619     *out << "Entering Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
00620 #endif
00621 
00622     int nq = Q.length();
00623     int xc = MVT::GetNumberVecs( X );
00624     int xr = MVT::GetVecLength( X );
00625     int rank;
00626 
00627     /* if the user doesn't want to store the coefficients, 
00628      * allocate some local memory for them 
00629      */
00630     if ( B == Teuchos::null ) {
00631       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00632     }
00633 
00634     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00635     if (this->_hasOp) {
00636       if (MX == Teuchos::null) {
00637         // we need to allocate space for MX
00638 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00639         *out << "Allocating MX...\n";
00640 #endif
00641         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00642         OPT::Apply(*(this->_Op),X,*MX);
00643         this->_OpCounter += MVT::GetNumberVecs(X);
00644       }
00645     }
00646     else {
00647       // Op == I  -->  MX = X (ignore it if the user passed it in)
00648       MX = Teuchos::rcpFromRef(X);
00649     }
00650 
00651     int mxc = MVT::GetNumberVecs( *MX );
00652     int mxr = MVT::GetVecLength( *MX );
00653 
00654     TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): X must be non-empty" );
00655 
00656     int numbas = 0;
00657     for (int i=0; i<nq; i++) {
00658       numbas += MVT::GetNumberVecs( *Q[i] );
00659     }
00660 
00661     // check size of B
00662     TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00663                         "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of B" );
00664     // check size of X and MX
00665     TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00666                         "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): MVT returned negative dimensions for X,MX" );
00667     // check size of X w.r.t. MX 
00668     TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00669                         "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of MX" );
00670     // check feasibility
00671     TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00672                         "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Orthogonality constraints not feasible" );
00673 
00674     // orthogonalize all of X against Q
00675 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00676     *out << "Orthogonalizing X against Q...\n";
00677 #endif
00678     projectMat(X,Q,C,MX,MQ);
00679 
00680     Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(xc,1);
00681     // start working
00682     rank = 0;
00683     int numTries = 10;   // each vector in X gets 10 random chances to escape degeneracy
00684     int oldrank = -1;
00685     do {
00686       int curxsize = xc - rank;
00687 
00688       // orthonormalize X, but quit if it is rank deficient
00689       // we can't let findBasis generated random vectors to complete the basis,
00690       // because it doesn't know about Q; we will do this ourselves below
00691 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00692       *out << "Attempting to find orthonormal basis for X...\n";
00693 #endif
00694       rank = findBasis(X,MX,*B,false,curxsize);
00695 
00696       if (oldrank != -1 && rank != oldrank) {
00697         // we had previously stopped before, after operating on vector oldrank
00698         // we saved its coefficients, augmented it with a random vector, and
00699         // then called findBasis() again, which proceeded to add vector oldrank
00700         // to the basis. 
00701         // now, restore the saved coefficients into B
00702         for (int i=0; i<xc; i++) {
00703           (*B)(i,oldrank) = oldCoeff(i,0);
00704         }
00705       }
00706 
00707       if (rank < xc) {
00708         if (rank != oldrank) {
00709           // we quit on this vector and will augment it with random below
00710           // this is the first time that we have quit on this vector
00711           // therefor, (*B)(:,rank) contains the actual coefficients of the 
00712           // input vectors with respect to the previous vectors in the basis
00713           // save these values, as (*B)(:,rank) will be overwritten by our next
00714           // call to findBasis()
00715           // we will restore it after we are done working on this vector
00716           for (int i=0; i<xc; i++) {
00717             oldCoeff(i,0) = (*B)(i,rank);
00718           }
00719         }
00720       }
00721 
00722       if (rank == xc) {
00723         // we are done
00724 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00725         *out << "Finished computing basis.\n";
00726 #endif
00727         break;
00728       }
00729       else {
00730         TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank, OrthoError,   
00731                             "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): basis lost rank; this shouldn't happen");
00732 
00733         if (rank != oldrank) {
00734           // we added a vector to the basis; reset the chance counter
00735           numTries = 10;
00736           // store old rank
00737           oldrank = rank;
00738         }
00739         else {
00740           // has this vector run out of chances to escape degeneracy?
00741           if (numTries <= 0) {
00742             break;
00743           }
00744         }
00745         // use one of this vector's chances
00746         numTries--;
00747 
00748         // randomize troubled direction
00749 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00750         *out << "Randomizing X[" << rank << "]...\n";
00751 #endif
00752         Teuchos::RCP<MV> curX, curMX;
00753         std::vector<int> ind(1);
00754         ind[0] = rank;
00755         curX  = MVT::CloneViewNonConst(X,ind);
00756         MVT::MvRandom(*curX);
00757         if (this->_hasOp) {
00758 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00759           *out << "Applying operator to random vector.\n";
00760 #endif
00761           curMX = MVT::CloneViewNonConst(*MX,ind);
00762           OPT::Apply( *(this->_Op), *curX, *curMX );
00763           this->_OpCounter += MVT::GetNumberVecs(*curX);
00764         }
00765 
00766         // orthogonalize against Q
00767         // if !this->_hasOp, the curMX will be ignored.
00768         // we don't care about these coefficients
00769         // on the contrary, we need to preserve the previous coeffs
00770         {
00771           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
00772           projectMat(*curX,Q,dummyC,curMX,MQ);
00773         }
00774       }
00775     } while (1);
00776 
00777     // this should never raise an exception; but our post-conditions oblige us to check
00778     TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00779                         "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Debug error in rank variable." );
00780 
00781 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00782     *out << "Leaving Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
00783 #endif
00784 
00785     return rank;
00786   }
00787 
00788 
00789 
00791   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
00792   // the rank is numvectors(X)
00793   template<class ScalarType, class MV, class OP>
00794   int BasicOrthoManager<ScalarType, MV, OP>::findBasis(
00795                 MV &X, Teuchos::RCP<MV> MX, 
00796                 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
00797                 bool completeBasis, int howMany ) const {
00798 
00799     // For the inner product defined by the operator Op or the identity (Op == 0)
00800     //   -> Orthonormalize X 
00801     // Modify MX accordingly
00802     //
00803     // Note that when Op is 0, MX is not referenced
00804     //
00805     // Parameter variables
00806     //
00807     // X  : Vectors to be orthonormalized
00808     //
00809     // MX : Image of the multivector X under the operator Op
00810     //
00811     // Op  : Pointer to the operator for the inner product
00812     //
00813 
00814 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00815     // Get a FancyOStream from out_arg or create a new one ...
00816     Teuchos::RCP<Teuchos::FancyOStream>
00817       out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
00818     out->setShowAllFrontMatter(false).setShowProcRank(true);
00819     *out << "Entering Anasazi::BasicOrthoManager::findBasis(...)\n";
00820 #endif
00821 
00822     const ScalarType ONE  = SCT::one();
00823     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
00824 
00825     int xc = MVT::GetNumberVecs( X );
00826 
00827     if (howMany == -1) {
00828       howMany = xc;
00829     }
00830 
00831     /*******************************************************
00832      *  If _hasOp == false, we will not reference MX below *
00833      *******************************************************/
00834     TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp == true && MX == Teuchos::null, std::logic_error,
00835         "Anasazi::BasicOrthoManager::findBasis(): calling routine did not specify MS.");
00836     TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error, 
00837                         "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" );
00838 
00839     /* xstart is which column we are starting the process with, based on howMany
00840      * columns before xstart are assumed to be Op-orthonormal already
00841      */
00842     int xstart = xc - howMany;
00843 
00844     for (int j = xstart; j < xc; j++) {
00845 
00846       // numX represents the number of currently orthonormal columns of X
00847       int numX = j;
00848       // j represents the index of the current column of X
00849       // these are different interpretations of the same value
00850 
00851       // 
00852       // set the lower triangular part of B to zero
00853       for (int i=j+1; i<xc; ++i) {
00854         B(i,j) = ZERO;
00855       }
00856 
00857       // Get a view of the vector currently being worked on.
00858       std::vector<int> index(1);
00859       index[0] = j;
00860       Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
00861       Teuchos::RCP<MV> MXj;
00862       if ((this->_hasOp)) {
00863         // MXj is a view of the current vector in MX
00864         MXj = MVT::CloneViewNonConst( *MX, index );
00865       }
00866       else {
00867         // MXj is a pointer to Xj, and MUST NOT be modified
00868         MXj = Xj;
00869       }
00870 
00871       // Get a view of the previous vectors.
00872       std::vector<int> prev_idx( numX );
00873       Teuchos::RCP<const MV> prevX, prevMX;
00874 
00875       if (numX > 0) {
00876         for (int i=0; i<numX; ++i) prev_idx[i] = i;
00877         prevX = MVT::CloneViewNonConst( X, prev_idx );
00878         if (this->_hasOp) {
00879           prevMX = MVT::CloneViewNonConst( *MX, prev_idx );
00880         }
00881       } 
00882 
00883       bool rankDef = true;
00884       /* numTrials>0 will denote that the current vector was randomized for the purpose
00885        * of finding a basis vector, and that the coefficients of that vector should
00886        * not be stored in B
00887        */
00888       for (int numTrials = 0; numTrials < 10; numTrials++) {
00889 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00890         *out << "Trial " << numTrials << " for vector " << j << "\n";
00891 #endif
00892 
00893         // Make storage for these Gram-Schmidt iterations.
00894         Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
00895         std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
00896 
00897         //
00898         // Save old MXj vector and compute Op-norm
00899         //
00900         Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
00901         MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,origNorm,MXj);
00902 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00903         *out << "origNorm = " << origNorm[0] << "\n";
00904 #endif
00905 
00906         if (numX > 0) {
00907           // Apply the first step of Gram-Schmidt
00908 
00909           // product <- prevX^T MXj
00910           MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,product,Teuchos::null,MXj);
00911 
00912           // Xj <- Xj - prevX prevX^T MXj   
00913           //     = Xj - prevX product
00914 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 
00915           *out << "Orthogonalizing X[" << j << "]...\n";
00916 #endif
00917           MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
00918 
00919           // Update MXj
00920           if (this->_hasOp) {
00921             // MXj <- Op*Xj_new
00922             //      = Op*(Xj_old - prevX prevX^T MXj)
00923             //      = MXj - prevMX product
00924 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 
00925             *out << "Updating MX[" << j << "]...\n";
00926 #endif
00927             MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
00928           }
00929 
00930           // Compute new Op-norm
00931           MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm,MXj);
00932           MagnitudeType product_norm = product.normOne();
00933           
00934 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00935           *out << "newNorm = " << newNorm[0] << "\n";
00936           *out << "prodoct_norm = " << product_norm << "\n";
00937 #endif
00938 
00939           // Check if a correction is needed.
00940           if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
00941 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00942             if (product_norm/newNorm[0] >= tol_) {
00943               *out << "product_norm/newNorm == " << product_norm/newNorm[0] << "... another step of Gram-Schmidt.\n";
00944             }
00945             else {
00946               *out << "eps*origNorm == " << eps_*origNorm[0] << "... another step of Gram-Schmidt.\n";
00947             }
00948 #endif
00949             // Apply the second step of Gram-Schmidt
00950             // This is the same as above
00951             Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
00952             MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,P2,Teuchos::null,MXj);
00953             product += P2;
00954 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 
00955             *out << "Orthogonalizing X[" << j << "]...\n";
00956 #endif
00957             MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
00958             if ((this->_hasOp)) {
00959 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 
00960               *out << "Updating MX[" << j << "]...\n";
00961 #endif
00962               MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
00963             }
00964             // Compute new Op-norms
00965             MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm2,MXj);
00966             product_norm = P2.normOne();
00967 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00968             *out << "newNorm2 = " << newNorm2[0] << "\n";
00969             *out << "product_norm = " << product_norm << "\n";
00970 #endif
00971             if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
00972               // we don't do another GS, we just set it to zero.
00973 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00974               if (product_norm/newNorm2[0] >= tol_) {
00975                 *out << "product_norm/newNorm2 == " << product_norm/newNorm2[0] << "... setting vector to zero.\n";
00976               }
00977               else if (newNorm[0] < newNorm2[0]) {
00978                 *out << "newNorm2 > newNorm... setting vector to zero.\n";
00979               }
00980               else {
00981                 *out << "eps*origNorm == " << eps_*origNorm[0] << "... setting vector to zero.\n";
00982               }
00983 #endif
00984               MVT::MvInit(*Xj,ZERO);
00985               if ((this->_hasOp)) {
00986 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00987                 *out << "Setting MX[" << j << "] to zero as well.\n";
00988 #endif
00989                 MVT::MvInit(*MXj,ZERO);
00990               }
00991             }
00992           } 
00993         } // if (numX > 0) do GS
00994 
00995         // save the coefficients, if we are working on the original vector and not a randomly generated one
00996         if (numTrials == 0) {
00997           for (int i=0; i<numX; i++) {
00998             B(i,j) = product(i,0);
00999           }
01000         }
01001 
01002         // Check if Xj has any directional information left after the orthogonalization.
01003         MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm,MXj);
01004         if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
01005 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01006           *out << "Normalizing X[" << j << "], norm(X[" << j << "]) = " << newNorm[0] << "\n";
01007 #endif
01008           // Normalize Xj.
01009           // Xj <- Xj / norm
01010           MVT::MvScale( *Xj, ONE/newNorm[0]);
01011           if (this->_hasOp) {
01012 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01013             *out << "Normalizing M*X[" << j << "]...\n";
01014 #endif
01015             // Update MXj.
01016             MVT::MvScale( *MXj, ONE/newNorm[0]);
01017           }
01018 
01019           // save it, if it corresponds to the original vector and not a randomly generated one
01020           if (numTrials == 0) {
01021             B(j,j) = newNorm[0];
01022           }
01023 
01024           // We are not rank deficient in this vector. Move on to the next vector in X.
01025           rankDef = false;
01026           break;
01027         }
01028         else {
01029 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01030           *out << "Not normalizing M*X[" << j << "]...\n";
01031 #endif
01032           // There was nothing left in Xj after orthogonalizing against previous columns in X.
01033           // X is rank deficient.
01034           // reflect this in the coefficients
01035           B(j,j) = ZERO;
01036 
01037           if (completeBasis) {
01038             // Fill it with random information and keep going.
01039 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01040             *out << "Inserting random vector in X[" << j << "]...\n";
01041 #endif
01042             MVT::MvRandom( *Xj );
01043             if (this->_hasOp) {
01044 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01045               *out << "Updating M*X[" << j << "]...\n";
01046 #endif
01047               OPT::Apply( *(this->_Op), *Xj, *MXj );
01048               this->_OpCounter += MVT::GetNumberVecs(*Xj);
01049             }
01050           }
01051           else {
01052             rankDef = true;
01053             break;
01054           }
01055         }
01056       }  // for (numTrials = 0; numTrials < 10; ++numTrials)
01057 
01058       // if rankDef == true, then quit and notify user of rank obtained
01059       if (rankDef == true) {
01060         TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError, 
01061                             "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" );
01062 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01063         *out << "Returning early, rank " << j << " from Anasazi::BasicOrthoManager::findBasis(...)\n";
01064 #endif
01065         return j;
01066       }
01067 
01068     } // for (j = 0; j < xc; ++j)
01069 
01070 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01071     *out << "Returning " << xc << " from Anasazi::BasicOrthoManager::findBasis(...)\n";
01072 #endif
01073     return xc;
01074   }
01075 
01076 } // namespace Anasazi
01077 
01078 #endif // ANASAZI_BASIC_ORTHOMANAGER_HPP
01079 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends