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