Anasazi Version of the Day
AnasaziICGSOrthoManager.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_ICSG_ORTHOMANAGER_HPP
00035 #define ANASAZI_ICSG_ORTHOMANAGER_HPP
00036 
00044 #include "AnasaziConfigDefs.hpp"
00045 #include "AnasaziMultiVecTraits.hpp"
00046 #include "AnasaziOperatorTraits.hpp"
00047 #include "AnasaziGenOrthoManager.hpp"
00048 #include "Teuchos_TimeMonitor.hpp"
00049 #include "Teuchos_LAPACK.hpp"
00050 #include "Teuchos_BLAS.hpp"
00051 #ifdef ANASAZI_ICGS_DEBUG
00052 #  include <Teuchos_FancyOStream.hpp>
00053 #endif
00054 
00055 namespace Anasazi {
00056 
00057   template<class ScalarType, class MV, class OP>
00058   class ICGSOrthoManager : public GenOrthoManager<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     ICGSOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null, int numIters = 2,
00073                       typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0,
00074                       typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 );
00075 
00076 
00078     ~ICGSOrthoManager() {}
00080 
00081 
00083 
00084 
00166     void projectGen( 
00167           MV &S,
00168           Teuchos::Array<Teuchos::RCP<const MV> > X,
00169           Teuchos::Array<Teuchos::RCP<const MV> > Y,
00170           bool isBiOrtho,
00171           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00172               = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00173           Teuchos::RCP<MV> MS                        = Teuchos::null,
00174           Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)),
00175           Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00176         ) const;
00177 
00178 
00270     int projectAndNormalizeGen (
00271           MV &S,
00272           Teuchos::Array<Teuchos::RCP<const MV> > X,
00273           Teuchos::Array<Teuchos::RCP<const MV> > Y,
00274           bool isBiOrtho,
00275           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00276               = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00277           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00278           Teuchos::RCP<MV> MS                        = Teuchos::null,
00279           Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)),
00280           Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00281         ) const;
00282 
00283 
00285 
00286 
00288 
00289 
00290 
00302     void projectMat ( 
00303           MV &X, 
00304           Teuchos::Array<Teuchos::RCP<const MV> > Q,
00305           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00306               = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00307           Teuchos::RCP<MV> MX                        = Teuchos::null,
00308           Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00309         ) const;
00310 
00311 
00320     int normalizeMat ( 
00321           MV &X, 
00322           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00323           Teuchos::RCP<MV> MX                                         = Teuchos::null
00324         ) const;
00325 
00326 
00335     int projectAndNormalizeMat ( 
00336           MV &X, 
00337           Teuchos::Array<Teuchos::RCP<const MV> >  Q,
00338           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00339               = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00340           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B                  = Teuchos::null, 
00341           Teuchos::RCP<MV> MX                        = Teuchos::null,
00342           Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00343         ) const;
00344 
00346 
00348 
00349 
00354     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00355     orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
00356 
00361     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00362     orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const;
00363 
00365 
00367 
00368 
00370     void setNumIters( int numIters ) {
00371       numIters_ = numIters;
00372       TEUCHOS_TEST_FOR_EXCEPTION(numIters_ < 1,std::invalid_argument,
00373           "Anasazi::ICGSOrthoManager::setNumIters(): input must be >= 1.");
00374     }
00375 
00377     int getNumIters() const { return numIters_; } 
00378 
00380 
00381   private:
00382     MagnitudeType eps_;
00383     MagnitudeType tol_;
00384     
00386     int numIters_;
00387   
00388     // ! Routine to find an orthonormal basis
00389     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00390                          Teuchos::SerialDenseMatrix<int,ScalarType> &B, 
00391                          bool completeBasis, int howMany = -1) const;
00392   };
00393 
00394 
00395 
00397   // Constructor
00398   template<class ScalarType, class MV, class OP>
00399   ICGSOrthoManager<ScalarType,MV,OP>::ICGSOrthoManager( Teuchos::RCP<const OP> Op,
00400                                                         int numIters,
00401                                                         typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps,
00402                                                         typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol) : 
00403     GenOrthoManager<ScalarType,MV,OP>(Op), eps_(eps), tol_(tol)
00404   {
00405     setNumIters(numIters);
00406     TEUCHOS_TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument,
00407         "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"eps\" must be non-negative.");
00408     if (eps_ == 0) {
00409       Teuchos::LAPACK<int,MagnitudeType> lapack;
00410       eps_ = lapack.LAMCH('E');
00411       eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.50);
00412     }
00413     TEUCHOS_TEST_FOR_EXCEPTION(
00414         tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()),
00415         std::invalid_argument,
00416         "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"tol\" must be in [0,1].");
00417   }
00418 
00419 
00420 
00422   // Compute the distance from orthonormality
00423   template<class ScalarType, class MV, class OP>
00424   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00425   ICGSOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const {
00426     const ScalarType ONE = SCT::one();
00427     int rank = MVT::GetNumberVecs(X);
00428     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00429     MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X,X,xTx,MX,MX);
00430     for (int i=0; i<rank; i++) {
00431       xTx(i,i) -= ONE;
00432     }
00433     return xTx.normFrobenius();
00434   }
00435 
00436 
00437 
00439   // Compute the distance from orthogonality
00440   template<class ScalarType, class MV, class OP>
00441   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00442   ICGSOrthoManager<ScalarType,MV,OP>::orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const {
00443     int r1 = MVT::GetNumberVecs(X1);
00444     int r2  = MVT::GetNumberVecs(X2);
00445     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
00446     MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X1,X2,xTx,MX1,MX2);
00447     return xTx.normFrobenius();
00448   }
00449 
00450 
00451 
00453   template<class ScalarType, class MV, class OP>
00454   void ICGSOrthoManager<ScalarType, MV, OP>::projectMat(
00455           MV &X, 
00456           Teuchos::Array<Teuchos::RCP<const MV> > Q,
00457           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00458           Teuchos::RCP<MV> MX,
00459           Teuchos::Array<Teuchos::RCP<const MV> > MQ
00460       ) const 
00461   {
00462     projectGen(X,Q,Q,true,C,MX,MQ,MQ);
00463   }
00464 
00465 
00466 
00468   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00469   template<class ScalarType, class MV, class OP>
00470   int ICGSOrthoManager<ScalarType, MV, OP>::normalizeMat(
00471         MV &X, 
00472         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00473         Teuchos::RCP<MV> MX) const 
00474   {
00475     // call findBasis(), with the instruction to try to generate a basis of rank numvecs(X)
00476     // findBasis() requires MX
00477 
00478     int xc = MVT::GetNumberVecs(X);
00479     ptrdiff_t xr = MVText::GetGlobalLength(X);
00480 
00481     // if Op==null, MX == X (via pointer)
00482     // Otherwise, either the user passed in MX or we will allocated and compute it
00483     if (this->_hasOp) {
00484       if (MX == Teuchos::null) {
00485         // we need to allocate space for MX
00486         MX = MVT::Clone(X,xc);
00487         OPT::Apply(*(this->_Op),X,*MX);
00488         this->_OpCounter += MVT::GetNumberVecs(X);
00489       }
00490     }
00491 
00492     // if the user doesn't want to store the coefficients, 
00493     // allocate some local memory for them 
00494     if ( B == Teuchos::null ) {
00495       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00496     }
00497 
00498     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00499     ptrdiff_t mxr = (this->_hasOp) ? MVText::GetGlobalLength( *MX )  : xr;
00500 
00501     // check size of C, B
00502     TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
00503                         "Anasazi::ICGSOrthoManager::normalizeMat(): X must be non-empty" );
00504     TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00505                         "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
00506     TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
00507                         "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
00508     TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument, 
00509                         "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
00510 
00511     return findBasis(X, MX, *B, true );
00512   }
00513 
00514 
00515 
00517   // Find an Op-orthonormal basis for span(X) - span(W)
00518   template<class ScalarType, class MV, class OP>
00519   int ICGSOrthoManager<ScalarType, MV, OP>::projectAndNormalizeMat(
00520           MV &X, 
00521           Teuchos::Array<Teuchos::RCP<const MV> >  Q,
00522           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00523           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00524           Teuchos::RCP<MV> MX,
00525           Teuchos::Array<Teuchos::RCP<const MV> > MQ
00526       ) const 
00527   {
00528     return projectAndNormalizeGen(X,Q,Q,true,C,B,MX,MQ,MQ);
00529   }
00530 
00531 
00532 
00534   template<class ScalarType, class MV, class OP>
00535   void ICGSOrthoManager<ScalarType, MV, OP>::projectGen(
00536           MV &S,
00537           Teuchos::Array<Teuchos::RCP<const MV> > X,
00538           Teuchos::Array<Teuchos::RCP<const MV> > Y,
00539           bool isBiortho,
00540           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00541           Teuchos::RCP<MV> MS,
00542           Teuchos::Array<Teuchos::RCP<const MV> > MX,
00543           Teuchos::Array<Teuchos::RCP<const MV> > MY
00544       ) const 
00545   {
00546     // For the inner product defined by the operator Op or the identity (Op == 0)
00547     //   -> Orthogonalize S against each Y[i], modifying it in the range of X[i]
00548     // Modify MS accordingly
00549     //
00550     // Note that when Op is 0, MS is not referenced
00551     //
00552     // Parameter variables
00553     //
00554     // S  : Multivector to be transformed
00555     //
00556     // MS : Image of the block vector S by the mass matrix
00557     //
00558     // X,Y: Bases to orthogonalize against/via.
00559     //
00560 
00561 #ifdef ANASAZI_ICGS_DEBUG
00562     // Get a FancyOStream from out_arg or create a new one ...
00563     Teuchos::RCP<Teuchos::FancyOStream>
00564       out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
00565     out->setShowAllFrontMatter(false).setShowProcRank(true);
00566     *out << "Entering Anasazi::ICGSOrthoManager::projectGen(...)\n";
00567 #endif
00568 
00569     const ScalarType     ONE = SCT::one();
00570     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
00571     Teuchos::LAPACK<int,ScalarType> lapack;
00572     Teuchos::BLAS<int,ScalarType> blas;
00573 
00574     int sc = MVT::GetNumberVecs( S );
00575     ptrdiff_t sr = MVText::GetGlobalLength( S );
00576     int numxy = X.length();
00577     TEUCHOS_TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument, 
00578         "Anasazi::ICGSOrthoManager::projectGen(): X and Y must contain the same number of multivectors.");
00579     std::vector<int> xyc(numxy);
00580     // short-circuit
00581     if (numxy == 0 || sc == 0 || sr == 0) {
00582 #ifdef ANASAZI_ICGS_DEBUG
00583       *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
00584 #endif
00585       return;
00586     }
00587     // if we don't have enough C, expand it with null references
00588     // if we have too many, resize to throw away the latter ones
00589     // if we have exactly as many as we have X,Y this call has no effect
00590     //
00591     // same for MX, MY
00592     C.resize(numxy);
00593     MX.resize(numxy);
00594     MY.resize(numxy);
00595 
00596     // check size of S w.r.t. common sense
00597     TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument, 
00598                         "Anasazi::ICGSOrthoManager::projectGen(): MVT returned negative dimensions for S." );
00599 
00600     // check size of MS
00601     if (this->_hasOp == true) {
00602       if (MS != Teuchos::null) {
00603         TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MS) != sr, std::invalid_argument, 
00604             "Anasazi::ICGSOrthoManager::projectGen(): MS length not consistent with S." );
00605         TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument,
00606             "Anasazi::ICGSOrthoManager::projectGen(): MS width not consistent with S." );
00607       }
00608     }
00609 
00610     // tally up size of all X,Y and check/allocate C
00611     ptrdiff_t sumxyc = 0;
00612     int MYmissing = 0;
00613     int MXmissing = 0;
00614     for (int i=0; i<numxy; i++) {
00615       if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
00616         TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*X[i]) != sr, std::invalid_argument, 
00617             "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "] length not consistent with S." );
00618         TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*Y[i]) != sr, std::invalid_argument, 
00619             "Anasazi::ICGSOrthoManager::projectGen(): Y[" << i << "] length not consistent with S." );
00620         TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument,
00621             "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "] and Y[" << i << "] widths not consistent." );
00622 
00623         xyc[i] = MVT::GetNumberVecs( *X[i] );
00624         TEUCHOS_TEST_FOR_EXCEPTION( sr < static_cast<ptrdiff_t>(xyc[i]), std::invalid_argument, 
00625             "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "],Y[" << i << "] have less rows than columns, and therefore cannot be full rank." );
00626         sumxyc += xyc[i];
00627 
00628         // check size of C[i]
00629         if ( C[i] == Teuchos::null ) {
00630           C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) );
00631         }
00632         else {
00633           TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument, 
00634               "Anasazi::ICGSOrthoManager::projectGen(): Size of Q not consistent with size of C." );
00635         }
00636         // check sizes of MX[i], MY[i] if present
00637         // if not present, count their absence
00638         if (MX[i] != Teuchos::null) {
00639           TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
00640               "Anasazi::ICGSOrthoManager::projectGen(): Size of MX[" << i << "] not consistent with size of X[" << i << "]." );
00641         }
00642         else {
00643           MXmissing += xyc[i];
00644         }
00645         if (MY[i] != Teuchos::null) {
00646           TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
00647               "Anasazi::ICGSOrthoManager::projectGen(): Size of MY[" << i << "] not consistent with size of Y[" << i << "]." );
00648         }
00649         else {
00650           MYmissing += xyc[i];
00651         }
00652       }
00653       else {
00654         // if one is null and the other is not... the user may have made a mistake
00655         TEUCHOS_TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument,
00656             "Anasazi::ICGSOrthoManager::projectGen(): " 
00657             << (X[i] == Teuchos::null ? "Y[" : "X[") << i << "] was provided but " 
00658             << (X[i] == Teuchos::null ? "X[" : "Y[") << i << "] was not.");
00659       }
00660     }
00661 
00662     // is this operation even feasible?
00663     TEUCHOS_TEST_FOR_EXCEPTION(sumxyc > sr, std::invalid_argument, 
00664         "Anasazi::ICGSOrthoManager::projectGen(): dimension of all X[i],Y[i] is "
00665         << sumxyc << ", but length of vectors is only " << sr << ". This is infeasible.");
00666 
00667 
00668     /******   DO NO MODIFY *MS IF _hasOp == false   
00669      * if _hasOp == false, we don't need MS, MX or MY
00670      *
00671      * if _hasOp == true, we need MS (for S M-norms) and
00672      * therefore, we must also update MS, regardless of whether
00673      * it gets returned to the user (i.e., whether the user provided it)
00674      * we may need to allocate and compute MX or MY
00675      *
00676      * let BXY denote:
00677      *    "X" - we have all M*X[i]
00678      *    "Y" - we have all M*Y[i]
00679      *    "B" - we have biorthogonality (for all X[i],Y[i])
00680      * Encode these as values 
00681      *   X = 1
00682      *   Y = 2
00683      *   B = 4
00684      * We have 8 possibilities, 0-7
00685      *
00686      * We must allocate storage and computer the following, lest 
00687      * innerProdMat do it for us:
00688      *  none (0) - allocate MX, for inv(<X,Y>) and M*S
00689      *
00690      * for the following, we will compute M*S manually before returning
00691      *   B(4)  BY(6)  Y(2)                    -->  updateMS = 1
00692      * for the following, we will update M*S as we go, using M*X
00693      *   XY(3)  X(1)  none(0)  BXY(7)  BX(5)  -->  updateMS = 2
00694      * these choices favor applications of M over allocation of memory
00695      *
00696      */
00697     int updateMS = -1;
00698     if (this->_hasOp) {
00699       int whichAlloc = 0;
00700       if (MXmissing == 0) {
00701         whichAlloc += 1;
00702       }
00703       if (MYmissing == 0) {
00704         whichAlloc += 2;
00705       }
00706       if (isBiortho) {
00707         whichAlloc += 4;
00708       }
00709 
00710       switch (whichAlloc) {
00711         case 2:
00712         case 4:
00713         case 6:
00714           updateMS = 1;
00715           break;
00716         case 0:
00717         case 1:
00718         case 3:
00719         case 5:
00720         case 7:
00721           updateMS = 2;
00722           break;
00723       }
00724 
00725       // produce MS
00726       if (MS == Teuchos::null) {
00727 #ifdef ANASAZI_ICGS_DEBUG
00728     *out << "Allocating MS...\n";
00729 #endif
00730         MS = MVT::Clone(S,MVT::GetNumberVecs(S));
00731         OPT::Apply(*(this->_Op),S,*MS);
00732         this->_OpCounter += MVT::GetNumberVecs(S);
00733       }
00734 
00735       // allocate the rest
00736       if (whichAlloc == 0) {
00737         // allocate and compute missing MX
00738         for (int i=0; i<numxy; i++) {
00739           if (MX[i] == Teuchos::null) {
00740 #ifdef ANASAZI_ICGS_DEBUG
00741             *out << "Allocating MX[" << i << "]...\n";
00742 #endif
00743             Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]);
00744             OPT::Apply(*(this->_Op),*X[i],*tmpMX);
00745             MX[i] = tmpMX;
00746             this->_OpCounter += xyc[i];
00747           }
00748         }
00749       }
00750     }
00751     else {
00752       // Op == I  -->  MS == S 
00753       MS = Teuchos::rcpFromRef(S);
00754       updateMS = 0;
00755     }
00756     TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error,
00757         "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
00758 
00759 
00761     // Perform the Gram-Schmidt transformation for a block of vectors
00763 
00764     // Compute Cholesky factorizations for the Y'*M*X
00765     // YMX stores the YMX (initially) and their Cholesky factorizations (utlimately)
00766     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > YMX(numxy);
00767     if (isBiortho == false) {
00768       for (int i=0; i<numxy; i++) {
00769 #ifdef ANASAZI_ICGS_DEBUG
00770         *out << "Computing YMX[" << i << "] and its Cholesky factorization...\n";
00771 #endif
00772         YMX[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],xyc[i]) );
00773         MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Y[i],*X[i],*YMX[i],MY[i],MX[i]);
00774 #ifdef ANASAZI_ICGS_DEBUG
00775         // YMX should be symmetric positive definite
00776         // the cholesky will check the positive definiteness, but it looks only at the upper half
00777         // we will check the symmetry by removing the upper half from the lower half, which should 
00778         // result in zeros
00779         // also, diagonal of YMX should be real; consider the complex part as error
00780         {
00781           MagnitudeType err = ZERO;
00782           for (int jj=0; jj<YMX[i]->numCols(); jj++) {
00783             err =+ SCT::magnitude(SCT::imag((*YMX[i])(jj,jj)));
00784             for (int ii=jj; ii<YMX[i]->numRows(); ii++) {
00785               err += SCT::magnitude( (*YMX[i])(ii,jj) - SCT::conjugate((*YMX[i])(jj,ii)) );
00786             }
00787           }
00788           *out << "Symmetry error in YMX[" << i << "] == " << err << "\n";
00789         }
00790 #endif
00791         // take the LU factorization
00792         int info;
00793         lapack.POTRF('U',YMX[i]->numRows(),YMX[i]->values(),YMX[i]->stride(),&info);
00794         TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
00795             "Anasazi::ICGSOrthoManager::projectGen(): Error computing Cholesky factorization of Y[i]^T * M * X[i] using POTRF: returned info " << info);
00796       }
00797     }
00798 
00799     // Compute the initial Op-norms
00800 #ifdef ANASAZI_ICGS_DEBUG
00801     std::vector<MagnitudeType> oldNorms(sc);
00802     MatOrthoManager<ScalarType,MV,OP>::normMat(S,oldNorms,MS);
00803     *out << "oldNorms = { ";
00804     std::copy(oldNorms.begin(), oldNorms.end(), std::ostream_iterator<MagnitudeType>(*out, " "));
00805     *out << "}\n";
00806 #endif
00807 
00808 
00809     // clear the C[i] and allocate Ccur
00810     Teuchos::Array<Teuchos::SerialDenseMatrix<int,ScalarType> > Ccur(numxy);
00811     for (int i=0; i<numxy; i++) {
00812       C[i]->putScalar(ZERO);
00813       Ccur[i].reshape(C[i]->numRows(),C[i]->numCols());
00814     }
00815 
00816     for (int iter=0; iter < numIters_; iter++) {
00817 #ifdef ANASAZI_ICGS_DEBUG
00818       *out << "beginning iteration " << iter+1 << "\n";
00819 #endif
00820 
00821       // Define the product Y[i]'*Op*S
00822       for (int i=0; i<numxy; i++) {
00823         // Compute Y[i]'*M*S
00824         MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Y[i],S,Ccur[i],MY[i],MS);
00825         if (isBiortho == false) {
00826           // C[i] = inv(YMX[i])*(Y[i]'*M*S)
00827           int info;
00828           lapack.POTRS('U',YMX[i]->numCols(),Ccur[i].numCols(),
00829               YMX[i]->values(),YMX[i]->stride(),
00830               Ccur[i].values(),Ccur[i].stride(), &info);
00831           TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
00832               "Anasazi::ICGSOrthoManager::projectGen(): Error code " << info << " from lapack::POTRS." );
00833         }
00834 
00835         // Multiply by X[i] and subtract the result in S
00836 #ifdef ANASAZI_ICGS_DEBUG
00837         *out << "Applying projector P_{X[" << i << "],Y[" << i << "]}...\n";
00838 #endif
00839         MVT::MvTimesMatAddMv( -ONE, *X[i], Ccur[i], ONE, S );
00840 
00841         // Accumulate coeffs into previous step
00842         *C[i] += Ccur[i];
00843 
00844         // Update MS as required
00845         if (updateMS == 1) {
00846 #ifdef ANASAZI_ICGS_DEBUG
00847           *out << "Updating MS...\n";
00848 #endif
00849           OPT::Apply( *(this->_Op), S, *MS);
00850           this->_OpCounter += sc;
00851         }
00852         else if (updateMS == 2) {
00853 #ifdef ANASAZI_ICGS_DEBUG
00854           *out << "Updating MS...\n";
00855 #endif
00856           MVT::MvTimesMatAddMv( -ONE, *MX[i], Ccur[i], ONE, *MS );
00857         }
00858       }
00859 
00860       // Compute new Op-norms
00861 #ifdef ANASAZI_ICGS_DEBUG
00862       std::vector<MagnitudeType> newNorms(sc);
00863       MatOrthoManager<ScalarType,MV,OP>::normMat(S,newNorms,MS);
00864       *out << "newNorms = { ";
00865       std::copy(newNorms.begin(), newNorms.end(), std::ostream_iterator<MagnitudeType>(*out, " "));
00866       *out << "}\n";
00867 #endif
00868     }
00869 
00870 #ifdef ANASAZI_ICGS_DEBUG
00871     *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
00872 #endif
00873   }
00874 
00875 
00876 
00878   // Find an Op-orthonormal basis for span(S) - span(Y)
00879   template<class ScalarType, class MV, class OP>
00880   int ICGSOrthoManager<ScalarType, MV, OP>::projectAndNormalizeGen(
00881           MV &S,
00882           Teuchos::Array<Teuchos::RCP<const MV> > X,
00883           Teuchos::Array<Teuchos::RCP<const MV> > Y,
00884           bool isBiortho,
00885           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00886           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00887           Teuchos::RCP<MV> MS,
00888           Teuchos::Array<Teuchos::RCP<const MV> > MX,
00889           Teuchos::Array<Teuchos::RCP<const MV> > MY
00890       ) const {
00891     // For the inner product defined by the operator Op or the identity (Op == 0)
00892     //   -> Orthogonalize S against each Y[i], modifying it in the range of X[i]
00893     // Modify MS accordingly
00894     // Then construct a M-orthonormal basis for the remaining part
00895     //
00896     // Note that when Op is 0, MS is not referenced
00897     //
00898     // Parameter variables
00899     //
00900     // S  : Multivector to be transformed
00901     //
00902     // MS : Image of the block vector S by the mass matrix
00903     //
00904     // X,Y: Bases to orthogonalize against/via.
00905     //
00906 
00907 #ifdef ANASAZI_ICGS_DEBUG
00908     // Get a FancyOStream from out_arg or create a new one ...
00909     Teuchos::RCP<Teuchos::FancyOStream>
00910       out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
00911     out->setShowAllFrontMatter(false).setShowProcRank(true);
00912     *out << "Entering Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
00913 #endif
00914 
00915     int rank;
00916     int sc = MVT::GetNumberVecs( S );
00917     ptrdiff_t sr = MVText::GetGlobalLength( S );
00918     int numxy = X.length();
00919     TEUCHOS_TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument, 
00920         "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X and Y must contain the same number of multivectors.");
00921     std::vector<int> xyc(numxy);
00922     // short-circuit
00923     if (sc == 0 || sr == 0) {
00924 #ifdef ANASAZI_ICGS_DEBUG
00925       *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
00926 #endif
00927       return 0;
00928     }
00929     // if we don't have enough C, expand it with null references
00930     // if we have too many, resize to throw away the latter ones
00931     // if we have exactly as many as we have X,Y this call has no effect
00932     //
00933     // same for MX, MY
00934     C.resize(numxy);
00935     MX.resize(numxy);
00936     MY.resize(numxy);
00937 
00938     // check size of S w.r.t. common sense
00939     TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument, 
00940                         "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MVT returned negative dimensions for S." );
00941 
00942     // check size of MS
00943     if (this->_hasOp == true) {
00944       if (MS != Teuchos::null) {
00945         TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MS) != sr, std::invalid_argument, 
00946             "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS length not consistent with S." );
00947         TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument,
00948             "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS width not consistent with S." );
00949       }
00950     }
00951 
00952     // tally up size of all X,Y and check/allocate C
00953     ptrdiff_t sumxyc = 0;
00954     int MYmissing = 0;
00955     int MXmissing = 0;
00956     for (int i=0; i<numxy; i++) {
00957       if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
00958         TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*X[i]) != sr, std::invalid_argument, 
00959             "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "] length not consistent with S." );
00960         TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*Y[i]) != sr, std::invalid_argument, 
00961             "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Y[" << i << "] length not consistent with S." );
00962         TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument,
00963             "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "] and Y[" << i << "] widths not consistent." );
00964 
00965         xyc[i] = MVT::GetNumberVecs( *X[i] );
00966         TEUCHOS_TEST_FOR_EXCEPTION( sr < static_cast<ptrdiff_t>(xyc[i]), std::invalid_argument, 
00967             "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "],Y[" << i << "] have less rows than columns, and therefore cannot be full rank." );
00968         sumxyc += xyc[i];
00969 
00970         // check size of C[i]
00971         if ( C[i] == Teuchos::null ) {
00972           C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) );
00973         }
00974         else {
00975           TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument, 
00976               "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of Q not consistent with size of C." );
00977         }
00978         // check sizes of MX[i], MY[i] if present
00979         // if not present, count their absence
00980         if (MX[i] != Teuchos::null) {
00981           TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
00982               "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MX[" << i << "] not consistent with size of X[" << i << "]." );
00983         }
00984         else {
00985           MXmissing += xyc[i];
00986         }
00987         if (MY[i] != Teuchos::null) {
00988           TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
00989               "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MY[" << i << "] not consistent with size of Y[" << i << "]." );
00990         }
00991         else {
00992           MYmissing += xyc[i];
00993         }
00994       }
00995       else {
00996         // if one is null and the other is not... the user may have made a mistake
00997         TEUCHOS_TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument,
00998             "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): " 
00999             << (X[i] == Teuchos::null ? "Y[" : "X[") << i << "] was provided but " 
01000             << (X[i] == Teuchos::null ? "X[" : "Y[") << i << "] was not.");
01001       }
01002     }
01003 
01004     // is this operation even feasible?
01005     TEUCHOS_TEST_FOR_EXCEPTION(sumxyc + sc > sr, std::invalid_argument, 
01006         "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): dimension of all X[i],Y[i] is "
01007         << sumxyc << " and requested " << sc << "-dimensional basis, but length of vectors is only " 
01008         << sr << ". This is infeasible.");
01009 
01010 
01011     /******   DO NO MODIFY *MS IF _hasOp == false   
01012      * if _hasOp == false, we don't need MS, MX or MY
01013      *
01014      * if _hasOp == true, we need MS (for S M-norms and normalization) and
01015      * therefore, we must also update MS, regardless of whether
01016      * it gets returned to the user (i.e., whether the user provided it)
01017      * we may need to allocate and compute MX or MY
01018      *
01019      * let BXY denote:
01020      *    "X" - we have all M*X[i]
01021      *    "Y" - we have all M*Y[i]
01022      *    "B" - we have biorthogonality (for all X[i],Y[i])
01023      * Encode these as values 
01024      *   X = 1
01025      *   Y = 2
01026      *   B = 4
01027      * We have 8 possibilities, 0-7
01028      *
01029      * We must allocate storage and computer the following, lest 
01030      * innerProdMat do it for us:
01031      *  none (0) - allocate MX, for inv(<X,Y>) and M*S
01032      *
01033      * for the following, we will compute M*S manually before returning
01034      *   B(4)  BY(6)  Y(2)                    -->  updateMS = 1
01035      * for the following, we will update M*S as we go, using M*X
01036      *   XY(3)  X(1)  none(0)  BXY(7)  BX(5)  -->  updateMS = 2
01037      * these choices favor applications of M over allocation of memory
01038      *
01039      */
01040     int updateMS = -1;
01041     if (this->_hasOp) {
01042       int whichAlloc = 0;
01043       if (MXmissing == 0) {
01044         whichAlloc += 1;
01045       }
01046       if (MYmissing == 0) {
01047         whichAlloc += 2;
01048       }
01049       if (isBiortho) {
01050         whichAlloc += 4;
01051       }
01052 
01053       switch (whichAlloc) {
01054         case 2:
01055         case 4:
01056         case 6:
01057           updateMS = 1;
01058           break;
01059         case 0:
01060         case 1:
01061         case 3:
01062         case 5:
01063         case 7:
01064           updateMS = 2;
01065           break;
01066       }
01067 
01068       // produce MS
01069       if (MS == Teuchos::null) {
01070 #ifdef ANASAZI_ICGS_DEBUG
01071     *out << "Allocating MS...\n";
01072 #endif
01073         MS = MVT::Clone(S,MVT::GetNumberVecs(S));
01074         OPT::Apply(*(this->_Op),S,*MS);
01075         this->_OpCounter += MVT::GetNumberVecs(S);
01076       }
01077 
01078       // allocate the rest
01079       if (whichAlloc == 0) {
01080         // allocate and compute missing MX
01081         for (int i=0; i<numxy; i++) {
01082           if (MX[i] == Teuchos::null) {
01083 #ifdef ANASAZI_ICGS_DEBUG
01084             *out << "Allocating MX[" << i << "]...\n";
01085 #endif
01086             Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]);
01087             OPT::Apply(*(this->_Op),*X[i],*tmpMX);
01088             MX[i] = tmpMX;
01089             this->_OpCounter += xyc[i];
01090           }
01091         }
01092       }
01093     }
01094     else {
01095       // Op == I  -->  MS == S 
01096       MS = Teuchos::rcpFromRef(S);
01097       updateMS = 0;
01098     }
01099     TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error,
01100         "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
01101 
01102 
01103     // if the user doesn't want to store the coefficients, 
01104     // allocate some local memory for them 
01105     if ( B == Teuchos::null ) {
01106       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(sc,sc) );
01107     }
01108     else {
01109       // check size of B
01110       TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != sc || B->numCols() != sc, std::invalid_argument, 
01111           "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of S must be consistent with size of B" );
01112     }
01113 
01114 
01115     // orthogonalize all of S against X,Y
01116     projectGen(S,X,Y,isBiortho,C,MS,MX,MY);
01117 
01118     Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(sc,1);
01119     // start working
01120     rank = 0;
01121     int numTries = 10;   // each vector in X gets 10 random chances to escape degeneracy
01122     int oldrank = -1;
01123     do {
01124       int curssize = sc - rank;
01125 
01126       // orthonormalize S, but quit if it is rank deficient
01127       // we can't let findBasis generated random vectors to complete the basis,
01128       // because it doesn't know about Q; we will do this ourselves below
01129 #ifdef ANASAZI_ICGS_DEBUG
01130       *out << "Attempting to find orthonormal basis for X...\n";
01131 #endif
01132       rank = findBasis(S,MS,*B,false,curssize);
01133 
01134       if (oldrank != -1 && rank != oldrank) {
01135         // we had previously stopped before, after operating on vector oldrank
01136         // we saved its coefficients, augmented it with a random vector, and
01137         // then called findBasis() again, which proceeded to add vector oldrank
01138         // to the basis. 
01139         // now, restore the saved coefficients into B
01140         for (int i=0; i<sc; i++) {
01141           (*B)(i,oldrank) = oldCoeff(i,0);
01142         }
01143       }
01144 
01145       if (rank < sc) {
01146         if (rank != oldrank) {
01147           // we quit on this vector and will augment it with random below
01148           // this is the first time that we have quit on this vector
01149           // therefor, (*B)(:,rank) contains the actual coefficients of the 
01150           // input vectors with respect to the previous vectors in the basis
01151           // save these values, as (*B)(:,rank) will be overwritten by our next
01152           // call to findBasis()
01153           // we will restore it after we are done working on this vector
01154           for (int i=0; i<sc; i++) {
01155             oldCoeff(i,0) = (*B)(i,rank);
01156           }
01157         }
01158       }
01159 
01160       if (rank == sc) {
01161         // we are done
01162 #ifdef ANASAZI_ICGS_DEBUG
01163         *out << "Finished computing basis.\n";
01164 #endif
01165         break;
01166       }
01167       else {
01168         TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank, OrthoError,   
01169                             "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): basis lost rank; this shouldn't happen");
01170 
01171         if (rank != oldrank) {
01172           // we added a basis vector from random info; reset the chance counter
01173           numTries = 10;
01174           // store old rank
01175           oldrank = rank;
01176         }
01177         else {
01178           // has this vector run out of chances to escape degeneracy?
01179           if (numTries <= 0) {
01180             break;
01181           }
01182         }
01183         // use one of this vector's chances
01184         numTries--;
01185 
01186         // randomize troubled direction
01187 #ifdef ANASAZI_ICGS_DEBUG
01188         *out << "Inserting random vector in X[" << rank << "]. Attempt " << 10-numTries << ".\n";
01189 #endif
01190         Teuchos::RCP<MV> curS, curMS;
01191         std::vector<int> ind(1);
01192         ind[0] = rank;
01193         curS  = MVT::CloneViewNonConst(S,ind);
01194         MVT::MvRandom(*curS);
01195         if (this->_hasOp) {
01196 #ifdef ANASAZI_ICGS_DEBUG
01197           *out << "Applying operator to random vector.\n";
01198 #endif
01199           curMS = MVT::CloneViewNonConst(*MS,ind);
01200           OPT::Apply( *(this->_Op), *curS, *curMS );
01201           this->_OpCounter += MVT::GetNumberVecs(*curS);
01202         }
01203 
01204         // orthogonalize against X,Y
01205         // if !this->_hasOp, the curMS will be ignored.
01206         // we don't care about these coefficients
01207         // on the contrary, we need to preserve the previous coeffs
01208         {
01209           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
01210           projectGen(*curS,X,Y,isBiortho,dummyC,curMS,MX,MY);
01211         }
01212       }
01213     } while (1);
01214 
01215     // this should never raise an exception; but our post-conditions oblige us to check
01216     TEUCHOS_TEST_FOR_EXCEPTION( rank > sc || rank < 0, std::logic_error, 
01217                         "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Debug error in rank variable." );
01218 
01219 #ifdef ANASAZI_ICGS_DEBUG
01220     *out << "Returning " << rank << " from Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
01221 #endif
01222 
01223     return rank;
01224   }
01225 
01226 
01227 
01229   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
01230   // the rank is numvectors(X)
01231   template<class ScalarType, class MV, class OP>
01232   int ICGSOrthoManager<ScalarType, MV, OP>::findBasis(
01233                 MV &X, Teuchos::RCP<MV> MX, 
01234                 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
01235                 bool completeBasis, int howMany ) const {
01236 
01237     // For the inner product defined by the operator Op or the identity (Op == 0)
01238     //   -> Orthonormalize X 
01239     // Modify MX accordingly
01240     //
01241     // Note that when Op is 0, MX is not referenced
01242     //
01243     // Parameter variables
01244     //
01245     // X  : Vectors to be orthonormalized
01246     //
01247     // MX : Image of the multivector X under the operator Op
01248     //
01249     // Op  : Pointer to the operator for the inner product
01250     //
01251 
01252 #ifdef ANASAZI_ICGS_DEBUG
01253     // Get a FancyOStream from out_arg or create a new one ...
01254     Teuchos::RCP<Teuchos::FancyOStream>
01255       out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
01256     out->setShowAllFrontMatter(false).setShowProcRank(true);
01257     *out << "Entering Anasazi::ICGSOrthoManager::findBasis(...)\n";
01258 #endif
01259 
01260     const ScalarType     ONE = SCT::one();
01261     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
01262 
01263     int xc = MVT::GetNumberVecs( X );
01264 
01265     if (howMany == -1) {
01266       howMany = xc;
01267     }
01268 
01269     /*******************************************************
01270      *  If _hasOp == false, we will not reference MX below *
01271      *******************************************************/
01272     TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp == true && MX == Teuchos::null, std::logic_error,
01273         "Anasazi::ICGSOrthoManager::findBasis(): calling routine did not specify MS.");
01274     TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error, 
01275                         "Anasazi::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
01276 
01277     /* xstart is which column we are starting the process with, based on howMany
01278      * columns before xstart are assumed to be Op-orthonormal already
01279      */
01280     int xstart = xc - howMany;
01281 
01282     for (int j = xstart; j < xc; j++) {
01283 
01284       // numX represents the number of currently orthonormal columns of X
01285       int numX = j;
01286       // j represents the index of the current column of X
01287       // these are different interpretations of the same value
01288 
01289       // 
01290       // set the lower triangular part of B to zero
01291       for (int i=j+1; i<xc; ++i) {
01292         B(i,j) = ZERO;
01293       }
01294 
01295       // Get a view of the vector currently being worked on.
01296       std::vector<int> index(1);
01297       index[0] = j;
01298       Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
01299       Teuchos::RCP<MV> MXj;
01300       if ((this->_hasOp)) {
01301         // MXj is a view of the current vector in MX
01302         MXj = MVT::CloneViewNonConst( *MX, index );
01303       }
01304       else {
01305         // MXj is a pointer to Xj, and MUST NOT be modified
01306         MXj = Xj;
01307       }
01308 
01309       // Get a view of the previous vectors.
01310       std::vector<int> prev_idx( numX );
01311       Teuchos::RCP<const MV> prevX, prevMX;
01312 
01313       if (numX > 0) {
01314         for (int i=0; i<numX; ++i) prev_idx[i] = i;
01315         prevX = MVT::CloneView( X, prev_idx );
01316         if (this->_hasOp) {
01317           prevMX = MVT::CloneView( *MX, prev_idx );
01318         }
01319       } 
01320 
01321       bool rankDef = true;
01322       /* numTrials>0 will denote that the current vector was randomized for the purpose
01323        * of finding a basis vector, and that the coefficients of that vector should
01324        * not be stored in B
01325        */
01326       for (int numTrials = 0; numTrials < 10; numTrials++) {
01327 #ifdef ANASAZI_ICGS_DEBUG
01328         *out << "Trial " << numTrials << " for vector " << j << "\n";
01329 #endif
01330 
01331         // Make storage for these Gram-Schmidt iterations.
01332         Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
01333         std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
01334 
01335         //
01336         // Save old MXj vector and compute Op-norm
01337         //
01338         Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
01339         MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,origNorm,MXj);
01340 #ifdef ANASAZI_ICGS_DEBUG
01341         *out << "origNorm = " << origNorm[0] << "\n";
01342 #endif
01343 
01344         if (numX > 0) {
01345           // Apply the first step of Gram-Schmidt
01346 
01347           // product <- prevX^T MXj
01348           MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,product,Teuchos::null,MXj);
01349 
01350           // Xj <- Xj - prevX prevX^T MXj   
01351           //     = Xj - prevX product
01352 #ifdef ANASAZI_ICGS_DEBUG
01353           *out << "Orthogonalizing X[" << j << "]...\n";
01354 #endif
01355           MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
01356 
01357           // Update MXj
01358           if (this->_hasOp) {
01359             // MXj <- Op*Xj_new
01360             //      = Op*(Xj_old - prevX prevX^T MXj)
01361             //      = MXj - prevMX product
01362 #ifdef ANASAZI_ICGS_DEBUG
01363             *out << "Updating MX[" << j << "]...\n";
01364 #endif
01365             MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
01366           }
01367 
01368           // Compute new Op-norm
01369           MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm,MXj);
01370           MagnitudeType product_norm = product.normOne();
01371           
01372 #ifdef ANASAZI_ICGS_DEBUG
01373           *out << "newNorm = " << newNorm[0] << "\n";
01374           *out << "prodoct_norm = " << product_norm << "\n";
01375 #endif
01376 
01377           // Check if a correction is needed.
01378           if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
01379 #ifdef ANASAZI_ICGS_DEBUG
01380             if (product_norm/newNorm[0] >= tol_) {
01381               *out << "product_norm/newNorm == " << product_norm/newNorm[0] << "... another step of Gram-Schmidt.\n";
01382             }
01383             else {
01384               *out << "eps*origNorm == " << eps_*origNorm[0] << "... another step of Gram-Schmidt.\n";
01385             }
01386 #endif
01387             // Apply the second step of Gram-Schmidt
01388             // This is the same as above
01389             Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
01390             MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,P2,Teuchos::null,MXj);
01391             product += P2;
01392 #ifdef ANASAZI_ICGS_DEBUG
01393             *out << "Orthogonalizing X[" << j << "]...\n";
01394 #endif
01395             MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
01396             if ((this->_hasOp)) {
01397 #ifdef ANASAZI_ICGS_DEBUG
01398               *out << "Updating MX[" << j << "]...\n";
01399 #endif
01400               MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
01401             }
01402             // Compute new Op-norms
01403             MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm2,MXj);
01404             product_norm = P2.normOne();
01405 #ifdef ANASAZI_ICGS_DEBUG
01406             *out << "newNorm2 = " << newNorm2[0] << "\n";
01407             *out << "product_norm = " << product_norm << "\n";
01408 #endif
01409             if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
01410               // we don't do another GS, we just set it to zero.
01411 #ifdef ANASAZI_ICGS_DEBUG
01412               if (product_norm/newNorm2[0] >= tol_) {
01413                 *out << "product_norm/newNorm2 == " << product_norm/newNorm2[0] << "... setting vector to zero.\n";
01414               }
01415               else if (newNorm[0] < newNorm2[0]) {
01416                 *out << "newNorm2 > newNorm... setting vector to zero.\n";
01417               }
01418               else {
01419                 *out << "eps*origNorm == " << eps_*origNorm[0] << "... setting vector to zero.\n";
01420               }
01421 #endif
01422               MVT::MvInit(*Xj,ZERO);
01423               if ((this->_hasOp)) {
01424 #ifdef ANASAZI_ICGS_DEBUG
01425                 *out << "Setting MX[" << j << "] to zero as well.\n";
01426 #endif
01427                 MVT::MvInit(*MXj,ZERO);
01428               }
01429             }
01430           } 
01431         } // if (numX > 0) do GS
01432 
01433         // save the coefficients, if we are working on the original vector and not a randomly generated one
01434         if (numTrials == 0) {
01435           for (int i=0; i<numX; i++) {
01436             B(i,j) = product(i,0);
01437           }
01438         }
01439 
01440         // Check if Xj has any directional information left after the orthogonalization.
01441         MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm,MXj);
01442         if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
01443 #ifdef ANASAZI_ICGS_DEBUG
01444           *out << "Normalizing X[" << j << "], norm(X[" << j << "]) = " << newNorm[0] << "\n";
01445 #endif
01446           // Normalize Xj.
01447           // Xj <- Xj / norm
01448           MVT::MvScale( *Xj, ONE/newNorm[0]);
01449           if (this->_hasOp) {
01450 #ifdef ANASAZI_ICGS_DEBUG
01451             *out << "Normalizing M*X[" << j << "]...\n";
01452 #endif
01453             // Update MXj.
01454             MVT::MvScale( *MXj, ONE/newNorm[0]);
01455           }
01456 
01457           // save it, if it corresponds to the original vector and not a randomly generated one
01458           if (numTrials == 0) {
01459             B(j,j) = newNorm[0];
01460           }
01461 
01462           // We are not rank deficient in this vector. Move on to the next vector in X.
01463           rankDef = false;
01464           break;
01465         }
01466         else {
01467 #ifdef ANASAZI_ICGS_DEBUG
01468           *out << "Not normalizing M*X[" << j << "]...\n";
01469 #endif
01470           // There was nothing left in Xj after orthogonalizing against previous columns in X.
01471           // X is rank deficient.
01472           // reflect this in the coefficients
01473           B(j,j) = ZERO;
01474 
01475           if (completeBasis) {
01476             // Fill it with random information and keep going.
01477 #ifdef ANASAZI_ICGS_DEBUG
01478             *out << "Inserting random vector in X[" << j << "]...\n";
01479 #endif
01480             MVT::MvRandom( *Xj );
01481             if (this->_hasOp) {
01482 #ifdef ANASAZI_ICGS_DEBUG
01483               *out << "Updating M*X[" << j << "]...\n";
01484 #endif
01485               OPT::Apply( *(this->_Op), *Xj, *MXj );
01486               this->_OpCounter += MVT::GetNumberVecs(*Xj);
01487             }
01488           }
01489           else {
01490             rankDef = true;
01491             break;
01492           }
01493         }
01494       }  // for (numTrials = 0; numTrials < 10; ++numTrials)
01495 
01496       // if rankDef == true, then quit and notify user of rank obtained
01497       if (rankDef == true) {
01498         TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError, 
01499                             "Anasazi::ICGSOrthoManager::findBasis(): Unable to complete basis" );
01500 #ifdef ANASAZI_ICGS_DEBUG
01501         *out << "Returning early, rank " << j << " from Anasazi::ICGSOrthoManager::findBasis(...)\n";
01502 #endif
01503         return j;
01504       }
01505 
01506     } // for (j = 0; j < xc; ++j)
01507 
01508 #ifdef ANASAZI_ICGS_DEBUG
01509     *out << "Returning " << xc << " from Anasazi::ICGSOrthoManager::findBasis(...)\n";
01510 #endif
01511     return xc;
01512   }
01513 
01514 } // namespace Anasazi
01515 
01516 #endif // ANASAZI_ICSG_ORTHOMANAGER_HPP
01517 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends