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

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