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