Anasazi Version of the Day
AnasaziTsqrOrthoManager.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2010) 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 
00031 
00032 #ifndef __AnasaziTsqrOrthoManager_hpp
00033 #define __AnasaziTsqrOrthoManager_hpp
00034 
00035 #include "AnasaziConfigDefs.hpp"
00036 #include "AnasaziMultiVecTraits.hpp"
00037 #include "AnasaziTsqrAdaptor.hpp"
00038 #include "AnasaziSVQBOrthoManager.hpp"
00039 #include "Teuchos_LAPACK.hpp"
00040 
00043 
00044 namespace Anasazi {
00045 
00048   class TsqrOrthoError : public OrthoError
00049   {
00050   public: 
00051     TsqrOrthoError (const std::string& what_arg) : 
00052       OrthoError (what_arg) {}
00053   };
00054 
00072   class TsqrOrthoFault : public OrthoError
00073   {
00074   public: 
00075     TsqrOrthoFault (const std::string& what_arg) : 
00076       OrthoError (what_arg) {}
00077   };
00078 
00079 
00080   class NonNullOperatorError : public OrthoError
00081   {
00082   public:
00083     NonNullOperatorError () : 
00084       OrthoError ("Sorry, TsqrOrthoManager doesn\'t work with a non-null Op "
00085       "argument.  I know this is bad class design, but it will "
00086       "have to do for now, since TsqrOrthoManager has to inherit "
00087       "from MatOrthoManager in order to work in Anasazi\'s "
00088       "eigensolvers.  If you want to solve this problem yourself, "
00089       "the thing to do is to have TsqrOrthoManager degrade to "
00090       "SVQBOrthoManager when this->getOp() != Teuchos::null.")
00091     {}
00092   };
00093 
00115   template< class ScalarType, class MV >
00116   class TsqrOrthoManagerImpl
00117   {
00118   private:
00119     typedef typename Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType;
00120     typedef Teuchos::ScalarTraits< ScalarType >    SCT;
00121     typedef Teuchos::ScalarTraits< MagnitudeType > SCTM;
00122     typedef MultiVecTraits< ScalarType, MV >       MVT;
00123 
00124     typedef typename Anasazi::TsqrAdaptor< ScalarType, MV >::adaptor_type tsqr_adaptor_type;
00125     typedef Teuchos::RCP< tsqr_adaptor_type > tsqr_adaptor_ptr;
00126 
00127   public:
00128 
00140     TsqrOrthoManagerImpl (const Teuchos::ParameterList& tsqrParams) :
00141       tsqrParams_ (tsqrParams),
00142       tsqrAdaptor_ (Teuchos::null),
00143       Q_ (Teuchos::null),
00144       eps_ (SCT::eps()), // ScalarTraits< ScalarType >::eps() returns MagnitudeType
00145       blockReorthogThreshold_ (MagnitudeType(1) / MagnitudeType(2)),
00146       relativeRankTolerance_ (MagnitudeType(100)*SCT::eps()),
00147       throwOnReorthogFault_ (true)
00148     {}
00149 
00150     // void 
00151     // setOp (Teuchos::RCP< const OP > Op)
00152     // {
00153     //   static_cast< MatOrthoManager< ScalarType, MV, OP >* const >(this)->setOp (Op);
00154     //   if (getOp() != Teuchos::null) 
00155     //  throw NonNullOperatorError();
00156     // }
00157 
00158     typedef Teuchos::RCP< MV >       mv_ptr;
00159     typedef Teuchos::RCP< const MV > const_mv_ptr;
00160     typedef Teuchos::Array< const_mv_ptr >                       const_prev_mvs_type;
00161     typedef Teuchos::SerialDenseMatrix< int, ScalarType >        serial_matrix_type;
00162     typedef Teuchos::RCP< serial_matrix_type >                   serial_matrix_ptr;
00163     typedef Teuchos::Array< Teuchos::RCP< serial_matrix_type > > prev_coeffs_type;
00164 
00165     void 
00166     innerProd (const MV& X, 
00167          const MV& Y, 
00168          serial_matrix_type& Z) const
00169     {
00170       MVT::MvTransMv (SCT::one(), X, Y, Z);
00171     }
00172 
00173     void
00174     norm (const MV& X, 
00175     std::vector< MagnitudeType >& normvec) const
00176     {
00177       const int ncols = MVT::GetNumberVecs (X);
00178 
00179       // mfh 20 Jul 2010: MatOrthoManager::normMat() computes norms
00180       // one column at a time, which results in too much
00181       // communication.  Instead, we compute the whole inner product,
00182       // which takes more computation but is the only option
00183       // available, given the current MVT interface.
00184       serial_matrix_type XTX (ncols, ncols); 
00185       innerProd (X, X, XTX);
00186 
00187       // Record the results.  Make sure normvec is long enough.
00188       if (normvec.size() < ncols)
00189   normvec.resize (ncols);
00190       for (int k = 0; k < ncols; ++k)
00191   normvec[k] = SCTM::squareroot (XTX(k,k));
00192     }
00193 
00199     void 
00200     project (MV& X, 
00201        const_prev_mvs_type Q,
00202        prev_coeffs_type C = Teuchos::tuple (serial_matrix_ptr (Teuchos::null))) const
00203     {
00204       int nrows_X, ncols_X, num_Q_blocks, ncols_Q_total;
00205       checkProjectionDims (nrows_X, ncols_X, num_Q_blocks, ncols_Q_total, X, Q);
00206       // Test for quick exit: any dimension of X is zero, or there are
00207       // zero Q blocks, or the total number of columns of the Q blocks
00208       // is zero.
00209       if (nrows_X == 0 || ncols_X == 0 || num_Q_blocks == 0 || ncols_Q_total == 0)
00210   return;
00211 
00212       // If we don't have enough C, expanding it creates null references.
00213       // If we have too many, resizing just throws away the later ones.
00214       // If we have exactly as many as we have Q, this call has no effect.
00215       C.resize (num_Q_blocks);
00216       for (int i = 0; i < num_Q_blocks; ++i) 
00217   {
00218     const int ncols_Q_i = MVT::GetNumberVecs (*Q[i]);
00219     // Create a new C[i] if necessary.
00220     if (C[i] == Teuchos::null)
00221       C[i] = Teuchos::rcp (new serial_matrix_type (ncols_Q_i, ncols_X));
00222   }
00223       rawProject (X, Q, C);
00224     }
00225 
00226 
00227     int 
00228     normalize (MV& X,
00229          serial_matrix_ptr B = Teuchos::null)
00230     {
00231       // Internal data used by this method require a specific MV
00232       // object for initialization (e.g., to get a Map / communicator,
00233       // and to initialize scratch space).  Thus, we delay (hence
00234       // "lazy") initialization until we get an X.
00235       lazyInit (X);
00236       
00237       // MVT returns int for this, even though the local_ordinal_type
00238       // of the MV may be some other type.
00239       const int ncols = MVT::GetNumberVecs (X);
00240 
00241       // TSQR's rank-revealing part doesn't work unless B is provided.
00242       // If B is not provided, allocate a temporary B for use in TSQR.
00243       // If it is provided, adjust dimensions as necessary.
00244       if (B == Teuchos::null)
00245   B = Teuchos::rcp (new serial_matrix_type (ncols, ncols));
00246       else
00247   B->shape (ncols, ncols);
00248       //
00249       // Compute rank-revealing decomposition (in this case, TSQR of X
00250       // followed by SVD of the R factor and appropriate updating of
00251       // the resulting Q factor) of X.  X is modified in place, and Q
00252       // contains the results.
00253       //
00254       // Compute TSQR and SVD of X.  Resulting orthogonal vectors go
00255       // into Q_, and coefficients (not necessarily upper triangular)
00256       // go into B.
00257       int rank;
00258       try {
00259   typedef typename tsqr_adaptor_type::factor_output_type 
00260     factor_output_type;
00261   factor_output_type factorOutput = tsqrAdaptor_->factor (X, *B);
00262   tsqrAdaptor_->explicitQ (X, factorOutput, *Q_);
00263   rank = tsqrAdaptor_->revealRank (*Q_, *B, relativeRankTolerance());
00264       } catch (std::exception& e) {
00265   throw OrthoError (e.what());
00266       }
00267       // Now we should copy Q_ back into X, but don't do it yet: if we
00268       // want to fill the last ncols-rank columns with random data, we
00269       // should do so in Q_, because it's fresh in the cache.
00270       if (false)
00271   {
00272     // If X did not have full (numerical rank), augment the last
00273     // ncols-rank columns of X with random data.
00274     const int ncolsToFill = ncols - rank;
00275     if (ncolsToFill > 0)
00276       {
00277         // ind: Indices of columns of X to fill with random data.
00278         std::vector< int > fillIndices (ncolsToFill);
00279         for (int j = 0; j < ncolsToFill; ++j)
00280     fillIndices[j] = j + rank;
00281 
00282         mv_ptr Q_null = MVT::CloneViewNonConst (*Q_, fillIndices);
00283         MVT::MvRandom (*Q_null);
00284         // Done with Q_null; tell Teuchos to deallocate it.
00285         Q_null = Teuchos::null;
00286       }
00287   }
00288       // MultiVecTraits (MVT) doesn't have a "deep copy from one MV
00289       // into an existing MV in place" method, but it does have two
00290       // methods which may be used to this effect:
00291       // 
00292       // 1. MvAddMv() (compute X := 1*Q_ + 0*X)
00293       //
00294       // 2. SetBlock() (Copy from A to mv by setting the index vector
00295       //    to [0, 1, ..., GetNumberVecs(mv)-1])
00296       //
00297       // MVT doesn't state the aliasing rules for MvAddMv(), but I'm
00298       // guessing it's OK for B and mv to alias one another (that is
00299       // the way AXPY works in the BLAS).
00300       MVT::MvAddMv (ScalarType(1), *Q_, ScalarType(0), X, X);
00301 
00302       // Don't deallocate B, even if the user provided Teuchos::null
00303       // as the B input (or the default B input took over).  The RCP's
00304       // destructor should work just fine.
00305       return rank;
00306     }
00307 
00308 
00309     int 
00310     projectAndNormalize (MV &X,
00311        const_prev_mvs_type Q,
00312        prev_coeffs_type C
00313        = Teuchos::tuple (serial_matrix_ptr (Teuchos::null)),
00314        serial_matrix_ptr B = Teuchos::null)
00315     {
00316       // if (_hasOp || _Op != Teuchos::null)
00317       //  throw NonNullOperatorError();
00318       // else if (MX != Teuchos::null)
00319       //  throw OrthoError("Sorry, TsqrOrthoManager::projectAndNormalizeMat() "
00320       //       "doesn\'t work with MX non-null yet");
00321 
00322       // Fetch dimensions of X and Q.
00323       int nrows_X, ncols_X, num_Q_blocks, ncols_Q_total;
00324       checkProjectionDims (nrows_X, ncols_X, num_Q_blocks, ncols_Q_total, X, Q);
00325 
00326       // Test for quick exit: any dimension of X is zero.
00327       if (nrows_X == 0 || ncols_X == 0)
00328   return 0;
00329 
00330       // If there are zero Q blocks or zero Q columns, just normalize!
00331       if (num_Q_blocks == 0 || ncols_Q_total == 0)
00332   return normalize (X, B);
00333       
00334       // If we don't have enough C, expanding it creates null references.
00335       // If we have too many, resizing just throws away the later ones.
00336       // If we have exactly as many as we have Q, this call has no effect.
00337       C.resize (num_Q_blocks);
00338       prev_coeffs_type newC (num_Q_blocks);
00339       for (int i = 0; i < num_Q_blocks; ++i) 
00340   {
00341     const int ncols_Q_i = MVT::GetNumberVecs (*Q[i]);
00342     // Create a new C[i] if necessary.
00343     if (C[i] == Teuchos::null)
00344       C[i] = Teuchos::rcp (new serial_matrix_type (ncols_Q_i, ncols_X));
00345     // Fill C[i] with zeros.
00346     C[i]->putScalar (ScalarType(0));
00347     // Create newC[i] as a clone of C[i].  (All that really
00348     // matters is that is has the same dimensions and is filled
00349     // with zeros; cloning C[i] accomplishes both in this case.)
00350     newC[i] = Teuchos::rcp (new serial_matrix_type (*C[i]));
00351   }
00352 
00353       // std::cerr << "Got past initialization of newC[0.." 
00354       //    << (num_Q_blocks-1) << "]" << std::endl;
00355 
00356       // Keep track of the column norms of X, both before and after
00357       // each orthogonalization pass.
00358       std::vector< MagnitudeType > normsBeforeFirstPass (ncols_X, MagnitudeType(0));
00359       std::vector< MagnitudeType > normsAfterFirstPass (ncols_X, MagnitudeType(0));
00360       std::vector< MagnitudeType > normsAfterSecondPass (ncols_X, MagnitudeType(0));
00361       MVT::MvNorm (X, normsBeforeFirstPass);
00362 
00363       // First BGS pass.  "Modified Gram-Schmidt" version of Block
00364       // Gram-Schmidt:
00365       //
00366       // \li \f$C^{\text{new}}_i := Q_i^* \cdot X\f$
00367       // \li \f$X := X - Q_i \cdot C^{\text{new}}_i\f$
00368       rawProject (X, Q, newC);
00369 
00370       // std::cerr << "Got past rawProject(X,Q,newC)" << std::endl;
00371 
00372       // Update the C matrices:
00373       //
00374       // \li \f$C_i := C_i + C^{\text{new}}_i\f$
00375       for (int i = 0; i < num_Q_blocks; ++i)
00376   *C[i] += *newC[i];
00377 
00378       // std::cerr << "Got past *C[i] += *newC[i]" << std::endl;
00379 
00380       // Normalize the matrix X.
00381       if (B == Teuchos::null)
00382   B = Teuchos::rcp (new serial_matrix_type (ncols_X, ncols_X));
00383       int rank = normalize (X, B);
00384 
00385       // std::cerr << "Got past normalize(X, B)" << std::endl;
00386       
00387       // Compute post-first-pass (pre-normalization) norms, using B.
00388       // normalize() doesn't guarantee in general that B is upper
00389       // triangular, so we compute norms using the entire column of B.
00390       Teuchos::BLAS< int, ScalarType > blas;
00391       for (int j = 0; j < ncols_X; ++j)
00392   {
00393     const ScalarType* const B_j = &(*B)(0,j);
00394     // Teuchos::BLAS returns a MagnitudeType result on
00395     // ScalarType inputs.
00396     normsAfterFirstPass[j] = blas.NRM2 (ncols_X, B_j, 1);
00397   }
00398       // Test whether any of the norms dropped below the
00399       // reorthogonalization threshold.
00400       bool reorthog = false;
00401       for (int j = 0; j < ncols_X; ++j)
00402   if (normsBeforeFirstPass[j] / normsAfterFirstPass[j] <= blockReorthogThreshold())
00403     {
00404       reorthog = true; 
00405       break;
00406     }
00407 
00408       // Perform another BGS pass if necessary.  "Twice is enough"
00409       // (Kahan's theorem) for a Krylov method, unless (using
00410       // Stewart's term) there is an "orthogonalization fault"
00411       // (indicated by reorthogFault).
00412       bool reorthogFault = false;
00413       // Indices of X at which there was an orthogonalization fault.
00414       std::vector< int > faultIndices (0);
00415       if (reorthog)
00416   {
00417     // Block Gram-Schmidt (again):
00418     //
00419     // \li \f$C^{\text{new}} = Q^* X\f$
00420     // \li \f$X := X - Q C^{\text{new}}\f$
00421     // \li \f$C := C + C^{\text{new}}\f$
00422     rawProject (X, Q, newC);
00423     for (int i = 0; i < num_Q_blocks; ++i)
00424       *C[i] += *newC[i];
00425 
00426     // Normalize the matrix X.
00427     serial_matrix_ptr B_new (new serial_matrix_type (ncols_X, ncols_X));
00428     rank = normalize (X, B_new);
00429     *B += *B_new;
00430 
00431     // Compute post-second-pass (pre-normalization) norms, using
00432     // B.  normalize() doesn't guarantee in general that B is
00433     // upper triangular, so we compute norms using the entire
00434     // column of B.
00435     for (int j = 0; j < ncols_X; ++j)
00436       {
00437         // FIXME Should we use B_new or B here?
00438         const ScalarType* const B_j = &(*B)(0,j);
00439         // Teuchos::BLAS returns a MagnitudeType result on
00440         // ScalarType inputs.
00441         normsAfterSecondPass[j] = blas.NRM2 (ncols_X, B_j, 1);
00442       }
00443     // Test whether any of the norms dropped below the
00444     // reorthogonalization threshold.  If so, it's an
00445     // orthogonalization fault, which requires expensive
00446     // recovery.
00447     reorthogFault = false;
00448     for (int j = 0; j < ncols_X; ++j)
00449       if (normsAfterSecondPass[j] / normsAfterFirstPass[j] <= blockReorthogThreshold())
00450         {
00451     reorthogFault = true; 
00452     faultIndices.push_back (j);
00453         }
00454   } // if (reorthog) // reorthogonalization pass
00455 
00456       if (reorthogFault)
00457   {
00458     if (throwOnReorthogFault_)
00459       {
00460         using std::endl;
00461         typedef std::vector<int>::size_type size_type;
00462         std::ostringstream os;
00463 
00464         os << "Orthogonalization fault at the following column(s) of X:" << endl;
00465         os << "Column\tNorm decrease factor" << endl;
00466         for (size_type k = 0; k < faultIndices.size(); ++k)
00467     {
00468       const int index = faultIndices[k];
00469       const MagnitudeType decreaseFactor = 
00470         normsAfterSecondPass[index] / normsAfterFirstPass[index];
00471       os << index << "\t" << decreaseFactor << endl;
00472     }
00473         throw TsqrOrthoFault (os.str());
00474       }
00475     else // Slowly reorthogonalize, one vector at a time, the offending vectors of X.
00476       {
00477         throw std::logic_error ("Not implemented yet");
00478 
00479         // for (int k = 0; k < faultIndices.size(); ++k)
00480         //  {
00481         //    // Get a view of the current column of X.
00482         //    std::vector<int> currentIndex (1, faultIndices[k]);
00483         //    Teuchos::RCP< MV > X_j = MVT::CloneViewNonConst (X, currentIndex);
00484 
00485         //    // Reorthogonalize X_j against all columns of each Q[i].
00486         //  }
00487       }
00488   }
00489       return rank;
00490     }
00491 
00496     MagnitudeType 
00497     orthonormError (const MV &X) const
00498     {
00499       const ScalarType ONE = SCT::one();
00500       const int ncols = MVT::GetNumberVecs(X);
00501       serial_matrix_type XTX (ncols, ncols);
00502       innerProd (X, X, XTX);
00503       for (int k = 0; k < ncols; ++k)
00504   XTX(k,k) -= ONE;
00505       return XTX.normFrobenius();
00506     }
00507 
00508     MagnitudeType 
00509     orthogError (const MV &X1, 
00510      const MV &X2) const
00511     {
00512       const int ncols_X1 = MVT::GetNumberVecs (X1);
00513       const int ncols_X2 = MVT::GetNumberVecs (X2);
00514       serial_matrix_type X1_T_X2 (ncols_X1, ncols_X2);
00515       innerProd (X1, X2, X1_T_X2);
00516       return X1_T_X2.normFrobenius();
00517     }
00518 
00522     MagnitudeType blockReorthogThreshold() const { return blockReorthogThreshold_; }
00525     MagnitudeType relativeRankTolerance() const { return relativeRankTolerance_; }
00526 
00527   private:
00530     Teuchos::ParameterList tsqrParams_;
00534     tsqr_adaptor_ptr tsqrAdaptor_;
00537     mv_ptr Q_;
00539     // Machine precision for ScalarType
00540     MagnitudeType eps_;
00543     MagnitudeType blockReorthogThreshold_;
00546     MagnitudeType relativeRankTolerance_;
00547     
00550     bool throwOnReorthogFault_;
00551 
00559     void
00560     lazyInit (const MV& X)
00561     {
00562       // The TSQR adaptor object requires a specific MV object for
00563       // initialization.  As long as subsequent MV objects use the
00564       // same communicator (e.g., the same Teuchos::Comm<int>), we
00565       // don't need to reinitialize the adaptor.
00566       //
00567       // FIXME (mfh 15 Jul 2010) If tsqrAdaptor_ has already been
00568       // initialized, check to make sure that X has the same
00569       // communicator as the multivector previously used to initialize
00570       // tsqrAdaptor_.
00571       if (tsqrAdaptor_ == Teuchos::null)
00572   tsqrAdaptor_ = Teuchos::rcp (new tsqr_adaptor_type (X, tsqrParams_));
00573 
00574       const int nrows = MVT::GetVecLength (X);
00575       const int ncols = MVT::GetNumberVecs (X);
00576 
00577       // Q_ is temporary workspace.  It must have the same dimensions
00578       // as X.  If not, we have to reallocate.  We also have to
00579       // allocate (not "re-") if we haven't allocated Q_ before.  (We
00580       // can't allocate Q_ until we have some X, so we need a
00581       // multivector as the "prototype.")
00582       if (Q_ == Teuchos::null || 
00583     nrows != MVT::GetVecLength(X) || 
00584     ncols != MVT::GetVecLength(X))
00585   // Allocate Q_ to have the same dimensions as X.  Contents of
00586   // X are not copied.
00587   Q_ = MVT::Clone (X, ncols);
00588     }
00589 
00590     void
00591     checkProjectionDims (int& nrows_X, 
00592        int& ncols_X, 
00593        int& num_Q_blocks,
00594        int& ncols_Q_total,
00595        const MV& X, 
00596        const_prev_mvs_type Q) const
00597     {
00598       // Test for quick exit
00599       num_Q_blocks = Q.length();
00600       nrows_X = MVT::GetVecLength (X);
00601       ncols_X = MVT::GetNumberVecs (X);
00602 
00603       //
00604       // Make sure that for each i, the dimensions of X and Q[i] are
00605       // compatible.
00606       //
00607       ncols_Q_total = 0; // total over all Q blocks
00608       for (int i = 0; i < num_Q_blocks; ++i)
00609   {
00610     const int nrows_Q = MVT::GetVecLength (*Q[i]);
00611     TEST_FOR_EXCEPTION( (nrows_Q != nrows_X), 
00612             std::invalid_argument,
00613             "Anasazi::TsqrOrthoManager::project(): "
00614             "Size of X not consistant with size of Q" );
00615     ncols_Q_total += MVT::GetNumberVecs (*Q[i]);
00616   }
00617     }
00618 
00619 
00622     void
00623     rawProject (MV& X, 
00624     const_prev_mvs_type Q,
00625     prev_coeffs_type C = Teuchos::tuple (serial_matrix_ptr (Teuchos::null))) const
00626     {
00627       int nrows_X, ncols_X, num_Q_blocks, ncols_Q_total;
00628       checkProjectionDims (nrows_X, ncols_X, num_Q_blocks, ncols_Q_total, X, Q);
00629       // Test for quick exit: any dimension of X is zero, or there are
00630       // zero Q blocks, or the total number of columns of the Q blocks
00631       // is zero.
00632       if (nrows_X == 0 || ncols_X == 0 || num_Q_blocks == 0 || ncols_Q_total == 0)
00633   return;
00634 
00635       // "Modified Gram-Schmidt" version of Block Gram-Schmidt.
00636       for (int i = 0; i < num_Q_blocks; ++i)
00637   {
00638     if (C[i] == Teuchos::null)
00639       {
00640         std::ostringstream os;
00641         os << "C[" << i << "] is null" << std::endl;
00642         throw std::logic_error (os.str());
00643       }
00644     innerProd (*Q[i], X, *C[i]);
00645     MVT::MvTimesMatAddMv (ScalarType(-1), *Q[i], *C[i], ScalarType(1), X);
00646   }
00647     }
00648   };
00649 
00650 
00656   template< class ScalarType, class MV >
00657   class TsqrOrthoManager : public OrthoManager< ScalarType, MV > {
00658   public:
00659     TsqrOrthoManager (const Teuchos::ParameterList& tsqrParams) :
00660       impl_ (tsqrParams)
00661     {}
00662 
00663     virtual ~TsqrOrthoManager() {}
00664 
00665     virtual void 
00666     innerProd (const MV &X, 
00667          const MV &Y, 
00668          Teuchos::SerialDenseMatrix<int,ScalarType>& Z) const
00669     {
00670       return impl_.innerProd (X, Y, Z);
00671     }
00672 
00673     virtual void 
00674     norm (const MV& X, 
00675     std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec) const
00676     {
00677       return impl_.norm (X, normvec);
00678     }
00679 
00680     virtual void 
00681     project (MV &X, 
00682        Teuchos::Array<Teuchos::RCP<const MV> > Q,
00683        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C 
00684        = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null))) const
00685     {
00686       return impl_.project (X, Q, C);
00687     }
00688 
00689     virtual int 
00690     normalize (MV &X, 
00691          Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null) const
00692     {
00693       return impl_.normalize (X, B);
00694     }
00695 
00696     virtual int 
00697     projectAndNormalize (MV &X, 
00698        Teuchos::Array<Teuchos::RCP<const MV> > Q,
00699        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C 
00700        = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00701        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null) const
00702     {
00703       return impl_.projectAndNormalize (X, Q, C, B);
00704     }
00705 
00706     virtual typename Teuchos::ScalarTraits< ScalarType >::magnitudeType 
00707     orthonormError (const MV &X) const
00708     {
00709       return impl_.orthonormError (X);
00710     }
00711 
00712     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00713     orthogError (const MV &X1, 
00714      const MV &X2) const 
00715     {
00716       return impl_.orthogError (X1, X2);
00717     }
00718 
00719   private:
00723     mutable TsqrOrthoManagerImpl< ScalarType, MV > impl_;
00724   };
00725 
00726 
00738   template< class ScalarType, class MV, class OP >
00739   class TsqrMatOrthoManager : public MatOrthoManager< ScalarType, MV, OP > {
00740   private:
00741     // This will be used to help C++ resolve getOp().  We can't call
00742     // getOp() directly, because C++ can't figure out that it belongs
00743     // to the base class MatOrthoManager.  (Remember that at this
00744     // point, we might not have specialized the specific base class
00745     // yet; it's just a template at the moment and not a "real
00746     // class.")
00747     typedef MatOrthoManager< ScalarType, MV, OP > base_type;
00748 
00749     typedef TsqrOrthoManagerImpl< ScalarType, MV > tsqr_type;
00750     typedef SVQBOrthoManager< ScalarType, MV, OP > svqb_type;
00751 
00752   public:
00753     typedef Teuchos::RCP< MV >       mv_ptr;
00754     typedef Teuchos::RCP< const MV > const_mv_ptr;
00755     typedef Teuchos::Array< const_mv_ptr >                       const_prev_mvs_type;
00756     typedef Teuchos::SerialDenseMatrix< int, ScalarType >        serial_matrix_type;
00757     typedef Teuchos::RCP< serial_matrix_type >                   serial_matrix_ptr;
00758     typedef Teuchos::Array< Teuchos::RCP< serial_matrix_type > > prev_coeffs_type;
00759     typedef typename Teuchos::ScalarTraits< ScalarType >::magnitudeType magnitude_type; 
00760 
00768     TsqrMatOrthoManager (const Teuchos::ParameterList& tsqrParams, 
00769        Teuchos::RCP< const OP > Op = Teuchos::null) :
00770       MatOrthoManager< ScalarType, MV, OP >(Op),
00771       tsqrParams_ (tsqrParams),
00772       pTsqr_ (Teuchos::null), // Lazy initialization
00773       pSvqb_ (Teuchos::null)  // Lazy initialization
00774     {}
00775 
00778     virtual ~TsqrMatOrthoManager() {}
00779 
00788     virtual void 
00789     setOp (Teuchos::RCP< const OP > Op) 
00790     {
00791       // We use this notation to help C++ resolve the name.
00792       // Otherwise, it won't know where to find setOp(), since this is
00793       // a member function of the base class which does not depend on
00794       // the template parameters.
00795       base_type::setOp (Op); // base class gets a copy of the Op too
00796       pSvqb_->setOp (Op);
00797     }
00799     virtual Teuchos::RCP< const OP > getOp () const { return base_type::getOp(); }
00800 
00801     virtual void 
00802     projectMat (MV &X, 
00803     const_prev_mvs_type Q,
00804     prev_coeffs_type C = Teuchos::tuple (serial_matrix_ptr (Teuchos::null)),
00805     mv_ptr MX = Teuchos::null,
00806     const_prev_mvs_type MQ = Teuchos::tuple (const_mv_ptr (Teuchos::null))) const
00807     {
00808       if (getOp() == Teuchos::null)
00809   {
00810     ensureTsqrInit ();
00811     pTsqr_->project (X, Q, C);
00812     // FIXME (mfh 20 Jul 2010) What about MX and MQ?
00813   }
00814       else
00815   {
00816     ensureSvqbInit ();
00817     pSvqb_->projectMat (X, Q, C, MX, MQ);
00818   }
00819     }
00820 
00821     virtual int 
00822     normalizeMat (MV &X, 
00823       serial_matrix_ptr B = Teuchos::null,
00824       mv_ptr MX = Teuchos::null) const
00825     {
00826       if (getOp() == Teuchos::null)
00827   {
00828     ensureTsqrInit ();
00829     return pTsqr_->normalize (X, B);
00830     // FIXME (mfh 20 Jul 2010) What about MX?
00831   }
00832       else
00833   {
00834     ensureSvqbInit ();
00835     return pSvqb_->normalizeMat (X, B, MX);
00836   }
00837     }
00838 
00839     virtual int 
00840     projectAndNormalizeMat (MV &X, 
00841           const_prev_mvs_type Q,
00842           prev_coeffs_type C = Teuchos::tuple (serial_matrix_ptr (Teuchos::null)),
00843           serial_matrix_ptr B = Teuchos::null,
00844           mv_ptr MX = Teuchos::null,
00845           const_prev_mvs_type MQ = Teuchos::tuple (const_mv_ptr (Teuchos::null))) const
00846     {
00847       if (getOp() == Teuchos::null)
00848   {
00849     ensureTsqrInit ();
00850     return pTsqr_->projectAndNormalize (X, Q, C, B); 
00851     // FIXME (mfh 20 Jul 2010) What about MX and MQ?
00852   }
00853       else
00854   {
00855     ensureSvqbInit ();
00856     return pSvqb_->projectAndNormalizeMat (X, Q, C, B, MX, MQ);
00857   }
00858     }
00859 
00860     virtual magnitude_type
00861     orthonormErrorMat (const MV &X, 
00862            const_mv_ptr MX = Teuchos::null) const
00863     {
00864       if (getOp() == Teuchos::null)
00865   {
00866     ensureTsqrInit ();
00867     return pTsqr_->orthonormError (X);
00868     // FIXME (mfh 20 Jul 2010) What about MX?
00869   }
00870       else
00871   {
00872     ensureSvqbInit ();
00873     return pSvqb_->orthonormErrorMat (X, MX);
00874   }
00875     }
00876 
00877     virtual magnitude_type
00878     orthogErrorMat (const MV &X, 
00879         const MV &Y,
00880         const_mv_ptr MX = Teuchos::null, 
00881         const_mv_ptr MY = Teuchos::null) const
00882     {
00883       if (getOp() == Teuchos::null)
00884   {
00885     ensureTsqrInit ();
00886     return pTsqr_->orthogError (X, Y);
00887     // FIXME (mfh 20 Jul 2010) What about MX and MY?
00888   }
00889       else
00890   {
00891     ensureSvqbInit ();
00892     return pSvqb_->orthogErrorMat (X, Y, MX, MY);
00893   }
00894     }
00895 
00896   private:
00897     void
00898     ensureTsqrInit () const
00899     {
00900       if (pTsqr_ == Teuchos::null)
00901   pTsqr_ = Teuchos::rcp (new tsqr_type (tsqrParams_));
00902     }
00903     void 
00904     ensureSvqbInit () const
00905     {
00906       if (pSvqb_ == Teuchos::null)
00907   pSvqb_ = Teuchos::rcp (new svqb_type (getOp()));
00908     }
00909 
00912     Teuchos::ParameterList tsqrParams_;
00916     mutable Teuchos::RCP< tsqr_type > pTsqr_;
00920     mutable Teuchos::RCP< svqb_type > pSvqb_;
00921   };
00922 
00923 } // namespace Anasazi
00924 
00925 #endif // __AnasaziTsqrOrthoManager_hpp
00926 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends