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