AnasaziBasicOrthoManager.hpp

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