Belos Version of the Day
BelosTsqrOrthoManagerImpl.hpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                  Copyright 2004 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00045 
00046 #ifndef __BelosTsqrOrthoManagerImpl_hpp
00047 #define __BelosTsqrOrthoManagerImpl_hpp
00048 
00049 #include "BelosConfigDefs.hpp" // HAVE_BELOS_TSQR
00050 #include "BelosMultiVecTraits.hpp"
00051 #include "BelosOrthoManager.hpp" // OrthoError, etc.
00052 
00053 #include "Teuchos_LAPACK.hpp"
00054 #include "Teuchos_ParameterList.hpp"
00055 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00056 #  include "Teuchos_TimeMonitor.hpp"
00057 #endif // BELOS_TEUCHOS_TIME_MONITOR
00058 #include <algorithm>
00059 
00062 
00063 namespace Belos {
00064 
00067   class TsqrOrthoError : public OrthoError
00068   {
00069   public: 
00070     TsqrOrthoError (const std::string& what_arg) : 
00071       OrthoError (what_arg) {}
00072   };
00073 
00092   class TsqrOrthoFault : public OrthoError
00093   {
00094   public: 
00095     TsqrOrthoFault (const std::string& what_arg) : 
00096       OrthoError (what_arg) {}
00097   };
00098 
00127   template<class Scalar, class MV>
00128   class TsqrOrthoManagerImpl
00129   {
00130   public:
00131     typedef Scalar scalar_type;
00132     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00133     typedef MV multivector_type;
00136     typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00137     typedef Teuchos::RCP<mat_type> mat_ptr;
00138 
00139   private:
00140     typedef Teuchos::ScalarTraits<Scalar> SCT;
00141     typedef Teuchos::ScalarTraits<magnitude_type> SCTM;
00142     typedef MultiVecTraits<Scalar, MV> MVT;
00143     typedef typename MVT::tsqr_adaptor_type tsqr_adaptor_type;
00144     typedef Teuchos::RCP<tsqr_adaptor_type> tsqr_adaptor_ptr;
00145 
00146   public:
00164     static Teuchos::RCP<const Teuchos::ParameterList> getDefaultParameters ();
00165 
00181     static Teuchos::RCP<const Teuchos::ParameterList> getFastParameters ();
00182 
00201     TsqrOrthoManagerImpl (const Teuchos::RCP<const Teuchos::ParameterList>& params,
00202         const std::string& label);
00203 
00205     void setLabel (const std::string& label) { label_ = label; }
00206 
00208     const std::string& getLabel () { return label_; }
00209 
00218     void 
00219     innerProd (const MV& X, const MV& Y, mat_type& Z) const
00220     {
00221       MVT::MvTransMv (SCT::one(), X, Y, Z);
00222     }
00223 
00240     void
00241     norm (const MV& X, std::vector<magnitude_type>& normVec) const
00242     {
00243       const int numCols = MVT::GetNumberVecs (X);
00244       // std::vector<T>::size_type is unsigned; int is signed.  Mixed
00245       // unsigned/signed comparisons trigger compiler warnings.
00246       if (normVec.size() < static_cast<size_t>(numCols))
00247   normVec.resize (numCols); // Resize normvec if necessary.
00248       MVT::MvNorm (X, normVec);
00249     }
00250 
00260     void 
00261     project (MV& X, 
00262        Teuchos::Array<mat_ptr> C, 
00263        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q);
00264 
00278     int normalize (MV& X, mat_ptr B);
00279 
00298     int 
00299     normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B)
00300     {
00301 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00302       Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
00303       Teuchos::TimeMonitor timerMonitorNormalize(*timerNormalize_);
00304 #endif // BELOS_TEUCHOS_TIME_MONITOR
00305 
00306       if (MVT::GetNumberVecs(X) == 1)
00307   {
00308     using Teuchos::Range1D;
00309     using Teuchos::RCP;
00310     using Teuchos::rcp;
00311 
00312     // Make space for the normalization coefficient
00313     if (B.is_null())
00314       B = rcp (new mat_type (1, 1));
00315     else 
00316       B->shape (1, 1);
00317     // Normalize X in place (faster than TSQR for one column)
00318     const int rank = normalizeOne (X, B);
00319     // Copy results to first column of Q
00320     RCP<MV> Q_0 = MVT::CloneViewNonConst (Q, Range1D(0,0));
00321     MVT::Assign (X, *Q_0);
00322     return rank;
00323   }
00324       else
00325   // "true" means the output vectors go into Q, and the contents
00326   // of X are overwritten with invalid values.
00327   return normalizeImpl (X, Q, B, true);
00328     }
00329 
00340     int 
00341     projectAndNormalize (MV &X,
00342        Teuchos::Array<mat_ptr> C,
00343        mat_ptr B,
00344        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
00345     {
00346       // "false" means we work on X in place.  The second argument is
00347       // not read or written in that case.
00348       return projectAndNormalizeImpl (X, X, false, C, B, Q);
00349     }
00350 
00370     int 
00371     projectAndNormalizeOutOfPlace (MV& X_in, 
00372            MV& X_out,
00373            Teuchos::Array<mat_ptr> C,
00374            mat_ptr B,
00375            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
00376     {
00377       // "true" means we work on X_in out of place, writing the
00378       // results into X_out.
00379       return projectAndNormalizeImpl (X_in, X_out, true, C, B, Q);
00380     }
00381     
00386     magnitude_type 
00387     orthonormError (const MV &X) const
00388     {
00389       const Scalar ONE = SCT::one();
00390       const int ncols = MVT::GetNumberVecs(X);
00391       mat_type XTX (ncols, ncols);
00392       innerProd (X, X, XTX);
00393       for (int k = 0; k < ncols; ++k)
00394   XTX(k,k) -= ONE;
00395       return XTX.normFrobenius();
00396     }
00397 
00399     magnitude_type 
00400     orthogError (const MV &X1, 
00401      const MV &X2) const
00402     {
00403       const int ncols_X1 = MVT::GetNumberVecs (X1);
00404       const int ncols_X2 = MVT::GetNumberVecs (X2);
00405       mat_type X1_T_X2 (ncols_X1, ncols_X2);
00406       innerProd (X1, X2, X1_T_X2);
00407       return X1_T_X2.normFrobenius();
00408     }
00409 
00413     magnitude_type blockReorthogThreshold() const { return blockReorthogThreshold_; }
00414 
00417     magnitude_type relativeRankTolerance() const { return relativeRankTolerance_; }
00418 
00419   private:
00421     Teuchos::RCP<const Teuchos::ParameterList> params_;
00423     std::string label_;
00425     tsqr_adaptor_ptr tsqrAdaptor_;
00432     Teuchos::RCP<MV> Q_;
00433 
00435     magnitude_type eps_;
00436 
00439     bool randomizeNullSpace_;
00441     bool reorthogonalizeBlocks_;
00444     bool throwOnReorthogFault_;
00446     magnitude_type blockReorthogThreshold_;
00448     magnitude_type relativeRankTolerance_;
00449 
00450 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00451 
00452     Teuchos::RCP<Teuchos::Time> timerOrtho_;
00454     Teuchos::RCP<Teuchos::Time> timerProject_;
00456     Teuchos::RCP<Teuchos::Time> timerNormalize_;
00457 
00466     static Teuchos::RCP<Teuchos::Time>
00467     makeTimer (const std::string& prefix, 
00468          const std::string& timerName)
00469     {
00470       const std::string timerLabel = 
00471   prefix.empty() ? timerName : (prefix + ": " + timerName);
00472       return Teuchos::TimeMonitor::getNewTimer (timerLabel);
00473     }
00474 #endif // BELOS_TEUCHOS_TIME_MONITOR
00475 
00477     void
00478     raiseReorthogFault (const std::vector<magnitude_type>& normsAfterFirstPass,
00479       const std::vector<magnitude_type>& normsAfterSecondPass,
00480       const std::vector<int>& faultIndices);
00481 
00487     bool 
00488     getBoolParamWithDefault (const Teuchos::RCP<const Teuchos::ParameterList>& params, 
00489            const std::string& key)
00490     {
00491       using Teuchos::Exceptions::InvalidParameterName;
00492       using Teuchos::Exceptions::InvalidParameterType;
00493       // Note: We don't try to catch an exception when reading from
00494       // the default parameters, since they key had better be there.
00495       // (I guess we could catch and rethrow std::logic_error...)
00496       Teuchos::RCP<const Teuchos::ParameterList> defaultParams = 
00497   getDefaultParameters();
00498       if (params.is_null())
00499   return defaultParams->get< bool >(key);
00500       else
00501   {
00502     try {
00503       // No validation is necessary or even sensible; validating
00504       // a boolean value would mean it's not really an option!
00505       return params->get< bool >(key);
00506     } catch (InvalidParameterName&) {
00507       return defaultParams->get< bool >(key);
00508     } catch (InvalidParameterType&) {
00509       std::ostringstream os;
00510       os << "The value of parameter \"" << key << "\" to Tsqr(Mat)"
00511         "OrthoManager(Impl) must be a boolean truth value.  Here "
00512         "is the documentation for that parameter:" << std::endl
00513          << defaultParams->getEntry(key).docString();
00514       throw std::invalid_argument (os.str());
00515     }
00516   }
00517     }
00518 
00524     magnitude_type
00525     getMagParamWithDefault (const Teuchos::RCP<const Teuchos::ParameterList>& params, 
00526           const std::string& key)
00527     {
00528       using Teuchos::Exceptions::InvalidParameterName;
00529       using Teuchos::Exceptions::InvalidParameterType;
00530       // Note: We don't try to catch an exception when reading from
00531       // the default parameters, since they key had better be there.
00532       // (I guess we could catch and rethrow std::logic_error...)
00533       Teuchos::RCP<const Teuchos::ParameterList> defaultParams = 
00534   getDefaultParameters();
00535       if (params.is_null())
00536   return defaultParams->get< magnitude_type >(key);
00537       else
00538   {
00539     try {
00540       const magnitude_type value = params->get< magnitude_type >(key);
00541       if (value >= magnitude_type(0)) // Validate
00542         return value;
00543       else
00544         {
00545     std::ostringstream os;
00546     os << "You specified " << key << " = " << value
00547        << ", but that parameter must be a nonnegative real "
00548        << "floating-point value.  Here is the documentation " 
00549        << "for that parameter:" << std::endl
00550        << defaultParams->getEntry(key).docString();
00551     throw std::invalid_argument (os.str());
00552         }
00553     } catch (InvalidParameterName&) { 
00554       return defaultParams->get< magnitude_type >(key);
00555     } catch (InvalidParameterType&) {
00556       std::ostringstream os;
00557       os << "The value of parameter \"" << key << "\" to Tsqr(Mat)OrthoMa"
00558         "nager(Impl) must be a nonnegative real floating-point value.  "
00559         "Here is the documentation for that parameter:" << std::endl
00560          << defaultParams->getEntry(key).docString();
00561       throw std::invalid_argument (os.str());
00562     }
00563   }
00564     }
00565 
00567     void
00568     readParams (const Teuchos::RCP<const Teuchos::ParameterList>& params)
00569     {
00570       randomizeNullSpace_ = 
00571   getBoolParamWithDefault (params, "randomizeNullSpace");
00572       reorthogonalizeBlocks_ = 
00573   getBoolParamWithDefault (params, "reorthogonalizeBlocks");
00574       throwOnReorthogFault_ = 
00575   getBoolParamWithDefault (params, "throwOnReorthogFault");
00576       blockReorthogThreshold_ = 
00577   getMagParamWithDefault (params, "blockReorthogThreshold");
00578       relativeRankTolerance_ = 
00579   getMagParamWithDefault (params, "relativeRankTolerance");
00580     }
00581 
00589     void
00590     lazyInit (const MV& X)
00591     {
00592       using Teuchos::Exceptions::InvalidParameter;
00593       using Teuchos::ParameterList;
00594       using Teuchos::parameterList;
00595       using Teuchos::rcpFromRef;
00596       using Teuchos::RCP;
00597       using Teuchos::rcp;
00598 
00599       // TSQR parameters.  They are stored as an RCP<const
00600       // ParameterList>, so we don't have to worry about the sublist
00601       // going out of scope.
00602       RCP<const ParameterList> tsqrParams;
00603       bool gotTsqrParams = false;
00604       try {
00605   tsqrParams = params_->get<RCP<const ParameterList> >("TsqrImpl");
00606   gotTsqrParams = true;
00607       } catch (InvalidParameter&) {
00608   // We'll try to fetch "TsqrImpl" as a sublist instead.
00609       }
00610       if (! gotTsqrParams)
00611   {
00612     try {
00613       // We don't need to make a deep copy, since the
00614       // tsqrAdaptor_ which receives this sublist will fall out
00615       // of scope at the same time as params_, and params_ will
00616       // protect the sublist from falling out of scope.
00617       tsqrParams = rcpFromRef (params_->sublist ("TsqrImpl"));
00618       gotTsqrParams = true;
00619     } catch (InvalidParameter&) {
00620       // We'll try to fetch "TsqrImpl" from the default parameter
00621       // list instead.
00622     }
00623   }
00624       if (! gotTsqrParams)
00625   {      
00626     RCP<const ParameterList> defaultParams = getDefaultParameters();
00627     try {
00628       tsqrParams = defaultParams->get<RCP<const ParameterList> >("TsqrImpl");
00629       gotTsqrParams = true;
00630     } catch (InvalidParameter&) {
00631       // We'll try to fetch "TsqrImpl" as a sublist instead.
00632     }
00633   }
00634       if (! gotTsqrParams) 
00635   {
00636     RCP<const ParameterList> defaultParams = getDefaultParameters();
00637     try {
00638       // We don't need to make a deep copy, since the
00639       // tsqrAdaptor_ which receives this sublist will fall out
00640       // of scope at the same time as params_, and params_ will
00641       // protect the sublist from falling out of scope.
00642       tsqrParams = rcpFromRef (defaultParams->sublist ("TsqrImpl"));
00643       gotTsqrParams = true;
00644     } catch (InvalidParameter&) {
00645       // We've tried to get TsqrImpl in all the ways we could.
00646       // If we haven't gotten it yet, or if we only got null,
00647       // we'll throw an exception below.
00648     }
00649   }
00650       TEST_FOR_EXCEPTION(!gotTsqrParams || tsqrParams.is_null(), 
00651        std::logic_error,
00652        "Belos::TsqrOrthoManagerImpl::lazyInit: Failed to find"
00653        " \"TsqrImpl\" parameter in either the input parameter"
00654        " list or in the default parameter list.  Please "
00655        "report this bug to the Belos developers.");
00656       // The TSQR adaptor object requires a specific MV object for
00657       // initialization.  As long as subsequent MV objects use the
00658       // same communicator (e.g., the same Teuchos::Comm<int>), we
00659       // don't need to reinitialize the adaptor.
00660       //
00661       // NOTE (mfh 15 Jul 2010) If tsqrAdaptor_ has already been
00662       // initialized, we really should check to make sure that X has
00663       // the same communicator as that of the the multivector with
00664       // which tsqrAdaptor_ was previously initialized.
00665       if (tsqrAdaptor_.is_null())
00666   tsqrAdaptor_ = rcp (new tsqr_adaptor_type (X, tsqrParams));
00667 
00668       // NOTE: We don't try to allocate the temporary workspace Q_
00669       // here.  This is done on demand in normalize().
00670     }
00671 
00681     void
00682     checkProjectionDims (int& ncols_X, 
00683        int& num_Q_blocks,
00684        int& ncols_Q_total,
00685        const MV& X, 
00686        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00687 
00698     void
00699     allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C, 
00700             Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
00701             const MV& X,
00702             const bool attemptToRecycle = true) const;
00703 
00712     int 
00713     projectAndNormalizeImpl (MV& X_in, 
00714            MV& X_out,
00715            const bool outOfPlace,
00716            Teuchos::Array<mat_ptr> C,
00717            mat_ptr B,
00718            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q);
00719 
00725     void
00726     rawProject (MV& X, 
00727     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
00728     Teuchos::ArrayView<mat_ptr> C) const;
00729 
00731     void
00732     rawProject (MV& X, 
00733     const Teuchos::RCP<const MV>& Q, 
00734     const mat_ptr& C) const;
00735 
00762     int rawNormalize (MV& X, MV& Q, mat_type& B);
00763 
00781     int normalizeOne (MV& X, mat_ptr B) const;
00782 
00810     int normalizeImpl (MV& X, MV& Q, mat_ptr B, const bool outOfPlace);
00811   };
00812 
00813 
00814   template<class Scalar, class MV>
00815   TsqrOrthoManagerImpl<Scalar, MV>::
00816   TsqrOrthoManagerImpl (const Teuchos::RCP<const Teuchos::ParameterList>& params,
00817       const std::string& label) :
00818     params_ (params.is_null() ? 
00819        getDefaultParameters() : 
00820        Teuchos::rcp_const_cast<const Teuchos::ParameterList> (Teuchos::parameterList (*params))),
00821     label_ (label),
00822     tsqrAdaptor_ (Teuchos::null),   // Initialized on demand
00823     Q_ (Teuchos::null),             // Scratch space for normalize()
00824     eps_ (SCTM::eps()),             // Machine precision
00825     randomizeNullSpace_ (true),     // Set later by readParams()
00826     reorthogonalizeBlocks_ (true),  // Set later by readParams()
00827     throwOnReorthogFault_ (false),  // Set later by readParams()
00828     blockReorthogThreshold_ (0),    // Set later by readParams()
00829     relativeRankTolerance_ (0)      // Set later by readParams()
00830   {
00831 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00832     timerOrtho_ = makeTimer (label, "All orthogonalization");
00833     timerProject_ = makeTimer (label, "Projection");
00834     timerNormalize_ = makeTimer (label, "Normalization");
00835 #endif // BELOS_TEUCHOS_TIME_MONITOR
00836 
00837     // Extract values for the parameters from the given parameter
00838     // list.  Use default values if none are provided.  The "TSQR"
00839     // sublist gets passed along to the underlying TSQR
00840     // implementation.
00841     readParams (params_);
00842   }
00843 
00844   template<class Scalar, class MV>
00845   void
00846   TsqrOrthoManagerImpl<Scalar, MV>::project (MV& X, 
00847                Teuchos::Array<mat_ptr> C,
00848                Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
00849   {
00850 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00851     // "Projection" only happens in rawProject(), so we only time
00852     // projection inside rawProject().
00853     //
00854     // If project() is called from projectAndNormalize(), the
00855     // TimeMonitor won't start timerOrtho_, because it is already
00856     // running in projectAndNormalize().
00857     Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
00858 #endif // BELOS_TEUCHOS_TIME_MONITOR
00859 
00860     // Internal data used by this method require a specific MV object
00861     // for initialization (e.g., to get a Map / communicator, and to
00862     // initialize scratch space).  Thus, we delay (hence "lazy")
00863     // initialization until we get an X.
00864     lazyInit (X);
00865 
00866     int ncols_X, num_Q_blocks, ncols_Q_total;
00867     checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X, Q);
00868     // Test for quick exit: any dimension of X is zero, or there are
00869     // zero Q blocks, or the total number of columns of the Q blocks
00870     // is zero.
00871     if (ncols_X == 0 || num_Q_blocks == 0 || ncols_Q_total == 0)
00872       return;
00873 
00874     // Make space for first-pass projection coefficients
00875     allocateProjectionCoefficients (C, Q, X, true);
00876 
00877     // We only use columnNormsBefore and compute pre-projection column
00878     // norms if doing block reorthogonalization.
00879     std::vector<magnitude_type> columnNormsBefore (ncols_X, magnitude_type(0));
00880     if (reorthogonalizeBlocks_)
00881       MVT::MvNorm (X, columnNormsBefore);
00882 
00883     // Project (first block orthogonalization step): 
00884     // C := Q^* X, X := X - Q C.
00885     rawProject (X, Q, C); 
00886 
00887     // If we are doing block reorthogonalization, reorthogonalize X if
00888     // necessary.
00889     if (reorthogonalizeBlocks_)
00890       {
00891   std::vector<magnitude_type> columnNormsAfter (ncols_X, magnitude_type(0));
00892   MVT::MvNorm (X, columnNormsAfter);
00893   
00894   // Relative block reorthogonalization threshold
00895   const magnitude_type relThres = blockReorthogThreshold();
00896   // Reorthogonalize X if any of its column norms decreased by a
00897   // factor more than the block reorthogonalization threshold.
00898   // Don't bother trying to subset the columns; that will make
00899   // the columns noncontiguous and thus hinder BLAS 3
00900   // optimizations.
00901   bool reorthogonalize = false;
00902   for (int j = 0; j < ncols_X; ++j)
00903     if (columnNormsAfter[j] < relThres * columnNormsBefore[j])
00904       {
00905         reorthogonalize = true;
00906         break;
00907       }
00908   if (reorthogonalize)
00909     {
00910       // Second-pass projection coefficients
00911       Teuchos::Array<mat_ptr> C2;
00912       allocateProjectionCoefficients (C2, Q, X, false);
00913 
00914       // Perform the second projection pass:
00915       // C2 = Q' X, X = X - Q*C2
00916       rawProject (X, Q, C2);
00917       // Update the projection coefficients
00918       for (int k = 0; k < num_Q_blocks; ++k)
00919         *C[k] += *C2[k];
00920     }
00921       }
00922   }
00923 
00924 
00925   template<class Scalar, class MV>
00926   int 
00927   TsqrOrthoManagerImpl<Scalar, MV>::normalize (MV& X, mat_ptr B)
00928   {
00929     using Teuchos::Range1D;
00930     using Teuchos::RCP;
00931 
00932 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00933     Teuchos::TimeMonitor timerMonitorNormalize(*timerNormalize_);
00934     // If normalize() is called internally -- i.e., called from
00935     // projectAndNormalize() -- the timer will not be started or 
00936     // stopped, because it is already running.  TimeMonitor handles
00937     // recursive invocation by doing nothing.
00938     Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
00939 #endif // BELOS_TEUCHOS_TIME_MONITOR
00940 
00941     // Internal data used by this method require a specific MV object
00942     // for initialization (in particular, to get a Map / communicator
00943     // object, and to initialize scratch space).  Thus, we delay
00944     // (hence "lazy") initialization until we get an X.
00945     lazyInit (X);
00946 
00947     // MVT returns int for this, even though the "local ordinal
00948     // type" of the MV may be some other type (for example,
00949     // Tpetra::MultiVector<double, int32_t, int64_t, ...>).
00950     const int numCols = MVT::GetNumberVecs (X);
00951 
00952     // This special case (for X having only one column) makes
00953     // TsqrOrthoManagerImpl equivalent to Modified Gram-Schmidt
00954     // orthogonalization with conditional full reorthogonalization,
00955     // if all multivector inputs have only one column.  It also
00956     // avoids allocating Q_ scratch space and copying data when we
00957     // don't need to invoke TSQR at all.
00958     if (numCols == 1)
00959       return normalizeOne (X, B);
00960 
00961     // We use Q_ as scratch space for the normalization, since TSQR
00962     // requires a scratch multivector (it can't factor in place).  Q_
00963     // should come from a vector space compatible with X's vector
00964     // space, and Q_ should have at least as many columns as X.
00965     // Otherwise, we have to reallocate.  We also have to allocate
00966     // (not "re-") Q_ if we haven't allocated it before.  (We can't
00967     // allocate Q_ until we have some X, so we need a multivector as
00968     // the "prototype.")
00969     //
00970     // NOTE (mfh 11 Jan 2011) We only increase the number of columsn
00971     // in Q_, never decrease.  This is OK for typical uses of TSQR,
00972     // but you might prefer different behavior in some cases.
00973     //
00974     // NOTE (mfh 10 Mar 2011) We should only reuse the scratch space
00975     // Q_ if X and Q_ have compatible data distributions.  However,
00976     // Belos' current MultiVecTraits interface does not let us check
00977     // for this.  Thus, we can only check whether X and Q_ have the
00978     // same number of rows.  This will behave correctly for the common
00979     // case in Belos that all multivectors with the same number of
00980     // rows have the same data distribution.
00981     //
00982     // The specific MV implementation may do more checks than this on
00983     // its own and throw an exception if X and Q_ are not compatible,
00984     // but it may not.  If you find that recycling the Q_ space causes
00985     // troubles, you may consider modifying the code below to
00986     // reallocate Q_ for every X that comes in.  
00987     if (Q_.is_null() || 
00988   MVT::GetVecLength(*Q_) != MVT::GetVecLength(X) ||
00989   numCols > MVT::GetNumberVecs (*Q_))
00990       Q_ = MVT::Clone (X, numCols);
00991 
00992     // normalizeImpl() wants the second MV argument to have the same
00993     // number of columns as X.  To ensure this, we pass it a view of
00994     // Q_ if Q_ has more columns than X.  (This is possible if we
00995     // previously called normalize() with a different multivector,
00996     // since we never reallocate Q_ if it has more columns than
00997     // necessary.)
00998     if (MVT::GetNumberVecs(*Q_) == numCols)
00999       return normalizeImpl (X, *Q_, B, false);
01000     else
01001       {
01002   RCP<MV> Q_view = MVT::CloneViewNonConst (*Q_, Range1D(0, numCols-1));
01003   return normalizeImpl (X, *Q_view, B, false);
01004       }
01005   }
01006 
01007   template<class Scalar, class MV>
01008   void
01009   TsqrOrthoManagerImpl<Scalar, MV>::
01010   allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C, 
01011           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
01012           const MV& X,
01013           const bool attemptToRecycle) const
01014   {
01015     const int num_Q_blocks = Q.size();
01016     const int ncols_X = MVT::GetNumberVecs (X);
01017     C.resize (num_Q_blocks);
01018     if (attemptToRecycle)
01019       {
01020   for (int i = 0; i < num_Q_blocks; ++i) 
01021     {
01022       const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
01023       // Create a new C[i] if necessary, otherwise resize if
01024       // necessary, otherwise fill with zeros.
01025       if (C[i].is_null())
01026         C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
01027       else 
01028         {
01029     mat_type& Ci = *C[i];
01030     if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X)
01031       Ci.shape (ncols_Qi, ncols_X);
01032     else
01033       Ci.putScalar (SCT::zero());
01034         }
01035     }
01036       }
01037     else
01038       {
01039   for (int i = 0; i < num_Q_blocks; ++i) 
01040     {
01041       const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
01042       C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
01043     }
01044       }
01045   }
01046 
01047 
01048 
01049   template<class Scalar, class MV>
01050   int 
01051   TsqrOrthoManagerImpl<Scalar, MV>::
01052   projectAndNormalizeImpl (MV& X_in, 
01053          MV& X_out, // Only written if outOfPlace==false.
01054          const bool outOfPlace,
01055          Teuchos::Array<mat_ptr> C,
01056          mat_ptr B,
01057          Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
01058   {
01059     using Teuchos::Range1D;
01060     using Teuchos::RCP;
01061     using Teuchos::rcp;
01062 
01063 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01064     // Projection is only timed in rawProject(), and normalization is
01065     // only timed in normalize() and normalizeOutOfPlace().
01066     Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
01067 #endif // BELOS_TEUCHOS_TIME_MONITOR
01068 
01069     if (outOfPlace)
01070       { // Make sure that X_out has at least as many columns as X_in.
01071   TEST_FOR_EXCEPTION(MVT::GetNumberVecs(X_out) < MVT::GetNumberVecs(X_in),
01072          std::invalid_argument, 
01073          "Belos::TsqrOrthoManagerImpl::"
01074          "projectAndNormalizeOutOfPlace(...):"
01075          "X_out has " << MVT::GetNumberVecs(X_out) 
01076          << " columns, but X_in has "
01077          << MVT::GetNumberVecs(X_in) << " columns.");
01078       }
01079     // Fetch dimensions of X_in and Q, and allocate space for first-
01080     // and second-pass projection coefficients (C resp. C2).
01081     int ncols_X, num_Q_blocks, ncols_Q_total;
01082     checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X_in, Q);
01083 
01084     // Test for quick exit: any dimension of X is zero.
01085     if (ncols_X == 0)
01086       return 0;
01087     // If there are zero Q blocks or zero Q columns, just normalize!
01088     if (num_Q_blocks == 0 || ncols_Q_total == 0)
01089       {
01090   if (outOfPlace)
01091     return normalizeOutOfPlace (X_in, X_out, B);
01092   else
01093     return normalize (X_in, B);
01094       }
01095 
01096     // The typical case is that the entries of C have been allocated
01097     // before, so we attempt to recycle the allocations.  The call
01098     // below will be correct regardless.
01099     allocateProjectionCoefficients (C, Q, X_in, true);
01100 
01101     // If we are doing block reorthogonalization, then compute the
01102     // column norms of X before projecting for the first time.  This
01103     // will help us decide whether we need to reorthogonalize X.
01104     std::vector<magnitude_type> normsBeforeFirstPass (ncols_X, SCTM::zero());
01105     if (reorthogonalizeBlocks_)
01106       MVT::MvNorm (X_in, normsBeforeFirstPass);
01107 
01108     // First (Modified) Block Gram-Schmidt pass, in place in X_in
01109     rawProject (X_in, Q, C);
01110 
01111     // Before normalizing, be sure to allocate a new B matrix here.
01112     // If we let the normalize() routine allocate, that storage will
01113     // go away at the end of normalize().  (This is because it passes
01114     // the RCP by value, not by reference.)
01115     if (B.is_null())
01116       B = rcp (new mat_type (ncols_X, ncols_X));
01117 
01118     // Rank of X(_in) after first projection pass.  If outOfPlace,
01119     // this overwrites X_in with invalid values, and the results go in
01120     // X_out.  Otherwise, it's done in place in X_in.
01121     const int firstPassRank = outOfPlace ? 
01122       normalizeOutOfPlace (X_in, X_out, B) : 
01123       normalize (X_in, B);
01124     int rank = firstPassRank; // Current rank of X
01125 
01126     // If X was not full rank after projection and randomizeNullSpace_
01127     // is true, then normalize(OutOfPlace)() replaced the null space
01128     // basis of X with random vectors, and orthogonalized them against
01129     // the column space basis of X.  However, we still need to
01130     // orthogonalize the random vectors against the Q[i], after which
01131     // we will need to renormalize them.
01132     //
01133     // If outOfPlace, then we need to work in X_out (where
01134     // normalizeOutOfPlace() wrote the normalized vectors).
01135     // Otherwise, we need to work in X_in.
01136     //
01137     // Note: We don't need to keep the new projection coefficients,
01138     // since they are multiplied by the "small" part of B
01139     // corresponding to the null space of the original X.
01140     if (firstPassRank < ncols_X && randomizeNullSpace_)
01141       {
01142   const int numNullSpaceCols = ncols_X - firstPassRank;
01143   const Range1D nullSpaceIndices (firstPassRank, ncols_X - 1);
01144 
01145   // Space for projection coefficients (will be thrown away)
01146   Teuchos::Array<mat_ptr> C_null (num_Q_blocks);
01147   for (int k = 0; k < num_Q_blocks; ++k)
01148     {
01149       const int numColsQk = MVT::GetNumberVecs(*Q[k]);
01150       C_null[k] = rcp (new mat_type (numColsQk, numNullSpaceCols));
01151     }
01152   // Space for normalization coefficients (will be thrown away).
01153   RCP<mat_type> B_null (new mat_type (numNullSpaceCols, numNullSpaceCols));
01154 
01155   int randomVectorsRank;
01156   if (outOfPlace)
01157     {
01158       // View of the null space basis columns of X.
01159       // normalizeOutOfPlace() wrote them into X_out.
01160       RCP<MV> X_out_null = MVT::CloneViewNonConst (X_out, nullSpaceIndices);
01161       // Use X_in for scratch space.  Copy X_out_null into the
01162       // last few columns of X_in (X_in_null) and do projections
01163       // in there.  (This saves us a copy wen we renormalize
01164       // (out of place) back into X_out.)
01165       RCP<MV> X_in_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
01166       MVT::Assign (*X_out_null, *X_in_null);
01167       // Project the new random vectors against the Q blocks, and
01168       // renormalize the result into X_out_null.  
01169       rawProject (*X_in_null, Q, C_null);
01170       randomVectorsRank = normalizeOutOfPlace (*X_in_null, *X_out_null, B_null);
01171     }
01172   else
01173     {
01174       // View of the null space columns of X.  
01175       // They live in X_in.
01176       RCP<MV> X_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
01177       // Project the new random vectors against the Q blocks,
01178       // and renormalize the result (in place).
01179       rawProject (*X_null, Q, C_null);
01180       randomVectorsRank = normalize (*X_null, B_null);
01181     }
01182   // While unusual, it is still possible for the random data not
01183   // to be full rank after projection and normalization.  In
01184   // that case, we could try another set of random data and
01185   // recurse as necessary, but instead for now we just raise an
01186   // exception.
01187   TEST_FOR_EXCEPTION(randomVectorsRank != numNullSpaceCols, 
01188          TsqrOrthoError, 
01189          "Belos::TsqrOrthoManagerImpl::projectAndNormalize"
01190          "OutOfPlace(): After projecting and normalizing the "
01191          "random vectors (used to replace the null space "
01192          "basis vectors from normalizing X), they have rank "
01193          << randomVectorsRank << ", but should have full "
01194          "rank " << numNullSpaceCols << ".");
01195       }
01196 
01197     // Whether or not X_in was full rank after projection, we still
01198     // might want to reorthogonalize against Q.
01199     if (reorthogonalizeBlocks_)
01200       {
01201   // We are only interested in the column space basis of X
01202   // resp. X_out.
01203   std::vector<magnitude_type> normsAfterFirstPass (firstPassRank, magnitude_type(0));
01204   std::vector<magnitude_type> normsAfterSecondPass (firstPassRank, magnitude_type(0));
01205 
01206   // Compute post-first-pass (pre-normalization) norms.  We
01207   // could have done that using MVT::MvNorm() on X_in after
01208   // projecting, but before the first normalization.  However,
01209   // that operation may be expensive.  It is also unnecessary:
01210   // after calling normalize(), the 2-norm of B(:,j) is the
01211   // 2-norm of X_in(:,j) before normalization, in exact
01212   // arithmetic.
01213   //
01214   // NOTE (mfh 06 Nov 2010) This is one way that combining
01215   // projection and normalization into a single kernel --
01216   // projectAndNormalize() -- pays off.  In project(), we have
01217   // to compute column norms of X before and after projection.
01218   // Here, we get them for free from the normalization
01219   // coefficients.
01220   Teuchos::BLAS<int, Scalar> blas;
01221   for (int j = 0; j < firstPassRank; ++j)
01222     {
01223       const Scalar* const B_j = &(*B)(0,j);
01224       // Teuchos::BLAS returns a magnitude_type result on
01225       // Scalar inputs.
01226       normsAfterFirstPass[j] = blas.NRM2 (firstPassRank, B_j, 1);
01227     }
01228   // Test whether any of the norms dropped below the
01229   // reorthogonalization threshold.
01230   bool reorthogonalize = false;
01231   for (int j = 0; j < firstPassRank; ++j)
01232     // If any column's norm decreased too much, mark this block
01233     // for reorthogonalization.  Note that this test will _not_
01234     // activate reorthogonalization if a column's norm before
01235     // the first project-and-normalize step was zero.  It _will_
01236     // activate reorthogonalization if the column's norm before
01237     // was not zero, but is zero now.
01238     if (normsAfterFirstPass[j] < blockReorthogThreshold() * normsBeforeFirstPass[j])
01239       {
01240         reorthogonalize = true; 
01241         break;
01242       }
01243 
01244   // Perform another Block Gram-Schmidt pass if necessary.
01245   // "Twice is enough" (Kahan's theorem) for a Krylov method,
01246   // unless (using Stewart's term) there is an
01247   // "orthogonalization fault" (indicated by reorthogFault).
01248   //
01249   // NOTE (mfh 07 Nov 2010) For now, we include the entire block
01250   // of X, including possible random data (that was already
01251   // projected and normalized above).  It might make more sense
01252   // just to process the first firstPassRank columns of X.
01253   // However, the resulting reorthogonalization should still be
01254   // correct regardless.
01255   bool reorthogFault = false;
01256   // Indices of X at which there was an orthogonalization fault.
01257   std::vector<int> faultIndices;
01258   if (reorthogonalize)
01259     {
01260       using Teuchos::Copy;
01261       using Teuchos::NO_TRANS;
01262 
01263       // If we're using out-of-place normalization, copy X_out
01264       // (results of first project and normalize pass) back into
01265       // X_in, for the second project and normalize pass.
01266       if (outOfPlace)
01267         MVT::Assign (X_out, X_in);
01268 
01269       // C2 is only used internally, so we know that we are
01270       // allocating fresh and not recycling allocations.
01271       Teuchos::Array<mat_ptr> C2;
01272       allocateProjectionCoefficients (C2, Q, X_in, false);
01273 
01274       // Block Gram-Schmidt (again).  Delay updating the block
01275       // coefficients until we have the new normalization
01276       // coefficients, which we need in order to do the update.
01277       rawProject (X_in, Q, C2);
01278       
01279       // Coefficients for (re)normalization of X_in
01280       RCP<mat_type> B2 (new mat_type (ncols_X, ncols_X));
01281 
01282       // Normalize X_in (into X_out, if working out of place)
01283       const int secondPassRank = outOfPlace ? 
01284         normalizeOutOfPlace (X_in, X_out, B2) : 
01285         normalize (X_in, B2);
01286       rank = secondPassRank; // Current rank of X
01287 
01288       // Update normalization coefficients.  We begin with
01289       // copying B, since the BLAS doesn't like aliased input
01290       // arguments.
01291       mat_type B_copy (Copy, *B, B->numRows(), B->numCols());
01292       // B := B2 * B
01293       const int err = B->multiply (NO_TRANS, NO_TRANS, SCT::one(), 
01294            *B2, B_copy, SCT::zero());
01295       TEST_FOR_EXCEPTION(err != 0, std::logic_error, 
01296              "Teuchos::SerialDenseMatrix::multiply "
01297              "returned err = " << err << " != 0");
01298       // Update the block coefficients from the projection step.
01299       // We use B_copy for this (a copy of the original B, the
01300       // first-pass normalization coefficients).
01301       for (int k = 0; k < num_Q_blocks; ++k)
01302         {
01303     mat_type C_k_copy (Copy, *C[k], C[k]->numRows(), C[k]->numCols());
01304 
01305     // C[k] := C2[k]*B_copy + C[k] (where B_copy is a copy
01306     // of the original B, the first-pass normalization
01307     // coefficients).
01308     const int err = C[k]->multiply (NO_TRANS, NO_TRANS, SCT::one(), 
01309             *C2[k], B_copy, SCT::one());
01310     TEST_FOR_EXCEPTION(err != 0, std::logic_error, 
01311            "Teuchos::SerialDenseMatrix::multiply "
01312            "returned err = " << err << " != 0");
01313         }
01314       // Compute post-second-pass (pre-normalization) norms,
01315       // using B2 (the coefficients from the second normalize()
01316       // invocation) in the same way as with B before.
01317       for (int j = 0; j < rank; ++j)
01318         {
01319     const Scalar* const B2_j = &(*B2)(0,j);
01320     normsAfterSecondPass[j] = blas.NRM2 (rank, B2_j, 1);
01321         }
01322       // Test whether any of the norms dropped below the
01323       // reorthogonalization threshold.  If so, it's an
01324       // orthogonalization fault, which requires expensive
01325       // recovery.
01326       reorthogFault = false;
01327       for (int j = 0; j < rank; ++j)
01328         {
01329     const magnitude_type relativeLowerBound = 
01330       blockReorthogThreshold() * normsAfterFirstPass[j];
01331     if (normsAfterSecondPass[j] < relativeLowerBound)
01332       {
01333         reorthogFault = true; 
01334         faultIndices.push_back (j);
01335       }
01336         }
01337     } // if (reorthogonalize) // reorthogonalization pass
01338 
01339   if (reorthogFault)
01340     {
01341       if (throwOnReorthogFault_)
01342         raiseReorthogFault (normsAfterFirstPass, normsAfterSecondPass, faultIndices);
01343       else 
01344         {
01345     // NOTE (mfh 19 Jan 2011) We could handle the
01346     // reorthogonalization fault here by slowly
01347     // reorthogonalizing, one vector at a time, the
01348     // offending vectors of X.  However, we choose not to
01349     // implement this for now.  If it becomes a problem,
01350     // let us know and we will prioritize implementing
01351     // this.
01352     TEST_FOR_EXCEPTION(true, std::logic_error, 
01353            "TsqrOrthoManagerImpl has not yet implemented"
01354            " recovery from an orthogonalization fault.");
01355         }
01356     }
01357       } // if (reorthogonalizeBlocks_)
01358     return rank;
01359   }
01360 
01361 
01362   template<class Scalar, class MV>
01363   void
01364   TsqrOrthoManagerImpl<Scalar, MV>::
01365   raiseReorthogFault (const std::vector<magnitude_type>& normsAfterFirstPass,
01366           const std::vector<magnitude_type>& normsAfterSecondPass,
01367           const std::vector<int>& faultIndices)
01368   {
01369     using std::endl;
01370     typedef std::vector<int>::size_type size_type;
01371     std::ostringstream os;
01372       
01373     os << "Orthogonalization fault at the following column(s) of X:" << endl;
01374     os << "Column\tNorm decrease factor" << endl;
01375     for (size_type k = 0; k < faultIndices.size(); ++k)
01376       {
01377   const int index = faultIndices[k];
01378   const magnitude_type decreaseFactor = 
01379     normsAfterSecondPass[index] / normsAfterFirstPass[index];
01380   os << index << "\t" << decreaseFactor << endl;
01381       }
01382     throw TsqrOrthoFault (os.str());
01383   }
01384 
01385   template<class Scalar, class MV>
01386   Teuchos::RCP<const Teuchos::ParameterList>
01387   TsqrOrthoManagerImpl<Scalar, MV>::getDefaultParameters ()
01388   {
01389     using Teuchos::ParameterList;
01390     using Teuchos::RCP;
01391     typedef Teuchos::ScalarTraits<magnitude_type> SCTM;
01392 
01393     // This part makes this class method non-reentrant.
01394     static RCP<ParameterList> params;
01395     if (! params.is_null())
01396       return params;
01397     params = Teuchos::parameterList();
01398 
01399     // 
01400     // TSQR parameters
01401     //
01402     params->set ("TsqrImpl", tsqr_adaptor_type::getDefaultParameters(), 
01403      "TSQR implementation parameters.");
01404     // 
01405     // Orthogonalization parameters
01406     //
01407     const bool defaultRandomizeNullSpace = true;
01408     params->set ("randomizeNullSpace", defaultRandomizeNullSpace, 
01409      "Whether to fill in null space vectors with random data.");
01410     //const bool defaultReorthogonalizeBlocks = false;
01411     const bool defaultReorthogonalizeBlocks = true;
01412     params->set ("reorthogonalizeBlocks", defaultReorthogonalizeBlocks,
01413      "Whether to do block reorthogonalization.");
01414     // This parameter corresponds to the "blk_tol_" parameter in
01415     // Belos' DGKSOrthoManager.  We choose the same default value.
01416     const magnitude_type defaultBlockReorthogThreshold = 
01417       magnitude_type(10) * SCTM::squareroot (SCTM::eps());
01418     params->set ("blockReorthogThreshold", defaultBlockReorthogThreshold, 
01419      "If reorthogonalizeBlocks==true, and if the norm of "
01420      "any column within a block decreases by this much or "
01421      "more after orthogonalization, we reorthogonalize.");
01422     // This parameter corresponds to the "sing_tol_" parameter in
01423     // Belos' DGKSOrthoManager.  We choose the same default value.
01424     const magnitude_type defaultRelativeRankTolerance = 
01425       magnitude_type(10) * SCTM::eps();
01426     // If the relative rank tolerance is zero, then we will always
01427     // declare blocks to be numerically full rank, as long as no
01428     // singular values are zero.
01429     params->set ("relativeRankTolerance", defaultRelativeRankTolerance,
01430      "Relative tolerance to determine the numerical rank of a "
01431      "block when normalizing.");
01432     const bool defaultThrowOnReorthogFault = true;
01433     // See Stewart's 2008 paper on block Gram-Schmidt for a
01434     // definition of "orthogonalization fault."
01435     params->set ("throwOnReorthogFault", defaultThrowOnReorthogFault,
01436      "Whether to throw an exception if an orthogonalization "
01437      "fault occurs.  This only matters if reorthogonalization "
01438      "is enabled (reorthogonalizeBlocks==true).");
01439     return params;
01440   }
01441 
01442 
01443   template<class Scalar, class MV>
01444   Teuchos::RCP<const Teuchos::ParameterList>
01445   TsqrOrthoManagerImpl<Scalar, MV>::getFastParameters ()
01446   {
01447     using Teuchos::ParameterList;
01448     using Teuchos::RCP;
01449     using Teuchos::rcp;
01450 
01451     // This part makes this class method non-reentrant.
01452     static RCP<ParameterList> params;
01453     if (params.is_null())
01454       {
01455   RCP<const ParameterList> defaultParams = getDefaultParameters();
01456   // Start with a clone of the default parameters
01457   params = rcp (new ParameterList (*defaultParams));
01458   
01459   // Disable reorthogonalization and randomization of the null
01460   // space basis.  Reorthogonalization tolerances don't matter,
01461   // since we aren't reorthogonalizing blocks in the fast
01462   // settings.  We can leave the default values.  Also,
01463   // (re)orthogonalization faults may only occur with
01464   // reorthogonalization, so we don't have to worry about the
01465   // "throwOnReorthogFault" setting.
01466   const bool randomizeNullSpace = false;
01467   params->set ("randomizeNullSpace", randomizeNullSpace);      
01468   const bool reorthogonalizeBlocks = false;
01469   params->set ("reorthogonalizeBlocks", reorthogonalizeBlocks);
01470       }
01471     return params;
01472   }
01473 
01474 
01475   template<class Scalar, class MV>
01476   int
01477   TsqrOrthoManagerImpl<Scalar, MV>::
01478   rawNormalize (MV& X, 
01479     MV& Q, 
01480     Teuchos::SerialDenseMatrix<int, Scalar>& B)
01481   {
01482     int rank;
01483     try {
01484       // This call only computes the QR factorization X = Q B.
01485       // It doesn't compute the rank of X.  That comes from
01486       // revealRank() below.
01487       tsqrAdaptor_->factorExplicit (X, Q, B);
01488       // This call will only modify *B if *B on input is not of full
01489       // numerical rank.
01490       rank = tsqrAdaptor_->revealRank (Q, B, relativeRankTolerance_);
01491     } catch (std::exception& e) {
01492       throw TsqrOrthoError (e.what()); // Toss the exception up the chain.
01493     }
01494     return rank;
01495   }
01496 
01497   template<class Scalar, class MV>
01498   int
01499   TsqrOrthoManagerImpl<Scalar, MV>::
01500   normalizeOne (MV& X, 
01501     Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > B) const
01502   {
01503     std::vector<magnitude_type> theNorm (1, SCTM::zero());
01504     MVT::MvNorm (X, theNorm);
01505     if (Teuchos::is_null(B))
01506       B = Teuchos::rcp (new mat_type (1, 1));
01507     (*B)(0,0) = theNorm[0];
01508     if (theNorm[0] == SCTM::zero())
01509       {
01510   // Make a view of the first column of Q, fill it with random
01511   // data, and normalize it.  Throw away the resulting norm.
01512   if (randomizeNullSpace_)
01513     {
01514       MVT::MvRandom(X);
01515       MVT::MvNorm (X, theNorm);
01516       if (theNorm[0] == SCTM::zero())
01517         throw TsqrOrthoError("normalizeOne(): a supposedly random "
01518            "vector has norm zero!");
01519       else
01520         {
01521     // NOTE (mfh 09 Nov 2010) I'm assuming that this
01522     // operation that implicitly converts from a
01523     // magnitude_type to a Scalar makes sense.
01524     const Scalar alpha = SCT::one() / theNorm[0];
01525     MVT::MvScale (X, alpha);
01526         }
01527     }
01528   return 0;
01529       }
01530     else
01531       {
01532   // NOTE (mfh 09 Nov 2010) I'm assuming that this operation
01533   // that implicitly converts from a magnitude_type to a Scalar
01534   // makes sense.
01535   const Scalar alpha = SCT::one() / theNorm[0];
01536   MVT::MvScale (X, alpha);
01537   return 1;
01538       }
01539   }
01540 
01541 
01542   template<class Scalar, class MV>
01543   void
01544   TsqrOrthoManagerImpl<Scalar, MV>::
01545   rawProject (MV& X, 
01546         Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
01547         Teuchos::ArrayView<Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > > C) const
01548   {
01549 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01550     Teuchos::TimeMonitor timerMonitorNormalize(*timerProject_);
01551 #endif // BELOS_TEUCHOS_TIME_MONITOR
01552 
01553     // "Modified Gram-Schmidt" version of Block Gram-Schmidt.
01554     const int num_Q_blocks = Q.size();
01555     for (int i = 0; i < num_Q_blocks; ++i)
01556       {
01557   // TEST_FOR_EXCEPTION(C[i].is_null(), std::logic_error,
01558   //       "TsqrOrthoManagerImpl::rawProject(): C[" 
01559   //       << i << "] is null");
01560   // TEST_FOR_EXCEPTION(Q[i].is_null(), std::logic_error,
01561   //       "TsqrOrthoManagerImpl::rawProject(): Q[" 
01562   //       << i << "] is null");
01563   mat_type& Ci = *C[i];
01564   const MV& Qi = *Q[i];
01565   innerProd (Qi, X, Ci);
01566   MVT::MvTimesMatAddMv (-SCT::one(), Qi, Ci, SCT::one(), X);
01567       }
01568   }
01569 
01570 
01571   template<class Scalar, class MV>
01572   void
01573   TsqrOrthoManagerImpl<Scalar, MV>::
01574   rawProject (MV& X, 
01575         const Teuchos::RCP<const MV>& Q, 
01576         const Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> >& C) const
01577   {
01578 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01579     Teuchos::TimeMonitor timerMonitorNormalize(*timerProject_);
01580 #endif // BELOS_TEUCHOS_TIME_MONITOR
01581 
01582     // Block Gram-Schmidt
01583     innerProd (*Q, X, *C);
01584     MVT::MvTimesMatAddMv (-SCT::one(), *Q, *C, SCT::one(), X);
01585   }
01586 
01587   template<class Scalar, class MV>
01588   int
01589   TsqrOrthoManagerImpl<Scalar, MV>::
01590   normalizeImpl (MV& X, 
01591      MV& Q, 
01592      Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > B, 
01593      const bool outOfPlace)
01594   {
01595     using Teuchos::Range1D;
01596     using Teuchos::RCP;
01597     using Teuchos::rcp;
01598     using Teuchos::tuple;
01599 
01600     using std::cerr;
01601     using std::endl;
01602 
01603     const int numCols = MVT::GetNumberVecs (X);
01604     if (numCols == 0)
01605       return 0; // Fast exit for an empty input matrix
01606 
01607     // We allow Q to have more columns than X, in which case we only
01608     // touch the first numCols columns of Q.
01609     TEST_FOR_EXCEPTION(MVT::GetNumberVecs(Q) < numCols, 
01610            std::invalid_argument, 
01611            "TsqrOrthoManagerImpl::normalizeImpl(X,Q,B): "
01612            "Q has " << MVT::GetNumberVecs(Q) << " columns, "
01613            "which is too few, since X has " << numCols 
01614            << " columns.");
01615     // TSQR wants a Q with the same number of columns as X, so have
01616     // it work on a nonconstant view of Q with the same number of
01617     // columns as X.
01618     RCP<MV> Q_view = MVT::CloneViewNonConst (Q, Range1D(0, numCols-1));
01619 
01620     // TSQR's rank-revealing part doesn't work unless B is provided.
01621     // If B is not provided, allocate a temporary B for use in TSQR.
01622     // If it is provided, adjust dimensions as necessary.  Adjusting
01623     // the dimensions of B may invalidate any data currently stored in
01624     // it, but that data will be overwritten anyway.
01625     if (B.is_null())
01626       B = rcp (new mat_type (numCols, numCols));
01627     else
01628       B->shape (numCols, numCols);
01629 
01630     if (false)
01631       {
01632   std::vector<magnitude_type> norms (numCols);
01633   MVT::MvNorm (X, norms);
01634   cerr << "Column norms of X before orthogonalization: ";
01635   for (typename std::vector<magnitude_type>::const_iterator iter = norms.begin(); 
01636        iter != norms.end(); ++iter)
01637     {
01638       cerr << *iter;
01639       if (iter+1 != norms.end())
01640         cerr << ", ";
01641     }
01642       }
01643 
01644     // Compute rank-revealing decomposition (in this case, TSQR of X
01645     // followed by SVD of the R factor and appropriate updating of the
01646     // resulting Q factor) of X.  X is modified in place and filled
01647     // with garbage, and Q_view contains the resulting explicit Q
01648     // factor.  Later, we will copy this back into X.
01649     //
01650     // The matrix *B will only be upper triangular if X is of full
01651     // numerical rank.
01652     const int rank = rawNormalize (X, *Q_view, *B);
01653 
01654     if (false)
01655       {
01656   std::vector<magnitude_type> norms (numCols);
01657   MVT::MvNorm (*Q_view, norms);
01658   cerr << "Column norms of Q_view after orthogonalization: ";
01659   for (typename std::vector<magnitude_type>::const_iterator iter = norms.begin(); 
01660        iter != norms.end(); ++iter)
01661     {
01662       cerr << *iter;
01663       if (iter+1 != norms.end())
01664         cerr << ", ";
01665     }
01666       }
01667 
01668     TEST_FOR_EXCEPTION(rank < 0 || rank > numCols, std::logic_error,
01669            "Belos::TsqrOrthoManagerImpl::normalizeImpl: "
01670            "rawNormalize() returned rank = " << rank << " for a "
01671            "matrix X with " << numCols << " columns.  Please report"
01672            " this bug to the Belos developers.");
01673     if (false && rank == 0)
01674       {
01675   // Sanity check: ensure that the columns of X are sufficiently
01676   // small for X to be reported as rank zero.
01677   const mat_type& B_ref = *B;
01678   std::vector<magnitude_type> norms (B_ref.numCols());
01679   for (typename mat_type::ordinalType j = 0; j < B_ref.numCols(); ++j)
01680     {
01681       typedef typename mat_type::scalarType mat_scalar_type;
01682       mat_scalar_type sumOfSquares = Teuchos::ScalarTraits<mat_scalar_type>::zero();
01683       for (typename mat_type::ordinalType i = 0; i <= j; ++i)
01684         {
01685     const mat_scalar_type B_ij = 
01686       Teuchos::ScalarTraits<mat_scalar_type>::magnitude (B_ref(i,j));
01687     sumOfSquares += B_ij*B_ij;
01688         }
01689       norms[j] = Teuchos::ScalarTraits<mat_scalar_type>::squareroot (sumOfSquares);
01690     }
01691   bool anyNonzero = false;
01692   typedef typename std::vector<magnitude_type>::const_iterator iter_type;
01693   for (iter_type it = norms.begin(); it != norms.end(); ++it)
01694     if (*it > relativeRankTolerance_)
01695       anyNonzero = true;
01696 
01697   using std::cerr;
01698   using std::endl;
01699   cerr << "Norms of columns of B after orthogonalization: ";
01700   for (typename mat_type::ordinalType j = 0; j < B_ref.numCols(); ++j)
01701     {
01702       cerr << norms[j];
01703       if (j != B_ref.numCols() - 1)
01704         cerr << ", ";
01705     }
01706   cerr << endl;
01707       }
01708 
01709     // If X is full rank or we don't want to replace its null space
01710     // basis with random vectors, then we're done.
01711     if (rank == numCols || ! randomizeNullSpace_)
01712       {
01713   // If we're supposed to be working in place in X, copy the
01714   // results back from Q_view into X.
01715   if (! outOfPlace)
01716     MVT::Assign (*Q_view, X);
01717   return rank;
01718       }
01719 
01720     // Replace the null space basis of X (in the last numCols-rank
01721     // columns of Q_view) with random data, project it against the
01722     // first rank columns of Q_view, and normalize.
01723     if (randomizeNullSpace_ && rank < numCols)
01724       {
01725   // Number of columns to fill with random data.
01726   const int nullSpaceNumCols = numCols - rank;
01727   // Inclusive range of indices of columns of X to fill with
01728   // random data.
01729   Range1D nullSpaceIndices (rank, numCols-1);
01730 
01731   // rawNormalize() wrote the null space basis vectors into
01732   // Q_view.  We continue to work in place in Q_view by
01733   // writing the random data there and (if there is a
01734   // nontrival column space) projecting in place against the
01735   // column space basis vectors (also in Q_view).
01736   RCP<MV> Q_null = MVT::CloneViewNonConst (*Q_view, nullSpaceIndices);
01737   // Replace the null space basis with random data.
01738   MVT::MvRandom (*Q_null); 
01739 
01740   // Make sure that the "random" data isn't all zeros.  This is
01741   // statistically nearly impossible, but we test for debugging
01742   // purposes.
01743   {
01744     std::vector<magnitude_type> norms (MVT::GetNumberVecs(*Q_null));
01745     MVT::MvNorm (*Q_null, norms);
01746 
01747     bool anyZero = false;
01748     typedef typename std::vector<magnitude_type>::const_iterator iter_type;
01749     for (iter_type it = norms.begin(); it != norms.end(); ++it)
01750       if (*it == SCTM::zero())
01751         anyZero = true;
01752 
01753     if (anyZero)
01754       {
01755         std::ostringstream os;
01756         os << "TsqrOrthoManagerImpl::normalizeImpl: "
01757     "We are being asked to randomize the null space, for a matrix "
01758     "with " << numCols << " columns and reported column rank "
01759      << rank << ".  The inclusive range of columns to fill with "
01760     "random data is [" << nullSpaceIndices.lbound() << "," 
01761      << nullSpaceIndices.ubound() << "].  After filling the null "
01762     "space vectors with random numbers, at least one of the vectors"
01763     " has norm zero.  Here are the norms of all the null space "
01764     "vectors: [";
01765         for (iter_type it = norms.begin(); it != norms.end(); ++it)
01766     {
01767       os << *it;
01768       if (it+1 != norms.end())
01769         os << ", ";
01770     }
01771         os << "].)  There is a tiny probability that this could happen "
01772     "randomly, but it is likely a bug.  Please report it to the "
01773     "Belos developers, especially if you are able to reproduce the "
01774     "behavior.";
01775         TEST_FOR_EXCEPTION(anyZero, TsqrOrthoError, os.str());
01776       }
01777   }
01778 
01779   if (rank > 0)
01780     {
01781       // Project the random data against the column space
01782       // basis of X, using a simple block projection ("Block
01783       // Classical Gram-Schmidt").  This is accurate because
01784       // we've already orthogonalized the column space basis
01785       // of X nearly to machine precision via a QR
01786       // factorization (TSQR) with accuracy comparable to
01787       // Householder QR.
01788       RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D(0, rank-1));
01789 
01790       // Temporary storage for projection coefficients.  We
01791       // don't need to keep them, since they represent the
01792       // null space basis (for which the coefficients are
01793       // logically zero).
01794       mat_ptr C_null (new mat_type (rank, nullSpaceNumCols));
01795       rawProject (*Q_null, Q_col, C_null);
01796     }
01797   // Normalize the projected random vectors, so that they are
01798   // mutually orthogonal (as well as orthogonal to the column
01799   // space basis of X).  We use X for the output of the
01800   // normalization: for out-of-place normalization (outOfPlace
01801   // == true), X is overwritten with "invalid values" anyway,
01802   // and for in-place normalization (outOfPlace == false), we
01803   // want the result to be in X anyway.
01804   RCP<MV> X_null = MVT::CloneViewNonConst (X, nullSpaceIndices);
01805   // Normalization coefficients for projected random vectors.
01806   // Will be thrown away.
01807   mat_type B_null (nullSpaceNumCols, nullSpaceNumCols);
01808   // Write the normalized vectors to X_null (in X).
01809   const int nullSpaceBasisRank = rawNormalize (*Q_null, *X_null, B_null);
01810     
01811   // It's possible, but unlikely, that X_null doesn't have
01812   // full rank (after the projection step).  We could
01813   // recursively fill in more random vectors until we finally
01814   // get a full rank matrix, but instead we just throw an
01815   // exception.
01816   //
01817   // NOTE (mfh 08 Nov 2010) Perhaps we should deal with this
01818   // case more elegantly.  Recursion might be one way to solve
01819   // it, but be sure to check that the recursion will terminate.
01820   // We could do this by supplying an additional argument to
01821   // rawNormalize, which is the null space basis rank from the
01822   // previous iteration.  The rank has to decrease each time, or
01823   // the recursion may go on forever.
01824   if (nullSpaceBasisRank < nullSpaceNumCols)
01825     {
01826       std::vector<magnitude_type> norms (MVT::GetNumberVecs(*X_null));
01827       MVT::MvNorm (*X_null, norms);
01828       std::ostringstream os;
01829       os << "TsqrOrthoManagerImpl::normalizeImpl: "
01830          << "We are being asked to randomize the null space, "
01831          << "for a matrix with " << numCols << " columns and "
01832          << "column rank " << rank << ".  After projecting and "
01833          << "normalizing the generated random vectors, they "
01834          << "only have rank " << nullSpaceBasisRank << ".  They"
01835          << " should have full rank " << nullSpaceNumCols 
01836          << ".  (The inclusive range of columns to fill with "
01837          << "random data is [" << nullSpaceIndices.lbound() 
01838          << "," << nullSpaceIndices.ubound() << "].  The "
01839          << "column norms of the resulting Q factor are: [";
01840       for (typename std::vector<magnitude_type>::size_type k = 0; 
01841      k < norms.size(); ++k)
01842         {
01843     os << norms[k];
01844     if (k != norms.size()-1)
01845       os << ", ";
01846         }
01847       os << "].)  There is a tiny probability that this could "
01848          << "happen randomly, but it is likely a bug.  Please "
01849          << "report it to the Belos developers, especially if "
01850          << "you are able to reproduce the behavior.";
01851 
01852       TEST_FOR_EXCEPTION(nullSpaceBasisRank < nullSpaceNumCols, 
01853                          TsqrOrthoError, os.str());
01854     }
01855   // If we're normalizing out of place, copy the X_null
01856   // vectors back into Q_null; the Q_col vectors are already
01857   // where they are supposed to be in that case.
01858   //
01859   // If we're normalizing in place, leave X_null alone (it's
01860   // already where it needs to be, in X), but copy Q_col back
01861   // into the first rank columns of X.
01862   if (outOfPlace)
01863     MVT::Assign (*X_null, *Q_null);
01864   else if (rank > 0)
01865     { // MVT::Assign() doesn't accept empty ranges of columns.
01866       RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D(0, rank-1));
01867       RCP<MV> X_col = MVT::CloneViewNonConst (X, Range1D(0, rank-1));
01868       MVT::Assign (*Q_col, *X_col);
01869     }
01870       }
01871     return rank;
01872   }
01873 
01874 
01875   template<class Scalar, class MV>
01876   void
01877   TsqrOrthoManagerImpl<Scalar, MV>::
01878   checkProjectionDims (int& ncols_X, 
01879            int& num_Q_blocks,
01880            int& ncols_Q_total,
01881            const MV& X, 
01882            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
01883   {
01884     // First assign to temporary values, so the function won't
01885     // commit any externally visible side effects unless it will
01886     // return normally (without throwing an exception).  (I'm being
01887     // cautious; MVT::GetNumberVecs() probably shouldn't have any
01888     // externally visible side effects, unless it is logging to a
01889     // file or something.)
01890     int the_ncols_X, the_num_Q_blocks, the_ncols_Q_total;
01891     the_num_Q_blocks = Q.size();
01892     the_ncols_X = MVT::GetNumberVecs (X);
01893 
01894     // Compute the total number of columns of all the Q[i] blocks.
01895     the_ncols_Q_total = 0;
01896     // You should be angry if your compiler doesn't support type
01897     // inference ("auto").  That's why I need this awful typedef.
01898     using Teuchos::ArrayView;
01899     using Teuchos::RCP;
01900     typedef typename ArrayView<RCP<const MV> >::const_iterator iter_type;
01901     for (iter_type it = Q.begin(); it != Q.end(); ++it)
01902       {
01903   const MV& Qi = **it;
01904   the_ncols_Q_total += MVT::GetNumberVecs (Qi);
01905       }
01906 
01907     // Commit temporary values to the output arguments.
01908     ncols_X = the_ncols_X;
01909     num_Q_blocks = the_num_Q_blocks;
01910     ncols_Q_total = the_ncols_Q_total;
01911   }
01912 
01913 
01914 } // namespace Belos
01915 
01916 #endif // __BelosTsqrOrthoManagerImpl_hpp
01917 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines