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 // #define ORTHO_DEBUG
00045 
00046 #include "AnasaziConfigDefs.hpp"
00047 #include "AnasaziMultiVecTraits.hpp"
00048 #include "AnasaziOperatorTraits.hpp"
00049 #include "AnasaziMatOrthoManager.hpp"
00050 #include "Teuchos_TimeMonitor.hpp"
00051 
00052 namespace Anasazi {
00053 
00054   template<class ScalarType, class MV, class OP>
00055   class BasicOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
00056 
00057   private:
00058     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00059     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
00060     typedef MultiVecTraits<ScalarType,MV>      MVT;
00061     typedef OperatorTraits<ScalarType,MV,OP>   OPT;
00062 
00063   public:
00064     
00066 
00067 
00068     BasicOrthoManager( Teuchos::RefCountPtr<const OP> Op = Teuchos::null, const MagnitudeType kappa = SCT::magnitude(1.5625) );
00069 
00070 
00072     ~BasicOrthoManager() {}
00074 
00075 
00077 
00078 
00080     void setKappa( const MagnitudeType kappa ) { kappa_ = kappa; }
00081 
00083     MagnitudeType getKappa() const { return kappa_; } 
00084 
00086 
00087 
00089 
00090 
00118     void project ( MV &X, Teuchos::RefCountPtr<MV> MX, 
00119                    Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00120                    Teuchos::Array<Teuchos::RefCountPtr<const MV> > Q) const;
00121 
00122 
00125     void project ( MV &X, 
00126                    Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00127                    Teuchos::Array<Teuchos::RefCountPtr<const MV> > Q) const {
00128       project(X,Teuchos::null,C,Q);
00129     }
00130 
00131 
00132  
00157     int normalize ( MV &X, Teuchos::RefCountPtr<MV> MX, 
00158                     Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
00159 
00160 
00163     int normalize ( MV &X, Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00164       return normalize(X,Teuchos::null,B);
00165     }
00166 
00167 
00200     int projectAndNormalize ( MV &X, Teuchos::RefCountPtr<MV> MX, 
00201                               Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00202                               Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00203                               Teuchos::Array<Teuchos::RefCountPtr<const MV> > Q) const;
00204 
00207     int projectAndNormalize ( MV &X, 
00208                               Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00209                               Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00210                               Teuchos::Array<Teuchos::RefCountPtr<const MV> > Q ) const {
00211       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00212     }
00213 
00215 
00217 
00218 
00222     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00223     orthonormError(const MV &X) const {
00224       return orthonormError(X,Teuchos::null);
00225     }
00226 
00231     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00232     orthonormError(const MV &X, Teuchos::RefCountPtr<const MV> MX) const;
00233 
00237     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00238     orthogError(const MV &X1, const MV &X2) const {
00239       return orthogError(X1,Teuchos::null,X2);
00240     }
00241 
00246     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00247     orthogError(const MV &X1, Teuchos::RefCountPtr<const MV> MX1, const MV &X2) const;
00248 
00250 
00251   private:
00252     
00254     MagnitudeType kappa_;
00255   
00256     // ! Routine to find an orthonormal basis for the 
00257     int findBasis(MV &X, Teuchos::RefCountPtr<MV> MX, 
00258                          Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > C, 
00259                          bool completeBasis, int howMany = -1 ) const;
00260 
00261     //
00262     // Internal timers
00263     //
00264     Teuchos::RefCountPtr<Teuchos::Time> timerReortho_;
00265     
00266   };
00267 
00268 
00270   // Constructor
00271   template<class ScalarType, class MV, class OP>
00272   BasicOrthoManager<ScalarType,MV,OP>::BasicOrthoManager( Teuchos::RefCountPtr<const OP> Op,
00273                                                           const MagnitudeType kappa         ) : 
00274     MatOrthoManager<ScalarType,MV,OP>(Op), 
00275     kappa_(kappa), 
00276     timerReortho_(Teuchos::TimeMonitor::getNewTimer("BasicOrthoManager::Re-orthogonalization"))
00277   {}
00278 
00279 
00281   // Compute the distance from orthonormality
00282   template<class ScalarType, class MV, class OP>
00283   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00284   BasicOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RefCountPtr<const MV> MX) const {
00285     const ScalarType ONE = SCT::one();
00286     int rank = MVT::GetNumberVecs(X);
00287     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00288     innerProd(X,X,MX,xTx);
00289     for (int i=0; i<rank; i++) {
00290       xTx(i,i) -= ONE;
00291     }
00292     return xTx.normFrobenius();
00293   }
00294 
00296   // Compute the distance from orthogonality
00297   template<class ScalarType, class MV, class OP>
00298   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00299   BasicOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RefCountPtr<const MV> MX1, const MV &X2) const {
00300     int r1 = MVT::GetNumberVecs(X1);
00301     int r2  = MVT::GetNumberVecs(X2);
00302     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
00303     innerProd(X2,X1,MX1,xTx);
00304     return xTx.normFrobenius();
00305   }
00306 
00308   // Find an Op-orthonormal basis for span(X) - span(W)
00309   template<class ScalarType, class MV, class OP>
00310   int BasicOrthoManager<ScalarType, MV, OP>::projectAndNormalize(
00311                                     MV &X, Teuchos::RefCountPtr<MV> MX, 
00312                                     Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00313                                     Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00314                                     Teuchos::Array<Teuchos::RefCountPtr<const MV> > Q ) const {
00315 
00316     int nq = Q.length();
00317     int xc = MVT::GetNumberVecs( X );
00318     int xr = MVT::GetVecLength( X );
00319     int rank;
00320 
00321     /* if the user doesn't want to store the coefficienets, 
00322      * allocate some local memory for them 
00323      */
00324     if ( B == Teuchos::null ) {
00325       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00326     }
00327 
00328     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00329     if (this->_hasOp) {
00330       if (MX == Teuchos::null) {
00331         // we need to allocate space for MX
00332         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00333         OPT::Apply(*(this->_Op),X,*MX);
00334       }
00335     }
00336     else {
00337       // Op == I  -->  MX = X (ignore it if the user passed it in)
00338       MX = Teuchos::rcp( &X, false );
00339     }
00340 
00341     int mxc = MVT::GetNumberVecs( *MX );
00342     int mxr = MVT::GetVecLength( *MX );
00343 
00344     // short-circuit
00345     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Anasazi::BasicOrthoManager::projectAndNormalize(): X must be non-empty" );
00346 
00347     int numbas = 0;
00348     for (int i=0; i<nq; i++) {
00349       numbas += MVT::GetNumberVecs( *Q[i] );
00350     }
00351 
00352     // check size of B
00353     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00354                         "Anasazi::BasicOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
00355     // check size of X and MX
00356     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00357                         "Anasazi::BasicOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
00358     // check size of X w.r.t. MX 
00359     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00360                         "Anasazi::BasicOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
00361     // check feasibility
00362     TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00363                         "Anasazi::BasicOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
00364 
00365     // orthogonalize all of X against Q
00366     project(X,MX,C,Q);
00367 
00368 
00369     Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(xc,1);
00370 
00371     // start working
00372     rank = 0;
00373     int numTries = 10;   // each vector in X gets 10 random chances to escape degeneracy
00374     int oldrank = -1;
00375     do {
00376       int curxsize = xc - rank;
00377 
00378       // orthonormalize X, but quit if it is rank deficient
00379       // we can't let findBasis generated random vectors to complete the basis,
00380       // because it doesn't know about Q; we will do this ourselves below
00381       rank = findBasis(X,MX,B,false,curxsize);
00382 
00383       if (rank < xc && numTries == 10) {
00384         // we quit on this vector, and for the first time;
00385         // save the coefficient information, because findBasis will overwrite it
00386         for (int i=0; i<xc; i++) {
00387           oldCoeff(i,0) = (*B)(i,rank);
00388         }
00389       }
00390 
00391       if (oldrank != -1 && rank != oldrank) {
00392         // we moved on; restore the previous coefficients
00393         for (int i=0; i<xc; i++) {
00394           (*B)(i,oldrank) = oldCoeff(i,0);
00395         }
00396       }
00397 
00398       if (rank == xc) {
00399         // we are done
00400         break;
00401       }
00402       else {
00403         TEST_FOR_EXCEPTION( rank < oldrank, OrthoError,   
00404                             "Anasazi::BasicOrthoManager::projectAndNormalize(): basis lost rank; this shouldn't happen");
00405 
00406         if (rank != oldrank) {
00407           // we added a basis vector from random info; reset the chance counter
00408           numTries = 10;
00409         }
00410 
00411         // store old rank
00412         oldrank = rank;
00413 
00414         // has this vector run out of chances to escape degeneracy?
00415         if (numTries <= 0) {
00416           break;
00417         }
00418         // use one of this vector's chances
00419         numTries--;
00420 
00421         // randomize troubled direction
00422 #ifdef ORTHO_DEBUG
00423         cout << "Random for column " << rank << endl;
00424 #endif
00425         Teuchos::RefCountPtr<MV> curX, curMX;
00426         std::vector<int> ind(1);
00427         ind[0] = rank;
00428         curX  = MVT::CloneView(X,ind);
00429         MVT::MvRandom(*curX);
00430         if (this->_hasOp) {
00431           curMX = MVT::CloneView(*MX,ind);
00432           OPT::Apply( *(this->_Op), *curX, *curMX );
00433         }
00434 
00435         // orthogonalize against Q
00436         // if !this->_hasOp, the curMX will be ignored.
00437         // we don't care about these coefficients; in fact, we need to preserve the previous coeffs
00438         project(*curX,curMX,Teuchos::null,Q);
00439       }
00440     } while (1);
00441 
00442     // this should not raise an exception; but our post-conditions oblige us to check
00443     TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00444                         "Anasazi::BasicOrthoManager::projectAndNormalize(): Debug error in rank variable." );
00445     return rank;
00446   }
00447 
00448 
00449 
00451   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00452   template<class ScalarType, class MV, class OP>
00453   int BasicOrthoManager<ScalarType, MV, OP>::normalize(
00454                                 MV &X, Teuchos::RefCountPtr<MV> MX, 
00455                                 Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00456     // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
00457     return findBasis(X, MX, B, true );
00458   }
00459 
00460 
00461 
00463   template<class ScalarType, class MV, class OP>
00464   void BasicOrthoManager<ScalarType, MV, OP>::project(
00465                           MV &X, Teuchos::RefCountPtr<MV> MX, 
00466                           Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00467                           Teuchos::Array<Teuchos::RefCountPtr<const MV> > Q) const {
00468     // For the inner product defined by the operator Op or the identity (Op == 0)
00469     //   -> Orthogonalize X against each Q[i]
00470     // Modify MX accordingly
00471     //
00472     // Note that when Op is 0, MX is not referenced
00473     //
00474     // Parameter variables
00475     //
00476     // X  : Vectors to be transformed
00477     //
00478     // MX : Image of the block vector X by the mass matrix
00479     //
00480     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independenty.
00481     //
00482 
00483     ScalarType    ONE  = SCT::one();
00484 
00485     int xc = MVT::GetNumberVecs( X );
00486     int xr = MVT::GetVecLength( X );
00487     int nq = Q.length();
00488     std::vector<int> qcs(nq);
00489     // short-circuit
00490     if (nq == 0 || xc == 0 || xr == 0) {
00491       return;
00492     }
00493     int qr = MVT::GetVecLength ( *Q[0] );
00494     // if we don't have enough C, expand it with null references
00495     // if we have too many, resize to throw away the latter ones
00496     // if we have exactly as many as we have Q, this call has no effect
00497     C.resize(nq);
00498 
00499 
00500     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00501     if (this->_hasOp) {
00502       if (MX == Teuchos::null) {
00503         // we need to allocate space for MX
00504         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00505         OPT::Apply(*(this->_Op),X,*MX);
00506       }
00507     }
00508     else {
00509       // Op == I  -->  MX = X (ignore it if the user passed it in)
00510       MX = Teuchos::rcp( &X, false );
00511     }
00512     int mxc = MVT::GetNumberVecs( *MX );
00513     int mxr = MVT::GetVecLength( *MX );
00514 
00515     // check size of X and Q w.r.t. common sense
00516     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00517                         "Anasazi::BasicOrthoManager::project(): MVT returned negative dimensions for X,MX" );
00518     // check size of X w.r.t. MX and Q
00519     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
00520                         "Anasazi::BasicOrthoManager::project(): Size of X not consistant with MX,Q" );
00521 
00522     // tally up size of all Q and check/allocate C
00523     int baslen = 0;
00524     for (int i=0; i<nq; i++) {
00525       TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 
00526                           "Anasazi::BasicOrthoManager::project(): Q lengths not mutually consistant" );
00527       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00528       TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
00529                           "Anasazi::BasicOrthoManager::project(): Q has less rows than columns" );
00530       baslen += qcs[i];
00531 
00532       // check size of C[i]
00533       if ( C[i] == Teuchos::null ) {
00534         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00535       }
00536       else {
00537         TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
00538                            "Anasazi::BasicOrthoManager::project(): Size of Q not consistant with size of C" );
00539       }
00540     }
00541 
00542     // Perform the Gram-Schmidt transformation for a block of vectors
00543 
00544     // Compute the initial Op-norms
00545     std::vector<ScalarType> oldDot( xc );
00546     MVT::MvDot( X, *MX, &oldDot );
00547 
00548     Teuchos::Array<Teuchos::RefCountPtr<MV> > MQ(nq);
00549     // Define the product Q^T * (Op*X)
00550     for (int i=0; i<nq; i++) {
00551       // Multiply Q' with MX
00552       innerProd(*Q[i],X,MX,*C[i]);
00553       // Multiply by Q and subtract the result in X
00554       MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
00555 
00556       // Update MX, with the least number of applications of Op as possible
00557       if (this->_hasOp) {
00558         if (xc <= qcs[i]) {
00559           OPT::Apply( *(this->_Op), X, *MX);
00560         }
00561         else {
00562           // this will possibly be used again below; don't delete it
00563           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
00564           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
00565           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
00566         }
00567       }
00568     }
00569 
00570     // Compute new Op-norms
00571     std::vector<ScalarType> newDot(xc);
00572     MVT::MvDot( X, *MX, &newDot );
00573 
00574     // determine (individually) whether to do another step of classical Gram-Schmidt
00575     for (int j = 0; j < xc; ++j) {
00576       
00577       if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
00578         Teuchos::TimeMonitor lcltimer( *timerReortho_ );
00579         for (int i=0; i<nq; i++) {
00580           Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
00581           
00582           // Apply another step of classical Gram-Schmidt
00583           innerProd(*Q[i],X,MX,C2);
00584           *C[i] += C2;
00585           MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
00586           
00587           // Update MX, with the least number of applications of Op as possible
00588           if (this->_hasOp) {
00589             if (MQ[i].get()) {
00590               // MQ was allocated and computed above; use it
00591               MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
00592             }
00593             else if (xc <= qcs[i]) {
00594               // MQ was not allocated and computed above; it was cheaper to use X before and it still is
00595               OPT::Apply( *(this->_Op), X, *MX);
00596             }
00597           }
00598         }
00599         break;
00600       } // if (kappa_*newDot[j] < oldDot[j])
00601     } // for (int j = 0; j < xc; ++j)
00602   }
00603 
00604 
00606   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
00607   // the rank is numvectors(X)
00608   template<class ScalarType, class MV, class OP>
00609   int BasicOrthoManager<ScalarType, MV, OP>::findBasis(
00610                 MV &X, Teuchos::RefCountPtr<MV> MX, 
00611                 Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00612                 bool completeBasis, int howMany ) const {
00613     // For the inner product defined by the operator Op or the identity (Op == 0)
00614     //   -> Orthonormalize X 
00615     // Modify MX accordingly
00616     //
00617     // Note that when Op is 0, MX is not referenced
00618     //
00619     // Parameter variables
00620     //
00621     // X  : Vectors to be orthonormalized
00622     //
00623     // MX : Image of the multivector X under the operator Op
00624     //
00625     // Op  : Pointer to the operator for the inner product
00626     //
00627     // TODO: add reference
00628     // kappa= Coefficient determining when to perform a second Gram-Schmidt step
00629     //        Default value = 1.5625 = (1.25)^2 (as suggested in Parlett's book)
00630     //
00631 
00632     const ScalarType ONE  = SCT::one();
00633     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
00634     const ScalarType EPS  = SCT::eps();
00635 
00636     int xc = MVT::GetNumberVecs( X );
00637     int xr = MVT::GetVecLength( X );
00638 
00639     if (howMany == -1) {
00640       howMany = xc;
00641     }
00642 
00643     /*******************************************************
00644      *  If _hasOp == false, we will not reference MX below *
00645      *******************************************************/
00646 
00647     // if Op==null, MX == X (via pointer)
00648     // Otherwise, either the user passed in MX or we will allocated and compute it
00649     if (this->_hasOp) {
00650       if (MX == Teuchos::null) {
00651         // we need to allocate space for MX
00652         MX = MVT::Clone(X,xc);
00653         OPT::Apply(*(this->_Op),X,*MX);
00654       }
00655     }
00656 
00657     /* if the user doesn't want to store the coefficienets, 
00658      * allocate some local memory for them 
00659      */
00660     if ( B == Teuchos::null ) {
00661       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00662     }
00663 
00664     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00665     int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX )  : xr;
00666 
00667     // check size of C, B
00668     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
00669                         "Anasazi::BasicOrthoManager::findBasis(): X must be non-empty" );
00670     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00671                         "Anasazi::BasicOrthoManager::findBasis(): Size of X not consistant with size of B" );
00672     TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
00673                         "Anasazi::BasicOrthoManager::findBasis(): Size of X not consistant with size of MX" );
00674     TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 
00675                         "Anasazi::BasicOrthoManager::findBasis(): Size of X not feasible for normalization" );
00676     TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
00677                         "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" );
00678 
00679     /* xstart is which column we are starting the process with, based on howMany
00680      * columns before xstart are assumed to be Op-orthonormal already
00681      */
00682     int xstart = xc - howMany;
00683 
00684     for (int j = xstart; j < xc; j++) {
00685 
00686       // numX is 
00687       // * number of currently orthonormal columns of X
00688       // * the index of the current column of X
00689       int numX = j;
00690 
00691       // Get a view of the vector currently being worked on.
00692       std::vector<int> index(1);
00693       index[0] = numX;
00694       Teuchos::RefCountPtr<MV> Xj = MVT::CloneView( X, index );
00695       Teuchos::RefCountPtr<MV> MXj;
00696       if ((this->_hasOp)) {
00697         // MXj is a view of the current vector in MX
00698         MXj = MVT::CloneView( *MX, index );
00699       }
00700       else {
00701         // MXj is a pointer to Xj, and MUST NOT be modified
00702         MXj = Xj;
00703       }
00704 
00705       // Get a view of the previous vectors.
00706       std::vector<int> prev_idx( numX );
00707       Teuchos::RefCountPtr<const MV> prevX, prevMX;
00708 
00709       if (numX > 0) {
00710         for (int i=0; i<numX; i++) {
00711           prev_idx[i] = i;
00712         }
00713         prevX = MVT::CloneView( X, prev_idx );
00714         if (this->_hasOp) {
00715           prevMX = MVT::CloneView( *MX, prev_idx );
00716         }
00717       } 
00718 
00719       bool rankDef = true;
00720       /* numTrials>0 will denote that the current vector was randomized for the purpose
00721        * of finding a basis vector, and that the coefficients of that vector should
00722        * not be stored in B
00723        */
00724       for (int numTrials = 0; numTrials < 10; numTrials++) {
00725 
00726         // Make storage for these Gram-Schmidt iterations.
00727         Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
00728         std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
00729 
00730         //
00731         // Save old MXj vector and compute Op-norm
00732         //
00733         Teuchos::RefCountPtr<MV> oldMXj = MVT::CloneCopy( *MXj ); 
00734         MVT::MvDot( *Xj, *MXj, &oldDot );
00735         // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
00736         TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
00737                             "Anasazi::BasicOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
00738 
00739         if (numX > 0) {
00740           // Apply the first step of Gram-Schmidt
00741 
00742           // product <- prevX^T MXj
00743           innerProd(*prevX,*Xj,MXj,product);
00744 
00745           // Xj <- Xj - prevX prevX^T MXj   
00746           //     = Xj - prevX product
00747           MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
00748 
00749           // Update MXj
00750           if (this->_hasOp) {
00751             // MXj <- Op*Xj_new
00752             //      = Op*(Xj_old - prevX prevX^T MXj)
00753             //      = MXj - prevMX product
00754             MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
00755           }
00756 
00757           // Compute new Op-norm
00758           MVT::MvDot( *Xj, *MXj, &newDot );
00759 
00760           // Check if a correction is needed.
00761           if ( SCT::magnitude(kappa_*newDot[0]) < SCT::magnitude(oldDot[0]) ) {
00762             // Apply the second step of Gram-Schmidt
00763             // This is the same as above
00764             Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
00765 
00766             innerProd(*prevX,*Xj,MXj,P2);
00767             product += P2;
00768             MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
00769             if ((this->_hasOp)) {
00770               MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
00771             }
00772           } // if (kappa_*newDot[0] < oldDot[0])
00773 
00774         } // if (numX > 0)
00775 
00776         // Compute Op-norm with old MXj
00777         MVT::MvDot( *Xj, *oldMXj, &newDot );
00778 
00779         // save the coefficients, if we are working on the original vector and not a randomly generated one
00780         if (numTrials == 0) {
00781           for (int i=0; i<numX; i++) {
00782             (*B)(i,j) = product(i,0);
00783           }
00784         }
00785 
00786         // Check if Xj has any directional information left after the orthogonalization.
00787 #ifdef ORTHO_DEBUG
00788         cout << "olddot: " << SCT::magnitude(oldDot[0]) << "    newdot: " << SCT::magnitude(newDot[0]);
00789 #endif
00790         if ( SCT::magnitude(newDot[0]) > SCT::magnitude(oldDot[0]*EPS*EPS) && SCT::real(newDot[0]) > ZERO ) {
00791 #ifdef ORTHO_DEBUG
00792           cout << " ACCEPTED" << endl;
00793 #endif
00794           // Normalize Xj.
00795           // Xj <- Xj / sqrt(newDot)
00796           ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
00797 
00798           MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
00799           if (this->_hasOp) {
00800             // Update MXj.
00801             MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
00802           }
00803 
00804           // save it, if it corresponds to the original vector and not a randomly generated one
00805           if (numTrials == 0) {
00806             (*B)(j,j) = diag;
00807           }
00808 
00809           // We are not rank deficient in this vector. Move on to the next vector in X.
00810           rankDef = false;
00811           break;
00812         }
00813         else {
00814 #ifdef ORTHO_DEBUG
00815           cout << " REJECTED" << endl;
00816 #endif
00817           // There was nothing left in Xj after orthogonalizing against previous columns in X.
00818           // X is rank deficient.
00819           // reflect this in the coefficients
00820           (*B)(j,j) = ZERO;
00821 
00822           if (completeBasis) {
00823             // Fill it with random information and keep going.
00824 #ifdef ORTHO_DEBUG
00825             cout << "Random for column " << j << endl;
00826 #endif
00827             MVT::MvRandom( *Xj );
00828             if (this->_hasOp) {
00829               OPT::Apply( *(this->_Op), *Xj, *MXj );
00830             }
00831           }
00832           else {
00833             rankDef = true;
00834             break;
00835           }
00836 
00837         } // if (norm > oldDot*EPS*EPS)
00838 
00839       }  // for (numTrials = 0; numTrials < 10; ++numTrials)
00840 
00841       // if rankDef == true, then quit and notify user of rank obtained
00842       if (rankDef == true) {
00843         (*B)(j,j) = ZERO;
00844         MVT::MvInit( *Xj, ZERO );
00845         if (this->_hasOp) {
00846           MVT::MvInit( *MXj, ZERO );
00847         }
00848         TEST_FOR_EXCEPTION( completeBasis, OrthoError, 
00849                             "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" );
00850         return j;
00851       }
00852 
00853     } // for (j = 0; j < xc; ++j)
00854 
00855     return xc;
00856   }
00857 
00858 } // namespace Anasazi
00859 
00860 #endif // ANASAZI_BASIC_ORTHOMANAGER_HPP
00861 

Generated on Thu Sep 18 12:31:33 2008 for Anasazi by doxygen 1.3.9.1