Anasazi Version of the Day
AnasaziSVQBOrthoManager.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 
00035 #ifndef ANASAZI_SVQB_ORTHOMANAGER_HPP
00036 #define ANASAZI_SVQB_ORTHOMANAGER_HPP
00037 
00047 #include "AnasaziConfigDefs.hpp"
00048 #include "AnasaziMultiVecTraits.hpp"
00049 #include "AnasaziOperatorTraits.hpp"
00050 #include "AnasaziMatOrthoManager.hpp"
00051 #include "Teuchos_LAPACK.hpp"
00052 
00053 namespace Anasazi {
00054 
00055   template<class ScalarType, class MV, class OP>
00056   class SVQBOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
00057 
00058   private:
00059     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00060     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
00061     typedef Teuchos::ScalarTraits<MagnitudeType>  SCTM;
00062     typedef MultiVecTraits<ScalarType,MV>      MVT;
00063     typedef OperatorTraits<ScalarType,MV,OP>   OPT;
00064     std::string dbgstr;
00065 
00066 
00067   public:
00068 
00070 
00071 
00072     SVQBOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null, bool debug = false );
00073 
00074 
00076     ~SVQBOrthoManager() {};
00078 
00079 
00080 
00082 
00083 
00084 
00123     void projectMat ( 
00124           MV &X, 
00125           Teuchos::Array<Teuchos::RCP<const MV> >  Q,
00126           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00127               = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00128           Teuchos::RCP<MV> MX                        = Teuchos::null,
00129           Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00130         ) const;
00131 
00132 
00171     int normalizeMat ( 
00172           MV &X, 
00173           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00174           Teuchos::RCP<MV> MX                                         = Teuchos::null
00175         ) const;
00176 
00177 
00237     int projectAndNormalizeMat ( 
00238           MV &X, 
00239           Teuchos::Array<Teuchos::RCP<const MV> >  Q,
00240           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00241               = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00242           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null, 
00243           Teuchos::RCP<MV>                                         MX = Teuchos::null,
00244           Teuchos::Array<Teuchos::RCP<const MV> >                  MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00245         ) const;
00246 
00248 
00250 
00251 
00256     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00257     orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
00258 
00263     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00264     orthogErrorMat(
00265           const MV &X, 
00266           const MV &Y,
00267           Teuchos::RCP<const MV> MX = Teuchos::null, 
00268           Teuchos::RCP<const MV> MY = Teuchos::null
00269         ) const;
00270 
00272 
00273   private:
00274     
00275     MagnitudeType eps_;
00276     bool debug_;
00277   
00278     // ! Routine to find an orthogonal/orthonormal basis
00279     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00280                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00281                    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00282                    Teuchos::Array<Teuchos::RCP<const MV> > Q,
00283                    bool normalize_in ) const;
00284   };
00285 
00286 
00288   // Constructor
00289   template<class ScalarType, class MV, class OP>
00290   SVQBOrthoManager<ScalarType,MV,OP>::SVQBOrthoManager( Teuchos::RCP<const OP> Op, bool debug) 
00291     : MatOrthoManager<ScalarType,MV,OP>(Op), dbgstr("                    *** "), debug_(debug) {
00292     
00293     Teuchos::LAPACK<int,MagnitudeType> lapack;
00294     eps_ = lapack.LAMCH('E');
00295     if (debug_) {
00296       std::cout << "eps_ == " << eps_ << std::endl;
00297     }
00298   }
00299 
00300 
00302   // Compute the distance from orthonormality
00303   template<class ScalarType, class MV, class OP>
00304   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00305   SVQBOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const {
00306     const ScalarType ONE = SCT::one();
00307     int rank = MVT::GetNumberVecs(X);
00308     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00309     innerProdMat(X,X,xTx,MX,MX);
00310     for (int i=0; i<rank; i++) {
00311       xTx(i,i) -= ONE;
00312     }
00313     return xTx.normFrobenius();
00314   }
00315 
00317   // Compute the distance from orthogonality
00318   template<class ScalarType, class MV, class OP>
00319   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00320   SVQBOrthoManager<ScalarType,MV,OP>::orthogErrorMat(
00321           const MV &X, 
00322           const MV &Y,
00323           Teuchos::RCP<const MV> MX,
00324           Teuchos::RCP<const MV> MY
00325       ) const {
00326     int r1 = MVT::GetNumberVecs(X);
00327     int r2 = MVT::GetNumberVecs(Y);
00328     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
00329     innerProdMat(X,Y,xTx,MX,MY);
00330     return xTx.normFrobenius();
00331   }
00332 
00333 
00334 
00336   template<class ScalarType, class MV, class OP>
00337   void SVQBOrthoManager<ScalarType, MV, OP>::projectMat(
00338           MV &X, 
00339           Teuchos::Array<Teuchos::RCP<const MV> >  Q,
00340           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00341           Teuchos::RCP<MV> MX,
00342           Teuchos::Array<Teuchos::RCP<const MV> > MQ) const {
00343     (void)MQ;
00344     findBasis(X,MX,C,Teuchos::null,Q,false);
00345   }
00346 
00347 
00348 
00350   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00351   template<class ScalarType, class MV, class OP>
00352   int SVQBOrthoManager<ScalarType, MV, OP>::normalizeMat(
00353           MV &X, 
00354           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00355           Teuchos::RCP<MV> MX) const {
00356     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C;
00357     Teuchos::Array<Teuchos::RCP<const MV> > Q;
00358     return findBasis(X,MX,C,B,Q,true);
00359   }
00360 
00361 
00362 
00364   // Find an Op-orthonormal basis for span(X) - span(W)
00365   template<class ScalarType, class MV, class OP>
00366   int SVQBOrthoManager<ScalarType, MV, OP>::projectAndNormalizeMat(
00367           MV &X, 
00368           Teuchos::Array<Teuchos::RCP<const MV> >  Q,
00369           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00370           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00371           Teuchos::RCP<MV> MX,
00372           Teuchos::Array<Teuchos::RCP<const MV> > MQ) const {
00373     (void)MQ;
00374     return findBasis(X,MX,C,B,Q,true);
00375   }
00376 
00377 
00378 
00379 
00381   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
00382   // the rank is numvectors(X)
00383   // 
00384   // Tracking the coefficients (C[i] and B) for this code is complicated by the fact that the loop
00385   // structure looks like
00386   // do
00387   //    project
00388   //    do
00389   //       ortho
00390   //    end
00391   // end
00392   // However, the recurrence for the coefficients is not complicated:
00393   // B = I
00394   // C = 0
00395   // do 
00396   //    project yields newC
00397   //    C = C + newC*B
00398   //    do 
00399   //       ortho yields newR
00400   //       B = newR*B
00401   //    end
00402   // end
00403   // This holds for each individual C[i] (which correspond to the list of bases we are orthogonalizing
00404   // against).
00405   //
00406   template<class ScalarType, class MV, class OP>
00407   int SVQBOrthoManager<ScalarType, MV, OP>::findBasis(
00408                 MV &X, Teuchos::RCP<MV> MX, 
00409                 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00410                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00411                 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00412                 bool normalize_in) const {
00413 
00414     const ScalarType ONE  = SCT::one();
00415     const MagnitudeType MONE = SCTM::one();
00416     const MagnitudeType ZERO = SCTM::zero();
00417 
00418     int numGS = 0,
00419         numSVQB = 0,
00420         numRand = 0;
00421 
00422     // get sizes of X,MX
00423     int xc = MVT::GetNumberVecs(X);
00424     int xr = MVT::GetVecLength( X );
00425 
00426     // get sizes of Q[i]
00427     int nq = Q.length();
00428     int qr = (nq == 0) ? 0 : MVT::GetVecLength(*Q[0]);
00429     int qsize = 0;
00430     std::vector<int> qcs(nq);
00431     for (int i=0; i<nq; i++) {
00432       qcs[i] = MVT::GetNumberVecs(*Q[i]);
00433       qsize += qcs[i];
00434     }
00435 
00436     if (normalize_in == true && qsize + xc > xr) {
00437       // not well-posed
00438       TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00439                           "Anasazi::SVQBOrthoManager::findBasis(): Orthogonalization constraints not feasible" );
00440     }
00441 
00442     // try to short-circuit as early as possible
00443     if (normalize_in == false && (qsize == 0 || xc == 0)) {
00444       // nothing to do
00445       return 0;
00446     }
00447     else if (normalize_in == true && (xc == 0 || xr == 0)) {
00448       // normalize requires X not empty
00449       TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00450                           "Anasazi::SVQBOrthoManager::findBasis(): X must be non-empty" );
00451     }
00452 
00453     // check that Q matches X
00454     TEST_FOR_EXCEPTION( qsize != 0 && qr != xr , std::invalid_argument, 
00455                         "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with size of Q" );
00456 
00457     /* If we don't have enough C, expanding it creates null references
00458      * If we have too many, resizing just throws away the later ones
00459      * If we have exactly as many as we have Q, this call has no effect
00460      */
00461     C.resize(nq);
00462     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > newC(nq);
00463     // check the size of the C[i] against the Q[i] and consistency between Q[i]
00464     for (int i=0; i<nq; i++) {
00465       // check size of Q[i]
00466       TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 
00467                           "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not mutually consistant" );
00468       TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
00469                           "Anasazi::SVQBOrthoManager::findBasis(): Q has less rows than columns" );
00470       // check size of C[i]
00471       if ( C[i] == Teuchos::null ) {
00472         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00473       }
00474       else {
00475         TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc, std::invalid_argument, 
00476                             "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not consistant with C" );
00477       }
00478       // clear C[i]
00479       C[i]->putScalar(ZERO);
00480       newC[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(*C[i]) );
00481     }
00482 
00483 
00485     // Allocate necessary storage
00486     // C were allocated above
00487     // Allocate MX and B (if necessary)
00488     // Set B = I
00489     if (normalize_in == true) {
00490       if ( B == Teuchos::null ) {
00491         B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00492       }
00493       TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00494                           "Anasazi::SVQBOrthoManager::findBasis(): Size of B not consistant with X" );
00495       // set B to I
00496       B->putScalar(ZERO);
00497       for (int i=0; i<xc; i++) {
00498         (*B)(i,i) = MONE;
00499       }
00500     }
00501     /******************************************
00502      *  If _hasOp == false, DO NOT MODIFY MX  *
00503      ****************************************** 
00504      * if Op==null, MX == X (via pointer)
00505      * Otherwise, either the user passed in MX or we will allocate and compute it
00506      *
00507      * workX will be a multivector of the same size as X, used to perform X*S when normalizing
00508      */
00509     Teuchos::RCP<MV> workX;
00510     if (normalize_in) {
00511       workX = MVT::Clone(X,xc);
00512     }
00513     if (this->_hasOp) {
00514       if (MX == Teuchos::null) {
00515         // we need to allocate space for MX
00516         MX = MVT::Clone(X,xc);
00517         OPT::Apply(*(this->_Op),X,*MX);
00518         this->_OpCounter += MVT::GetNumberVecs(X);
00519       }
00520     }
00521     else {
00522       MX = Teuchos::rcpFromRef(X);
00523     }
00524     std::vector<ScalarType> normX(xc), invnormX(xc);
00525     Teuchos::SerialDenseMatrix<int,ScalarType> XtMX(xc,xc), workU(1,1);
00526     Teuchos::LAPACK<int,ScalarType> lapack;
00527     /**********************************************************************
00528      * allocate storage for eigenvectors,eigenvalues of X^T Op X, and for
00529      * the work space needed to compute this xc-by-xc eigendecomposition
00530      **********************************************************************/
00531     std::vector<ScalarType> work;
00532     std::vector<MagnitudeType> lambda, lambdahi, rwork;
00533     if (normalize_in) {
00534       // get size of work from ILAENV
00535       int lwork = lapack.ILAENV(1,"hetrd","VU",xc,-1,-1,-1);
00536       // lwork >= (nb+1)*n for complex
00537       // lwork >= (nb+2)*n for real
00538       TEST_FOR_EXCEPTION( lwork < 0, OrthoError, 
00539                           "Anasazi::SVQBOrthoManager::findBasis(): Error code from ILAENV" );
00540 
00541       lwork = (lwork+2)*xc;
00542       work.resize(lwork);
00543       // size of rwork is max(1,3*xc-2)
00544       lwork = (3*xc-2 > 1) ? 3*xc - 2 : 1;
00545       rwork.resize(lwork);
00546       // size of lambda is xc
00547       lambda.resize(xc);
00548       lambdahi.resize(xc);
00549       workU.reshape(xc,xc);
00550     }
00551 
00552     // test sizes of X,MX
00553     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00554     int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX )  : xr;
00555     TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
00556                         "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with MX" );
00557 
00558     // sentinel to continue the outer loop (perform another projection step)
00559     bool doGramSchmidt = true;          
00560     // variable for testing orthonorm/orthog 
00561     MagnitudeType tolerance = MONE/SCTM::squareroot(eps_);
00562 
00563     // outer loop
00564     while (doGramSchmidt) {
00565 
00567       // perform projection
00568       if (qsize > 0) {
00569 
00570         numGS++;
00571 
00572         // Compute the norms of the vectors
00573         {
00574           std::vector<MagnitudeType> normX_mag(xc);
00575           normMat(X,normX_mag,MX);
00576           for (int i=0; i<xc; ++i) {
00577             normX[i] = normX_mag[i];
00578             invnormX[i] = (normX_mag[i] == ZERO) ? ZERO : MONE/normX_mag[i]; 
00579           }
00580         }
00581         // normalize the vectors
00582         MVT::MvScale(X,invnormX);
00583         if (this->_hasOp) {
00584           MVT::MvScale(*MX,invnormX);
00585         }
00586         // check that vectors are normalized now
00587         if (debug_) {
00588           std::vector<MagnitudeType> nrm2(xc);
00589           std::cout << dbgstr << "max post-scale norm: (with/without MX) : ";
00590           MagnitudeType maxpsnw = ZERO, maxpsnwo = ZERO;
00591           normMat(X,nrm2,MX);
00592           for (int i=0; i<xc; i++) {
00593             maxpsnw = (nrm2[i] > maxpsnw ? nrm2[i] : maxpsnw);
00594           }
00595           this->norm(X,nrm2);
00596           for (int i=0; i<xc; i++) {
00597             maxpsnwo = (nrm2[i] > maxpsnwo ? nrm2[i] : maxpsnwo);
00598           }
00599           std::cout << "(" << maxpsnw << "," << maxpsnwo << ")" << std::endl;
00600         }
00601         // project the vectors onto the Qi
00602         for (int i=0; i<nq; i++) {
00603           innerProdMat(*Q[i],X,*newC[i],Teuchos::null,MX);
00604         }
00605         // remove the components in Qi from X
00606         for (int i=0; i<nq; i++) {
00607           MVT::MvTimesMatAddMv(-ONE,*Q[i],*newC[i],ONE,X);
00608         }
00609         // un-scale the vectors 
00610         MVT::MvScale(X,normX);
00611         // Recompute the vectors in MX
00612         if (this->_hasOp) {
00613           OPT::Apply(*(this->_Op),X,*MX);
00614           this->_OpCounter += MVT::GetNumberVecs(X);
00615         }
00616 
00617         //          
00618         // Compute largest column norm of 
00619         //            (  C[0]   )  
00620         //        C = (  ....   )
00621         //            ( C[nq-1] )
00622         MagnitudeType maxNorm = ZERO;
00623         for  (int j=0; j<xc; j++) {
00624           MagnitudeType sum = ZERO;
00625           for (int k=0; k<nq; k++) {
00626             for (int i=0; i<qcs[k]; i++) {
00627               sum += SCT::magnitude((*newC[k])(i,j))*SCT::magnitude((*newC[k])(i,j));
00628             }
00629           }
00630           maxNorm = (sum > maxNorm) ? sum : maxNorm;
00631         }
00632               
00633         // do we perform another GS?
00634         if (maxNorm < 0.36) {
00635           doGramSchmidt = false;
00636         }
00637 
00638         // unscale newC to reflect the scaling of X
00639         for (int k=0; k<nq; k++) {
00640           for (int j=0; j<xc; j++) {
00641             for (int i=0; i<qcs[k]; i++) {
00642               (*newC[k])(i,j) *= normX[j];
00643             }
00644           }
00645         }
00646         // accumulate into C
00647         if (normalize_in) {
00648           // we are normalizing
00649           int info;
00650           for (int i=0; i<nq; i++) {
00651             info = C[i]->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,*newC[i],*B,ONE);
00652             TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
00653           }
00654         }
00655         else {
00656           // not normalizing
00657           for (int i=0; i<nq; i++) {
00658             (*C[i]) += *newC[i];
00659           }
00660         }
00661       }
00662       else { // qsize == 0... don't perform projection
00663         // don't do any more outer loops; all we need is to call the normalize code below
00664         doGramSchmidt = false;
00665       }
00666 
00667 
00669       // perform normalization
00670       if (normalize_in) {
00671 
00672         MagnitudeType condT = tolerance;
00673         
00674         while (condT >= tolerance) {
00675 
00676           numSVQB++;
00677 
00678           // compute X^T Op X
00679           innerProdMat(X,X,XtMX,MX,MX);
00680 
00681           // compute scaling matrix for XtMX: D^{.5} and D^{-.5} (D-half  and  D-half-inv)
00682           std::vector<MagnitudeType> Dh(xc), Dhi(xc);
00683           for (int i=0; i<xc; i++) {
00684             Dh[i]  = SCT::magnitude(SCT::squareroot(XtMX(i,i)));
00685             Dhi[i] = (Dh[i] == ZERO ? ZERO : MONE/Dh[i]);
00686           }
00687           // scale XtMX :   S = D^{-.5} * XtMX * D^{-.5}
00688           for (int i=0; i<xc; i++) {
00689             for (int j=0; j<xc; j++) {
00690               XtMX(i,j) *= Dhi[i]*Dhi[j];
00691             }
00692           }
00693 
00694           // compute the eigenvalue decomposition of S=U*Lambda*U^T (using upper part)
00695           int info;
00696           lapack.HEEV('V', 'U', xc, XtMX.values(), XtMX.stride(), &lambda[0], &work[0], work.size(), &rwork[0], &info);
00697           std::stringstream os;
00698           os << "Anasazi::SVQBOrthoManager::findBasis(): Error code " << info << " from HEEV";
00699           TEST_FOR_EXCEPTION( info != 0, OrthoError, 
00700                               os.str() );
00701           if (debug_) {
00702             std::cout << dbgstr << "eigenvalues of XtMX: (";
00703             for (int i=0; i<xc-1; i++) {
00704               std::cout << lambda[i] << ",";
00705             }
00706             std::cout << lambda[xc-1] << ")" << std::endl;
00707           }
00708 
00709           // remember, HEEV orders the eigenvalues from smallest to largest
00710           // examine condition number of Lambda, compute Lambda^{-.5}
00711           MagnitudeType maxLambda = lambda[xc-1],
00712                         minLambda = lambda[0];
00713           int iZeroMax = -1;
00714           for (int i=0; i<xc; i++) {
00715             if (lambda[i] <= 10*eps_*maxLambda) {      // finish: this was eps_*eps_*maxLambda
00716               iZeroMax = i;
00717               lambda[i]  = ZERO;
00718               lambdahi[i] = ZERO;
00719             }
00720             /*
00721             else if (lambda[i] < eps_*maxLambda) {
00722               lambda[i]  = SCTM::squareroot(eps_*maxLambda);
00723               lambdahi[i] = MONE/lambda[i];
00724             }
00725             */
00726             else {
00727               lambda[i] = SCTM::squareroot(lambda[i]);
00728               lambdahi[i] = MONE/lambda[i];
00729             }
00730           }
00731 
00732           // compute X * D^{-.5} * U * Lambda^{-.5} and new Op*X
00733           //
00734           // copy X into workX
00735           std::vector<int> ind(xc);
00736           for (int i=0; i<xc; i++) {ind[i] = i;}
00737           MVT::SetBlock(X,ind,*workX);
00738           //
00739           // compute D^{-.5}*U*Lambda^{-.5} into workU
00740           workU.assign(XtMX);
00741           for (int j=0; j<xc; j++) {
00742             for (int i=0; i<xc; i++) {
00743               workU(i,j) *= Dhi[i]*lambdahi[j];
00744             }
00745           }
00746           // compute workX * workU into X
00747           MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,X);
00748           //
00749           // note, it seems important to apply Op exactly for large condition numbers.
00750           // for small condition numbers, we can update MX "implicitly"
00751           // this trick reduces the number of applications of Op
00752           if (this->_hasOp) {
00753             if (maxLambda >= tolerance * minLambda) {
00754               // explicit update of MX
00755               OPT::Apply(*(this->_Op),X,*MX);
00756               this->_OpCounter += MVT::GetNumberVecs(X);
00757             }
00758             else {
00759               // implicit update of MX
00760               // copy MX into workX
00761               MVT::SetBlock(*MX,ind,*workX);
00762               //
00763               // compute workX * workU into MX
00764               MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,*MX);
00765             }
00766           }
00767 
00768           // accumulate new B into previous B
00769           // B = Lh * U^H * Dh * B
00770           for (int j=0; j<xc; j++) {
00771             for (int i=0; i<xc; i++) {
00772               workU(i,j) = Dh[i] * (*B)(i,j);
00773             }
00774           }
00775           info = B->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,XtMX,workU,ZERO);
00776           TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
00777           for (int j=0; j<xc ;j++) {
00778             for (int i=0; i<xc; i++) {
00779               (*B)(i,j) *= lambda[i];
00780             }
00781           }
00782 
00783           // check iZeroMax (rank indicator)
00784           if (iZeroMax >= 0) {
00785             if (debug_) {
00786               std::cout << dbgstr << "augmenting multivec with " << iZeroMax+1 << " random directions" << std::endl;
00787             }
00788 
00789             numRand++;
00790             // put random info in the first iZeroMax+1 vectors of X,MX
00791             std::vector<int> ind2(iZeroMax+1);
00792             for (int i=0; i<iZeroMax+1; i++) {
00793               ind2[i] = i;
00794             }
00795             Teuchos::RCP<MV> Xnull,MXnull;
00796             Xnull = MVT::CloneViewNonConst(X,ind2);
00797             MVT::MvRandom(*Xnull);
00798             if (this->_hasOp) {
00799               MXnull = MVT::CloneViewNonConst(*MX,ind2);
00800               OPT::Apply(*(this->_Op),*Xnull,*MXnull);
00801               this->_OpCounter += MVT::GetNumberVecs(*Xnull);
00802               MXnull = Teuchos::null;
00803             }
00804             Xnull = Teuchos::null;
00805             condT = tolerance;
00806             doGramSchmidt = true;
00807             break; // break from while(condT > tolerance)
00808           }
00809 
00810           condT = SCTM::magnitude(maxLambda / minLambda);
00811           if (debug_) {
00812             std::cout << dbgstr << "condT: " << condT << std::endl;
00813           }
00814           
00815         } // end while (condT >= tolerance)
00816 
00817         if ((doGramSchmidt == false) && (condT > SCTM::squareroot(tolerance))) {
00818           doGramSchmidt = true;
00819         }
00820       }
00821       // end if(normalize_in)
00822        
00823     } // end while (doGramSchmidt)
00824 
00825     if (debug_) {
00826       std::cout << dbgstr << "(numGS,numSVQB,numRand)                : " 
00827            << "(" << numGS 
00828            << "," << numSVQB 
00829            << "," << numRand 
00830            << ")" << std::endl;
00831     }
00832     return xc;
00833   }
00834 
00835 
00836 } // namespace Anasazi
00837 
00838 #endif // ANASAZI_SVQB_ORTHOMANAGER_HPP
00839 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends