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 #ifndef __BelosTsqrOrthoManagerImpl_hpp
00046 #define __BelosTsqrOrthoManagerImpl_hpp
00047 
00048 #include "BelosConfigDefs.hpp" // HAVE_BELOS_TSQR
00049 #include "BelosMultiVecTraits.hpp"
00050 #include "BelosOrthoManager.hpp" // OrthoError, etc.
00051 
00052 #include "Teuchos_as.hpp"
00053 #include "Teuchos_LAPACK.hpp"
00054 #include "Teuchos_ParameterList.hpp"
00055 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
00056 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00057 #  include "Teuchos_TimeMonitor.hpp"
00058 #endif // BELOS_TEUCHOS_TIME_MONITOR
00059 #include <algorithm>
00060 #include <functional> // std::binary_function
00061 
00062 namespace Belos {
00063 
00067   class TsqrOrthoError : public OrthoError {
00068   public: 
00069     TsqrOrthoError (const std::string& what_arg) : 
00070       OrthoError (what_arg) {}
00071   };
00072 
00092   class TsqrOrthoFault : public OrthoError {
00093   public: 
00094     TsqrOrthoFault (const std::string& what_arg) : 
00095       OrthoError (what_arg) {}
00096   };
00097 
00130   template<class Scalar>
00131   class ReorthogonalizationCallback : 
00132     public std::binary_function<Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>, 
00133         Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>,
00134         void>
00135   {
00136   public:
00138     typedef Scalar scalar_type;
00143     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00144 
00149     virtual void
00150     operator() (Teuchos::ArrayView<magnitude_type> normsBeforeFirstPass,
00151     Teuchos::ArrayView<magnitude_type> normsAfterFirstPass) = 0;
00152   };
00153 
00154 
00187   template<class Scalar, class MV>
00188   class TsqrOrthoManagerImpl : 
00189     public Teuchos::ParameterListAcceptorDefaultBase {
00190   public:
00191     typedef Scalar scalar_type;
00192     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00193     typedef MV multivector_type;
00196     typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00197     typedef Teuchos::RCP<mat_type> mat_ptr;
00198 
00199   private:
00200     typedef Teuchos::ScalarTraits<Scalar> SCT;
00201     typedef Teuchos::ScalarTraits<magnitude_type> SCTM;
00202     typedef MultiVecTraits<Scalar, MV> MVT;
00203     typedef MultiVecTraitsExt<Scalar, MV> MVText;
00204     typedef typename MVT::tsqr_adaptor_type tsqr_adaptor_type;
00205 
00206   public:
00214     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters () const;
00215 
00217     void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params);
00218 
00229     Teuchos::RCP<const Teuchos::ParameterList> getFastParameters ();
00230 
00247     TsqrOrthoManagerImpl (const Teuchos::RCP<Teuchos::ParameterList>& params,
00248         const std::string& label);
00249 
00254     TsqrOrthoManagerImpl (const std::string& label);
00255     
00275     void 
00276     setReorthogonalizationCallback (const Teuchos::RCP<ReorthogonalizationCallback<Scalar> >& callback)
00277     {
00278       reorthogCallback_ = callback;
00279     }
00280 
00288     void setLabel (const std::string& label) { 
00289       if (label != label_) {
00290   label_ = label; 
00291 
00292 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00293   clearTimer (label, "All orthogonalization");
00294   clearTimer (label, "Projection");
00295   clearTimer (label, "Normalization");
00296 
00297   timerOrtho_ = makeTimer (label, "All orthogonalization");
00298   timerProject_ = makeTimer (label, "Projection");
00299   timerNormalize_ = makeTimer (label, "Normalization");
00300 #endif // BELOS_TEUCHOS_TIME_MONITOR  
00301       }
00302     }
00303 
00305     const std::string& getLabel () const { return label_; }
00306 
00315     void 
00316     innerProd (const MV& X, const MV& Y, mat_type& Z) const
00317     {
00318       MVT::MvTransMv (SCT::one(), X, Y, Z);
00319     }
00320 
00338     void
00339     norm (const MV& X, std::vector<magnitude_type>& normVec) const;
00340 
00350     void 
00351     project (MV& X, 
00352        Teuchos::Array<mat_ptr> C, 
00353        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q);
00354 
00368     int normalize (MV& X, mat_ptr B);
00369 
00388     int 
00389     normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B);
00390 
00404     int 
00405     projectAndNormalize (MV &X,
00406        Teuchos::Array<mat_ptr> C,
00407        mat_ptr B,
00408        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
00409     {
00410       // "false" means we work on X in place.  The second argument is
00411       // not read or written in that case.
00412       return projectAndNormalizeImpl (X, X, false, C, B, Q);
00413     }
00414 
00434     int 
00435     projectAndNormalizeOutOfPlace (MV& X_in, 
00436            MV& X_out,
00437            Teuchos::Array<mat_ptr> C,
00438            mat_ptr B,
00439            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
00440     {
00441       // "true" means we work on X_in out of place, writing the
00442       // results into X_out.
00443       return projectAndNormalizeImpl (X_in, X_out, true, C, B, Q);
00444     }
00445     
00450     magnitude_type 
00451     orthonormError (const MV &X) const
00452     {
00453       const Scalar ONE = SCT::one();
00454       const int ncols = MVT::GetNumberVecs(X);
00455       mat_type XTX (ncols, ncols);
00456       innerProd (X, X, XTX);
00457       for (int k = 0; k < ncols; ++k) {
00458   XTX(k,k) -= ONE;
00459       }
00460       return XTX.normFrobenius();
00461     }
00462 
00464     magnitude_type 
00465     orthogError (const MV &X1, 
00466      const MV &X2) const
00467     {
00468       const int ncols_X1 = MVT::GetNumberVecs (X1);
00469       const int ncols_X2 = MVT::GetNumberVecs (X2);
00470       mat_type X1_T_X2 (ncols_X1, ncols_X2);
00471       innerProd (X1, X2, X1_T_X2);
00472       return X1_T_X2.normFrobenius();
00473     }
00474 
00478     magnitude_type blockReorthogThreshold() const { return blockReorthogThreshold_; }
00479 
00482     magnitude_type relativeRankTolerance() const { return relativeRankTolerance_; }
00483 
00484   private:
00486     Teuchos::RCP<Teuchos::ParameterList> params_;
00487 
00489     mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
00490 
00492     std::string label_;
00493 
00495     tsqr_adaptor_type tsqrAdaptor_;
00496 
00506     Teuchos::RCP<MV> Q_;
00507 
00509     magnitude_type eps_;
00510 
00514     bool randomizeNullSpace_;
00515 
00521     bool reorthogonalizeBlocks_;
00522 
00526     bool throwOnReorthogFault_;
00527 
00529     magnitude_type blockReorthogThreshold_;
00530 
00532     magnitude_type relativeRankTolerance_;
00533 
00540     bool forceNonnegativeDiagonal_;
00541 
00542 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00543 
00544     Teuchos::RCP<Teuchos::Time> timerOrtho_;
00545 
00547     Teuchos::RCP<Teuchos::Time> timerProject_;
00548 
00550     Teuchos::RCP<Teuchos::Time> timerNormalize_;
00551 #endif // BELOS_TEUCHOS_TIME_MONITOR
00552 
00554     Teuchos::RCP<ReorthogonalizationCallback<Scalar> > reorthogCallback_;
00555 
00556 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00557 
00558 
00559 
00560 
00561 
00562 
00563 
00564 
00565     static Teuchos::RCP<Teuchos::Time>
00566     makeTimer (const std::string& prefix, 
00567          const std::string& timerName)
00568     {
00569       const std::string timerLabel = 
00570   prefix.empty() ? timerName : (prefix + ": " + timerName);
00571       return Teuchos::TimeMonitor::getNewCounter (timerLabel);
00572     }
00573 
00579     void
00580     clearTimer (const std::string& prefix, 
00581     const std::string& timerName)
00582     {
00583       const std::string timerLabel = 
00584   prefix.empty() ? timerName : (prefix + ": " + timerName);
00585       Teuchos::TimeMonitor::clearCounter (timerLabel);
00586     }
00587 #endif // BELOS_TEUCHOS_TIME_MONITOR
00588 
00590     void
00591     raiseReorthogFault (const std::vector<magnitude_type>& normsAfterFirstPass,
00592       const std::vector<magnitude_type>& normsAfterSecondPass,
00593       const std::vector<int>& faultIndices);
00594 
00604     void
00605     checkProjectionDims (int& ncols_X, 
00606        int& num_Q_blocks,
00607        int& ncols_Q_total,
00608        const MV& X, 
00609        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00610 
00621     void
00622     allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C, 
00623             Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
00624             const MV& X,
00625             const bool attemptToRecycle = true) const;
00626 
00635     int 
00636     projectAndNormalizeImpl (MV& X_in, 
00637            MV& X_out,
00638            const bool outOfPlace,
00639            Teuchos::Array<mat_ptr> C,
00640            mat_ptr B,
00641            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q);
00642 
00648     void
00649     rawProject (MV& X, 
00650     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
00651     Teuchos::ArrayView<mat_ptr> C) const;
00652 
00654     void
00655     rawProject (MV& X, 
00656     const Teuchos::RCP<const MV>& Q, 
00657     const mat_ptr& C) const;
00658 
00685     int rawNormalize (MV& X, MV& Q, mat_type& B);
00686 
00704     int normalizeOne (MV& X, mat_ptr B) const;
00705 
00733     int normalizeImpl (MV& X, MV& Q, mat_ptr B, const bool outOfPlace);
00734   };
00735 
00736   template<class Scalar, class MV>
00737   void
00738   TsqrOrthoManagerImpl<Scalar, MV>::
00739   setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params)
00740   {
00741     using Teuchos::ParameterList;
00742     using Teuchos::parameterList;
00743     using Teuchos::RCP;
00744     using Teuchos::sublist;
00745     typedef magnitude_type M; // abbreviation.
00746 
00747     RCP<const ParameterList> defaultParams = getValidParameters ();
00748     // Sublist of TSQR implementation parameters; to get below.
00749     RCP<ParameterList> tsqrParams;
00750 
00751     RCP<ParameterList> theParams;
00752     if (params.is_null()) {
00753       theParams = parameterList (*defaultParams);
00754     } else {
00755       theParams = params;
00756 
00757       // Don't call validateParametersAndSetDefaults(); we prefer to
00758       // ignore parameters that we don't recognize, at least for now.
00759       // However, we do fill in missing parameters with defaults.
00760 
00761       randomizeNullSpace_ = 
00762   theParams->get<bool> ("randomizeNullSpace", 
00763             defaultParams->get<bool> ("randomizeNullSpace"));
00764       reorthogonalizeBlocks_ = 
00765   theParams->get<bool> ("reorthogonalizeBlocks", 
00766             defaultParams->get<bool> ("reorthogonalizeBlocks"));
00767       throwOnReorthogFault_ = 
00768   theParams->get<bool> ("throwOnReorthogFault", 
00769             defaultParams->get<bool> ("throwOnReorthogFault"));
00770       blockReorthogThreshold_ = 
00771   theParams->get<M> ("blockReorthogThreshold",
00772          defaultParams->get<M> ("blockReorthogThreshold"));
00773       relativeRankTolerance_ = 
00774   theParams->get<M> ("relativeRankTolerance", 
00775          defaultParams->get<M> ("relativeRankTolerance"));
00776       forceNonnegativeDiagonal_ = 
00777   theParams->get<bool> ("forceNonnegativeDiagonal", 
00778             defaultParams->get<bool> ("forceNonnegativeDiagonal"));
00779 
00780       // Get the sublist of TSQR implementation parameters.  Use the
00781       // default sublist if one isn't provided.
00782       if (! theParams->isSublist ("TSQR implementation")) {
00783   theParams->set ("TSQR implementation", 
00784       defaultParams->sublist ("TSQR implementation"));
00785       }
00786       tsqrParams = sublist (theParams, "TSQR implementation", true);
00787     }
00788 
00789     // Send the TSQR implementation parameters to the TSQR adaptor.
00790     tsqrAdaptor_.setParameterList (tsqrParams);
00791 
00792     // Save the input parameter list.
00793     setMyParamList (theParams);
00794   }
00795  
00796   template<class Scalar, class MV>
00797   TsqrOrthoManagerImpl<Scalar, MV>::
00798   TsqrOrthoManagerImpl (const Teuchos::RCP<Teuchos::ParameterList>& params,
00799       const std::string& label) :
00800     label_ (label),
00801     Q_ (Teuchos::null),               // Initialized on demand
00802     eps_ (SCTM::eps()),               // Machine precision
00803     randomizeNullSpace_ (true),
00804     reorthogonalizeBlocks_ (true),
00805     throwOnReorthogFault_ (false),
00806     blockReorthogThreshold_ (0),
00807     relativeRankTolerance_ (0),
00808     forceNonnegativeDiagonal_ (false)
00809   {
00810     setParameterList (params); // This also sets tsqrAdaptor_'s parameters.
00811 
00812 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00813     timerOrtho_ = makeTimer (label, "All orthogonalization");
00814     timerProject_ = makeTimer (label, "Projection");
00815     timerNormalize_ = makeTimer (label, "Normalization");
00816 #endif // BELOS_TEUCHOS_TIME_MONITOR
00817   }
00818 
00819   template<class Scalar, class MV>
00820   TsqrOrthoManagerImpl<Scalar, MV>::
00821   TsqrOrthoManagerImpl (const std::string& label) :
00822     label_ (label),
00823     Q_ (Teuchos::null),               // Initialized on demand
00824     eps_ (SCTM::eps()),               // Machine precision
00825     randomizeNullSpace_ (true),
00826     reorthogonalizeBlocks_ (true),
00827     throwOnReorthogFault_ (false),
00828     blockReorthogThreshold_ (0),
00829     relativeRankTolerance_ (0), 
00830     forceNonnegativeDiagonal_ (false) 
00831   {
00832     setParameterList (Teuchos::null); // Set default parameters.
00833 
00834 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00835     timerOrtho_ = makeTimer (label, "All orthogonalization");
00836     timerProject_ = makeTimer (label, "Projection");
00837     timerNormalize_ = makeTimer (label, "Normalization");
00838 #endif // BELOS_TEUCHOS_TIME_MONITOR
00839   }
00840 
00841   template<class Scalar, class MV>
00842   void
00843   TsqrOrthoManagerImpl<Scalar, MV>::
00844   norm (const MV& X, std::vector<magnitude_type>& normVec) const
00845   {
00846     const int numCols = MVT::GetNumberVecs (X);
00847     // std::vector<T>::size_type is unsigned; int is signed.  Mixed
00848     // unsigned/signed comparisons trigger compiler warnings.
00849     if (normVec.size() < static_cast<size_t>(numCols))
00850       normVec.resize (numCols); // Resize normvec if necessary.
00851     MVT::MvNorm (X, normVec);
00852   }
00853 
00854   template<class Scalar, class MV>
00855   void
00856   TsqrOrthoManagerImpl<Scalar, MV>::project (MV& X, 
00857                Teuchos::Array<mat_ptr> C,
00858                Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
00859   {
00860 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00861     // "Projection" only happens in rawProject(), so we only time
00862     // projection inside rawProject().  However, we count the time
00863     // spend in project() as part of the whole orthogonalization.
00864     //
00865     // If project() is called from projectAndNormalize(), the
00866     // TimeMonitor won't start timerOrtho_, because it is already
00867     // running in projectAndNormalize().
00868     Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
00869 #endif // BELOS_TEUCHOS_TIME_MONITOR
00870 
00871     int ncols_X, num_Q_blocks, ncols_Q_total;
00872     checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X, Q);
00873     // Test for quick exit: any dimension of X is zero, or there are
00874     // zero Q blocks, or the total number of columns of the Q blocks
00875     // is zero.
00876     if (ncols_X == 0 || num_Q_blocks == 0 || ncols_Q_total == 0)
00877       return;
00878 
00879     // Make space for first-pass projection coefficients
00880     allocateProjectionCoefficients (C, Q, X, true);
00881 
00882     // We only use columnNormsBefore and compute pre-projection column
00883     // norms if doing block reorthogonalization.
00884     std::vector<magnitude_type> columnNormsBefore (ncols_X, magnitude_type(0));
00885     if (reorthogonalizeBlocks_)
00886       MVT::MvNorm (X, columnNormsBefore);
00887 
00888     // Project (first block orthogonalization step): 
00889     // C := Q^* X, X := X - Q C.
00890     rawProject (X, Q, C); 
00891 
00892     // If we are doing block reorthogonalization, reorthogonalize X if
00893     // necessary.
00894     if (reorthogonalizeBlocks_) {
00895       std::vector<magnitude_type> columnNormsAfter (ncols_X, magnitude_type(0));
00896       MVT::MvNorm (X, columnNormsAfter);
00897   
00898       // Relative block reorthogonalization threshold.
00899       const magnitude_type relThres = blockReorthogThreshold();
00900       // Reorthogonalize X if any of its column norms decreased by a
00901       // factor more than the block reorthogonalization threshold.
00902       // Don't bother trying to subset the columns; that will make the
00903       // columns noncontiguous and thus hinder BLAS 3 optimizations.
00904       bool reorthogonalize = false;
00905       for (int j = 0; j < ncols_X; ++j) {
00906   if (columnNormsAfter[j] < relThres * columnNormsBefore[j]) {
00907     reorthogonalize = true;
00908     break;
00909   }
00910       }
00911       if (reorthogonalize) {
00912   // Notify the caller via callback about the need for
00913   // reorthogonalization.
00914   if (! reorthogCallback_.is_null()) {
00915     reorthogCallback_->operator() (Teuchos::arrayViewFromVector (columnNormsBefore), 
00916            Teuchos::arrayViewFromVector (columnNormsAfter));
00917   }
00918   // Second-pass projection coefficients
00919   Teuchos::Array<mat_ptr> C2;
00920   allocateProjectionCoefficients (C2, Q, X, false);
00921 
00922   // Perform the second projection pass:
00923   // C2 = Q' X, X = X - Q*C2
00924   rawProject (X, Q, C2);
00925   // Update the projection coefficients
00926   for (int k = 0; k < num_Q_blocks; ++k)
00927     *C[k] += *C2[k];
00928       }
00929     }
00930   }
00931 
00932 
00933   template<class Scalar, class MV>
00934   int 
00935   TsqrOrthoManagerImpl<Scalar, MV>::normalize (MV& X, mat_ptr B)
00936   {
00937     using Teuchos::Range1D;
00938     using Teuchos::RCP;
00939 
00940 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00941     Teuchos::TimeMonitor timerMonitorNormalize(*timerNormalize_);
00942     // If normalize() is called internally -- i.e., called from
00943     // projectAndNormalize() -- the timer will not be started or 
00944     // stopped, because it is already running.  TimeMonitor handles
00945     // recursive invocation by doing nothing.
00946     Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
00947 #endif // BELOS_TEUCHOS_TIME_MONITOR
00948 
00949     // MVT returns int for this, even though the "local ordinal
00950     // type" of the MV may be some other type (for example,
00951     // Tpetra::MultiVector<double, int32_t, int64_t, ...>).
00952     const int numCols = MVT::GetNumberVecs (X);
00953 
00954     // This special case (for X having only one column) makes
00955     // TsqrOrthoManagerImpl equivalent to Modified Gram-Schmidt
00956     // orthogonalization with conditional full reorthogonalization,
00957     // if all multivector inputs have only one column.  It also
00958     // avoids allocating Q_ scratch space and copying data when we
00959     // don't need to invoke TSQR at all.
00960     if (numCols == 1) {
00961       return normalizeOne (X, B);
00962     }
00963 
00964     // We use Q_ as scratch space for the normalization, since TSQR
00965     // requires a scratch multivector (it can't factor in place).  Q_
00966     // should come from a vector space compatible with X's vector
00967     // space, and Q_ should have at least as many columns as X.
00968     // Otherwise, we have to reallocate.  We also have to allocate
00969     // (not "re-") Q_ if we haven't allocated it before.  (We can't
00970     // allocate Q_ until we have some X, so we need a multivector as
00971     // the "prototype.")
00972     //
00973     // NOTE (mfh 11 Jan 2011) We only increase the number of columsn
00974     // in Q_, never decrease.  This is OK for typical uses of TSQR,
00975     // but you might prefer different behavior in some cases.
00976     //
00977     // NOTE (mfh 10 Mar 2011) We should only reuse the scratch space
00978     // Q_ if X and Q_ have compatible data distributions.  However,
00979     // Belos' current MultiVecTraits interface does not let us check
00980     // for this.  Thus, we can only check whether X and Q_ have the
00981     // same number of rows.  This will behave correctly for the common
00982     // case in Belos that all multivectors with the same number of
00983     // rows have the same data distribution.
00984     //
00985     // The specific MV implementation may do more checks than this on
00986     // its own and throw an exception if X and Q_ are not compatible,
00987     // but it may not.  If you find that recycling the Q_ space causes
00988     // troubles, you may consider modifying the code below to
00989     // reallocate Q_ for every X that comes in.  
00990     if (Q_.is_null() || 
00991   MVText::GetGlobalLength(*Q_) != MVText::GetGlobalLength(X) ||
00992   numCols > MVT::GetNumberVecs (*Q_)) {
00993       Q_ = MVT::Clone (X, numCols);
00994     }
00995 
00996     // normalizeImpl() wants the second MV argument to have the same
00997     // number of columns as X.  To ensure this, we pass it a view of
00998     // Q_ if Q_ has more columns than X.  (This is possible if we
00999     // previously called normalize() with a different multivector,
01000     // since we never reallocate Q_ if it has more columns than
01001     // necessary.)
01002     if (MVT::GetNumberVecs(*Q_) == numCols) {
01003       return normalizeImpl (X, *Q_, B, false);
01004     } else {
01005       RCP<MV> Q_view = MVT::CloneViewNonConst (*Q_, Range1D(0, numCols-1));
01006       return normalizeImpl (X, *Q_view, B, false);
01007     }
01008   }
01009 
01010   template<class Scalar, class MV>
01011   void
01012   TsqrOrthoManagerImpl<Scalar, MV>::
01013   allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C, 
01014           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
01015           const MV& X,
01016           const bool attemptToRecycle) const
01017   {
01018     const int num_Q_blocks = Q.size();
01019     const int ncols_X = MVT::GetNumberVecs (X);
01020     C.resize (num_Q_blocks);
01021     if (attemptToRecycle)
01022       {
01023   for (int i = 0; i < num_Q_blocks; ++i) 
01024     {
01025       const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
01026       // Create a new C[i] if necessary, otherwise resize if
01027       // necessary, otherwise fill with zeros.
01028       if (C[i].is_null())
01029         C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
01030       else 
01031         {
01032     mat_type& Ci = *C[i];
01033     if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X)
01034       Ci.shape (ncols_Qi, ncols_X);
01035     else
01036       Ci.putScalar (SCT::zero());
01037         }
01038     }
01039       }
01040     else
01041       {
01042   for (int i = 0; i < num_Q_blocks; ++i) 
01043     {
01044       const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
01045       C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
01046     }
01047       }
01048   }
01049 
01050   template<class Scalar, class MV>
01051   int 
01052   TsqrOrthoManagerImpl<Scalar, MV>::
01053   normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B)
01054   {
01055 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01056     Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
01057     Teuchos::TimeMonitor timerMonitorNormalize(*timerNormalize_);
01058 #endif // BELOS_TEUCHOS_TIME_MONITOR
01059 
01060     const int numVecs = MVT::GetNumberVecs(X);
01061     if (numVecs == 0) {
01062       return 0; // Nothing to do.
01063     } else if (numVecs == 1) {
01064       // Special case for a single column; scale and copy over.
01065       using Teuchos::Range1D;
01066       using Teuchos::RCP;
01067       using Teuchos::rcp;
01068 
01069       // Normalize X in place (faster than TSQR for one column).
01070       const int rank = normalizeOne (X, B);
01071       // Copy results to first column of Q.
01072       RCP<MV> Q_0 = MVT::CloneViewNonConst (Q, Range1D(0,0));
01073       MVT::Assign (X, *Q_0);
01074       return rank;
01075     } else {
01076       // The "true" argument to normalizeImpl() means the output
01077       // vectors go into Q, and the contents of X are overwritten with
01078       // invalid values.
01079       return normalizeImpl (X, Q, B, true);
01080     }
01081   }
01082 
01083   template<class Scalar, class MV>
01084   int 
01085   TsqrOrthoManagerImpl<Scalar, MV>::
01086   projectAndNormalizeImpl (MV& X_in, 
01087          MV& X_out, // Only written if outOfPlace==false.
01088          const bool outOfPlace,
01089          Teuchos::Array<mat_ptr> C,
01090          mat_ptr B,
01091          Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
01092   {
01093     using Teuchos::Range1D;
01094     using Teuchos::RCP;
01095     using Teuchos::rcp;
01096 
01097 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01098     // Projection is only timed in rawProject(), and normalization is
01099     // only timed in normalize() and normalizeOutOfPlace().
01100     Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
01101 #endif // BELOS_TEUCHOS_TIME_MONITOR
01102 
01103     if (outOfPlace) {
01104       // Make sure that X_out has at least as many columns as X_in.
01105       TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(X_out) < MVT::GetNumberVecs(X_in),
01106        std::invalid_argument, 
01107        "Belos::TsqrOrthoManagerImpl::"
01108        "projectAndNormalizeImpl(..., outOfPlace=true, ...):"
01109        "X_out has " << MVT::GetNumberVecs(X_out) 
01110        << " columns, but X_in has "
01111        << MVT::GetNumberVecs(X_in) << " columns.");
01112     }
01113     // Fetch dimensions of X_in and Q, and allocate space for first-
01114     // and second-pass projection coefficients (C resp. C2).
01115     int ncols_X, num_Q_blocks, ncols_Q_total;
01116     checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X_in, Q);
01117 
01118     // Test for quick exit: if any dimension of X is zero.
01119     if (ncols_X == 0) {
01120       return 0;
01121     }
01122     // If there are zero Q blocks or zero Q columns, just normalize!
01123     if (num_Q_blocks == 0 || ncols_Q_total == 0) {
01124       if (outOfPlace) {
01125   return normalizeOutOfPlace (X_in, X_out, B);
01126       } else {
01127   return normalize (X_in, B);
01128       }
01129     }
01130 
01131     // The typical case is that the entries of C have been allocated
01132     // before, so we attempt to recycle the allocations.  The call
01133     // below will reallocate if it cannot recycle.
01134     allocateProjectionCoefficients (C, Q, X_in, true);
01135 
01136     // If we are doing block reorthogonalization, then compute the
01137     // column norms of X before projecting for the first time.  This
01138     // will help us decide whether we need to reorthogonalize X.
01139     std::vector<magnitude_type> normsBeforeFirstPass (ncols_X, SCTM::zero());
01140     if (reorthogonalizeBlocks_) {
01141       MVT::MvNorm (X_in, normsBeforeFirstPass);
01142     }
01143 
01144     // First (Modified) Block Gram-Schmidt pass, in place in X_in.
01145     rawProject (X_in, Q, C);
01146 
01147     // Make space for the normalization coefficients.  This will
01148     // either be a freshly allocated matrix (if B is null), or a view
01149     // of the appropriately sized upper left submatrix of *B (if B is
01150     // not null).
01151     //
01152     // Note that if we let the normalize() routine allocate (in the
01153     // case that B is null), that storage will go away at the end of
01154     // normalize().  (This is because it passes the RCP by value, not
01155     // by reference.)
01156     mat_ptr B_out;
01157     if (B.is_null()) {
01158       B_out = rcp (new mat_type (ncols_X, ncols_X));
01159     } else {
01160       // Make sure that B is no smaller than numCols x numCols.
01161       TEUCHOS_TEST_FOR_EXCEPTION(B->numRows() < ncols_X || B->numCols() < ncols_X,
01162        std::invalid_argument,
01163        "normalizeOne: Input matrix B must be at "
01164        "least " << ncols_X << " x " << ncols_X 
01165        << ", but is instead " << B->numRows()
01166        << " x " << B->numCols() << ".");
01167       // Create a view of the ncols_X by ncols_X upper left
01168       // submatrix of *B.  TSQR will write the normalization
01169       // coefficients there.
01170       B_out = rcp (new mat_type (Teuchos::View, *B, ncols_X, ncols_X));
01171     }
01172 
01173     // Rank of X(_in) after first projection pass.  If outOfPlace,
01174     // this overwrites X_in with invalid values, and the results go in
01175     // X_out.  Otherwise, it's done in place in X_in.
01176     const int firstPassRank = outOfPlace ? 
01177       normalizeOutOfPlace (X_in, X_out, B_out) : 
01178       normalize (X_in, B_out);
01179     if (B.is_null()) {
01180       // The input matrix B is null, so assign B_out to it.  If B was
01181       // not null on input, then B_out is a view of *B, so we don't
01182       // have to do anything here.  Note that SerialDenseMatrix uses
01183       // raw pointers to store data and represent views, so we have to
01184       // be careful about scope.
01185       B = B_out;
01186     }
01187     int rank = firstPassRank; // Current rank of X.
01188 
01189     // If X was not full rank after projection and randomizeNullSpace_
01190     // is true, then normalize(OutOfPlace)() replaced the null space
01191     // basis of X with random vectors, and orthogonalized them against
01192     // the column space basis of X.  However, we still need to
01193     // orthogonalize the random vectors against the Q[i], after which
01194     // we will need to renormalize them.
01195     //
01196     // If outOfPlace, then we need to work in X_out (where
01197     // normalizeOutOfPlace() wrote the normalized vectors).
01198     // Otherwise, we need to work in X_in.
01199     //
01200     // Note: We don't need to keep the new projection coefficients,
01201     // since they are multiplied by the "small" part of B
01202     // corresponding to the null space of the original X.
01203     if (firstPassRank < ncols_X && randomizeNullSpace_) {
01204       const int numNullSpaceCols = ncols_X - firstPassRank;
01205       const Range1D nullSpaceIndices (firstPassRank, ncols_X - 1);
01206 
01207       // Space for projection coefficients (will be thrown away)
01208       Teuchos::Array<mat_ptr> C_null (num_Q_blocks);
01209       for (int k = 0; k < num_Q_blocks; ++k) {
01210   const int numColsQk = MVT::GetNumberVecs(*Q[k]);
01211   C_null[k] = rcp (new mat_type (numColsQk, numNullSpaceCols));
01212       }
01213       // Space for normalization coefficients (will be thrown away).
01214       RCP<mat_type> B_null (new mat_type (numNullSpaceCols, numNullSpaceCols));
01215 
01216       int randomVectorsRank;
01217       if (outOfPlace) {
01218   // View of the null space basis columns of X.
01219   // normalizeOutOfPlace() wrote them into X_out.
01220   RCP<MV> X_out_null = MVT::CloneViewNonConst (X_out, nullSpaceIndices);
01221   // Use X_in for scratch space.  Copy X_out_null into the
01222   // last few columns of X_in (X_in_null) and do projections
01223   // in there.  (This saves us a copy wen we renormalize
01224   // (out of place) back into X_out.)
01225   RCP<MV> X_in_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
01226   MVT::Assign (*X_out_null, *X_in_null);
01227   // Project the new random vectors against the Q blocks, and
01228   // renormalize the result into X_out_null.  
01229   rawProject (*X_in_null, Q, C_null);
01230   randomVectorsRank = normalizeOutOfPlace (*X_in_null, *X_out_null, B_null);
01231       } else {
01232   // View of the null space columns of X.  
01233   // They live in X_in.
01234   RCP<MV> X_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
01235   // Project the new random vectors against the Q blocks,
01236   // and renormalize the result (in place).
01237   rawProject (*X_null, Q, C_null);
01238   randomVectorsRank = normalize (*X_null, B_null);
01239       }
01240       // While unusual, it is still possible for the random data not
01241       // to be full rank after projection and normalization.  In that
01242       // case, we could try another set of random data and recurse as
01243       // necessary, but instead for now we just raise an exception.
01244       TEUCHOS_TEST_FOR_EXCEPTION(randomVectorsRank != numNullSpaceCols, 
01245        TsqrOrthoError, 
01246        "Belos::TsqrOrthoManagerImpl::projectAndNormalize"
01247        "OutOfPlace(): After projecting and normalizing the "
01248        "random vectors (used to replace the null space "
01249        "basis vectors from normalizing X), they have rank "
01250        << randomVectorsRank << ", but should have full "
01251        "rank " << numNullSpaceCols << ".");
01252     }
01253 
01254     // Whether or not X_in was full rank after projection, we still
01255     // might want to reorthogonalize against Q.
01256     if (reorthogonalizeBlocks_) {
01257       // We are only interested in the column space basis of X
01258       // resp. X_out.
01259       std::vector<magnitude_type> 
01260   normsAfterFirstPass (firstPassRank, SCTM::zero());
01261       std::vector<magnitude_type> 
01262   normsAfterSecondPass (firstPassRank, SCTM::zero());
01263 
01264       // Compute post-first-pass (pre-normalization) norms.  We could
01265       // have done that using MVT::MvNorm() on X_in after projecting,
01266       // but before the first normalization.  However, that operation
01267       // may be expensive.  It is also unnecessary: after calling
01268       // normalize(), the 2-norm of B(:,j) is the 2-norm of X_in(:,j)
01269       // before normalization, in exact arithmetic.
01270       //
01271       // NOTE (mfh 06 Nov 2010) This is one way that combining
01272       // projection and normalization into a single kernel --
01273       // projectAndNormalize() -- pays off.  In project(), we have to
01274       // compute column norms of X before and after projection.  Here,
01275       // we get them for free from the normalization coefficients.
01276       Teuchos::BLAS<int, Scalar> blas;
01277       for (int j = 0; j < firstPassRank; ++j) {
01278   const Scalar* const B_j = &(*B_out)(0,j);
01279   // Teuchos::BLAS::NRM2 returns a magnitude_type result on
01280   // Scalar inputs.
01281   normsAfterFirstPass[j] = blas.NRM2 (firstPassRank, B_j, 1);
01282       }
01283       // Test whether any of the norms dropped below the
01284       // reorthogonalization threshold.
01285       bool reorthogonalize = false;
01286       for (int j = 0; j < firstPassRank; ++j) {
01287   // If any column's norm decreased too much, mark this block
01288   // for reorthogonalization.  Note that this test will _not_
01289   // activate reorthogonalization if a column's norm before the
01290   // first project-and-normalize step was zero.  It _will_
01291   // activate reorthogonalization if the column's norm before
01292   // was not zero, but is zero now.
01293   const magnitude_type curThreshold = 
01294     blockReorthogThreshold() * normsBeforeFirstPass[j];
01295   if (normsAfterFirstPass[j] < curThreshold) {
01296     reorthogonalize = true; 
01297     break;
01298   }
01299       }
01300 
01301       // Notify the caller via callback about the need for
01302       // reorthogonalization.
01303       if (! reorthogCallback_.is_null()) {
01304   using Teuchos::arrayViewFromVector;
01305   (*reorthogCallback_) (arrayViewFromVector (normsBeforeFirstPass),
01306             arrayViewFromVector (normsAfterFirstPass));
01307       }
01308 
01309       // Perform another Block Gram-Schmidt pass if necessary.  "Twice
01310       // is enough" (Kahan's theorem) for a Krylov method, unless
01311       // (using Stewart's term) there is an "orthogonalization fault"
01312       // (indicated by reorthogFault).
01313       //
01314       // NOTE (mfh 07 Nov 2010) For now, we include the entire block
01315       // of X, including possible random data (that was already
01316       // projected and normalized above).  It might make more sense
01317       // just to process the first firstPassRank columns of X.
01318       // However, the resulting reorthogonalization should still be
01319       // correct regardless.
01320       bool reorthogFault = false;
01321       // Indices of X at which there was an orthogonalization fault.
01322       std::vector<int> faultIndices;
01323       if (reorthogonalize) {
01324   using Teuchos::Copy;
01325   using Teuchos::NO_TRANS;
01326 
01327   // If we're using out-of-place normalization, copy X_out
01328   // (results of first project and normalize pass) back into
01329   // X_in, for the second project and normalize pass.
01330   if (outOfPlace) {
01331     MVT::Assign (X_out, X_in);
01332   }
01333 
01334   // C2 is only used internally, so we know that we are
01335   // allocating fresh and not recycling allocations.  Stating
01336   // this lets us save time checking dimensions.
01337   Teuchos::Array<mat_ptr> C2;
01338   allocateProjectionCoefficients (C2, Q, X_in, false);
01339 
01340   // Block Gram-Schmidt (again).  Delay updating the block
01341   // coefficients until we have the new normalization
01342   // coefficients, which we need in order to do the update.
01343   rawProject (X_in, Q, C2);
01344       
01345   // Coefficients for (re)normalization of X_in.
01346   RCP<mat_type> B2 (new mat_type (ncols_X, ncols_X));
01347 
01348   // Normalize X_in (into X_out, if working out of place).
01349   const int secondPassRank = outOfPlace ? 
01350     normalizeOutOfPlace (X_in, X_out, B2) : 
01351     normalize (X_in, B2);
01352   rank = secondPassRank; // Current rank of X
01353 
01354   // Update normalization coefficients.  We begin with copying
01355   // B_out, since the BLAS' _GEMM routine doesn't let us alias
01356   // its input and output arguments.
01357   mat_type B_copy (Copy, *B_out, B_out->numRows(), B_out->numCols());
01358   // B_out := B2 * B_out (where input B_out is in B_copy).
01359   const int err = B_out->multiply (NO_TRANS, NO_TRANS, SCT::one(), 
01360            *B2, B_copy, SCT::zero());
01361   TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::logic_error, 
01362          "Teuchos::SerialDenseMatrix::multiply "
01363          "returned err = " << err << " != 0");
01364   // Update the block coefficients from the projection step.  We
01365   // use B_copy for this (a copy of B_out, the first-pass
01366   // normalization coefficients).
01367   for (int k = 0; k < num_Q_blocks; ++k) {
01368     mat_type C_k_copy (Copy, *C[k], C[k]->numRows(), C[k]->numCols());
01369 
01370     // C[k] := C2[k]*B_copy + C[k].
01371     const int err1 = C[k]->multiply (NO_TRANS, NO_TRANS, SCT::one(), 
01372             *C2[k], B_copy, SCT::one());
01373     TEUCHOS_TEST_FOR_EXCEPTION(err1 != 0, std::logic_error, 
01374            "Teuchos::SerialDenseMatrix::multiply "
01375            "returned err = " << err1 << " != 0");
01376   }
01377   // Compute post-second-pass (pre-normalization) norms, using
01378   // B2 (the coefficients from the second normalization step) in
01379   // the same way as with B_out before.
01380   for (int j = 0; j < rank; ++j) {
01381     const Scalar* const B2_j = &(*B2)(0,j);
01382     normsAfterSecondPass[j] = blas.NRM2 (rank, B2_j, 1);
01383   }
01384   // Test whether any of the norms dropped below the
01385   // reorthogonalization threshold.  If so, it's an
01386   // orthogonalization fault, which requires expensive recovery.
01387   reorthogFault = false;
01388   for (int j = 0; j < rank; ++j) {
01389     const magnitude_type relativeLowerBound = 
01390       blockReorthogThreshold() * normsAfterFirstPass[j];
01391     if (normsAfterSecondPass[j] < relativeLowerBound) {
01392       reorthogFault = true; 
01393       faultIndices.push_back (j);
01394     }
01395   }
01396       } // if (reorthogonalize) // reorthogonalization pass
01397 
01398       if (reorthogFault) {
01399   if (throwOnReorthogFault_) {
01400     raiseReorthogFault (normsAfterFirstPass, 
01401             normsAfterSecondPass, 
01402             faultIndices);
01403   } else {
01404     // NOTE (mfh 19 Jan 2011) We could handle the fault here by
01405     // slowly reorthogonalizing, one vector at a time, the
01406     // offending vectors of X.  However, we choose not to
01407     // implement this for now.  If it becomes a problem, let us
01408     // know and we will prioritize implementing this.
01409     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, 
01410            "TsqrOrthoManagerImpl has not yet implemented"
01411            " recovery from an orthogonalization fault.");
01412   }
01413       }
01414     } // if (reorthogonalizeBlocks_)
01415     return rank;
01416   }
01417 
01418 
01419   template<class Scalar, class MV>
01420   void
01421   TsqrOrthoManagerImpl<Scalar, MV>::
01422   raiseReorthogFault (const std::vector<magnitude_type>& normsAfterFirstPass,
01423           const std::vector<magnitude_type>& normsAfterSecondPass,
01424           const std::vector<int>& faultIndices)
01425   {
01426     using std::endl;
01427     typedef std::vector<int>::size_type size_type;
01428     std::ostringstream os;
01429       
01430     os << "Orthogonalization fault at the following column(s) of X:" << endl;
01431     os << "Column\tNorm decrease factor" << endl;
01432     for (size_type k = 0; k < faultIndices.size(); ++k) {
01433       const int index = faultIndices[k];
01434       const magnitude_type decreaseFactor = 
01435   normsAfterSecondPass[index] / normsAfterFirstPass[index];
01436       os << index << "\t" << decreaseFactor << endl;
01437     }
01438     throw TsqrOrthoFault (os.str());
01439   }
01440 
01441   template<class Scalar, class MV>
01442   Teuchos::RCP<const Teuchos::ParameterList>
01443   TsqrOrthoManagerImpl<Scalar, MV>::getValidParameters () const
01444   {
01445     using Teuchos::ParameterList;
01446     using Teuchos::parameterList;
01447     using Teuchos::RCP;
01448 
01449     if (defaultParams_.is_null()) {
01450       RCP<ParameterList> params = parameterList ("TsqrOrthoManagerImpl");
01451       //
01452       // TSQR parameters (set as a sublist).
01453       //
01454       params->set ("TSQR implementation", *(tsqrAdaptor_.getValidParameters()),
01455        "TSQR implementation parameters.");
01456       // 
01457       // Orthogonalization parameters
01458       //
01459       const bool defaultRandomizeNullSpace = true;
01460       params->set ("randomizeNullSpace", defaultRandomizeNullSpace, 
01461        "Whether to fill in null space vectors with random data.");
01462 
01463       const bool defaultReorthogonalizeBlocks = true;
01464       params->set ("reorthogonalizeBlocks", defaultReorthogonalizeBlocks,
01465        "Whether to do block reorthogonalization as necessary.");
01466 
01467       // This parameter corresponds to the "blk_tol_" parameter in
01468       // Belos' DGKSOrthoManager.  We choose the same default value.
01469       const magnitude_type defaultBlockReorthogThreshold = 
01470   magnitude_type(10) * SCTM::squareroot (SCTM::eps());
01471       params->set ("blockReorthogThreshold", defaultBlockReorthogThreshold, 
01472        "If reorthogonalizeBlocks==true, and if the norm of "
01473        "any column within a block decreases by this much or "
01474        "more after orthogonalization, we reorthogonalize.");
01475 
01476       // This parameter corresponds to the "sing_tol_" parameter in
01477       // Belos' DGKSOrthoManager.  We choose the same default value.
01478       const magnitude_type defaultRelativeRankTolerance = 
01479   Teuchos::as<magnitude_type>(10) * SCTM::eps();
01480 
01481       // If the relative rank tolerance is zero, then we will always
01482       // declare blocks to be numerically full rank, as long as no
01483       // singular values are zero.
01484       params->set ("relativeRankTolerance", defaultRelativeRankTolerance,
01485        "Relative tolerance to determine the numerical rank of a "
01486        "block when normalizing.");
01487 
01488       // See Stewart's 2008 paper on block Gram-Schmidt for a definition
01489       // of "orthogonalization fault."
01490       const bool defaultThrowOnReorthogFault = true;
01491       params->set ("throwOnReorthogFault", defaultThrowOnReorthogFault,
01492        "Whether to throw an exception if an orthogonalization "
01493        "fault occurs.  This only matters if reorthogonalization "
01494        "is enabled (reorthogonalizeBlocks==true).");
01495 
01496       const bool defaultForceNonnegativeDiagonal = false;
01497       params->set ("forceNonnegativeDiagonal", defaultForceNonnegativeDiagonal,
01498        "Whether to force the R factor produced by the normalization "
01499        "step to have a nonnegative diagonal.");
01500 
01501       defaultParams_ = params;
01502     }
01503     return defaultParams_;
01504   }
01505 
01506   template<class Scalar, class MV>
01507   Teuchos::RCP<const Teuchos::ParameterList>
01508   TsqrOrthoManagerImpl<Scalar, MV>::getFastParameters ()
01509   {
01510     using Teuchos::ParameterList;
01511     using Teuchos::RCP;
01512     using Teuchos::rcp;
01513 
01514     RCP<const ParameterList> defaultParams = getValidParameters();
01515     // Start with a clone of the default parameters.
01516     RCP<ParameterList> params = rcp (new ParameterList (*defaultParams));
01517   
01518     // Disable reorthogonalization and randomization of the null
01519     // space basis.  Reorthogonalization tolerances don't matter,
01520     // since we aren't reorthogonalizing blocks in the fast
01521     // settings.  We can leave the default values.  Also,
01522     // (re)orthogonalization faults may only occur with
01523     // reorthogonalization, so we don't have to worry about the
01524     // "throwOnReorthogFault" setting.
01525     const bool randomizeNullSpace = false;
01526     params->set ("randomizeNullSpace", randomizeNullSpace);      
01527     const bool reorthogonalizeBlocks = false;
01528     params->set ("reorthogonalizeBlocks", reorthogonalizeBlocks);
01529 
01530     return params;
01531   }
01532 
01533   template<class Scalar, class MV>
01534   int
01535   TsqrOrthoManagerImpl<Scalar, MV>::
01536   rawNormalize (MV& X, 
01537     MV& Q, 
01538     Teuchos::SerialDenseMatrix<int, Scalar>& B)
01539   {
01540     int rank;
01541     try {
01542       // This call only computes the QR factorization X = Q B.
01543       // It doesn't compute the rank of X.  That comes from
01544       // revealRank() below.
01545       tsqrAdaptor_.factorExplicit (X, Q, B, forceNonnegativeDiagonal_);
01546       // This call will only modify *B if *B on input is not of full
01547       // numerical rank.
01548       rank = tsqrAdaptor_.revealRank (Q, B, relativeRankTolerance_);
01549     } catch (std::exception& e) {
01550       throw TsqrOrthoError (e.what()); // Toss the exception up the chain.
01551     }
01552     return rank;
01553   }
01554 
01555   template<class Scalar, class MV>
01556   int
01557   TsqrOrthoManagerImpl<Scalar, MV>::
01558   normalizeOne (MV& X, 
01559     Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > B) const
01560   {
01561     // Make space for the normalization coefficient.  This will either
01562     // be a freshly allocated matrix (if B is null), or a view of the
01563     // 1x1 upper left submatrix of *B (if B is not null).
01564     mat_ptr B_out;
01565     if (B.is_null()) {
01566       B_out = rcp (new mat_type (1, 1));
01567     } else {
01568       const int numRows = B->numRows();
01569       const int numCols = B->numCols();
01570       TEUCHOS_TEST_FOR_EXCEPTION(numRows < 1 || numCols < 1, 
01571        std::invalid_argument,
01572        "normalizeOne: Input matrix B must be at "
01573        "least 1 x 1, but is instead " << numRows
01574        << " x " << numCols << ".");
01575       // Create a view of the 1x1 upper left submatrix of *B.
01576       B_out = rcp (new mat_type (Teuchos::View, *B, 1, 1));
01577     }
01578 
01579     // Compute the norm of X, and write the result to B_out.
01580     std::vector<magnitude_type> theNorm (1, SCTM::zero());
01581     MVT::MvNorm (X, theNorm);
01582     (*B_out)(0,0) = theNorm[0];
01583 
01584     if (B.is_null()) {
01585       // The input matrix B is null, so assign B_out to it.  If B was
01586       // not null on input, then B_out is a view of *B, so we don't
01587       // have to do anything here.  Note that SerialDenseMatrix uses
01588       // raw pointers to store data and represent views, so we have to
01589       // be careful about scope.
01590       B = B_out;
01591     }
01592 
01593     // Scale X by its norm, if its norm is zero.  Otherwise, do the
01594     // right thing based on whether the user wants us to fill the null
01595     // space with random vectors.
01596     if (theNorm[0] == SCTM::zero()) {
01597       // Make a view of the first column of Q, fill it with random
01598       // data, and normalize it.  Throw away the resulting norm.
01599       if (randomizeNullSpace_) {
01600   MVT::MvRandom(X);
01601   MVT::MvNorm (X, theNorm);
01602   if (theNorm[0] == SCTM::zero()) {
01603     // It is possible that a random vector could have all zero
01604     // entries, but unlikely.  We could try again, but it's also
01605     // possible that multiple tries could result in zero
01606     // vectors.  We choose instead to give up.
01607     throw TsqrOrthoError("normalizeOne: a supposedly random "
01608              "vector has norm zero!");
01609   } else {
01610     // NOTE (mfh 09 Nov 2010) I'm assuming that dividing a
01611     // Scalar by a magnitude_type is defined and that it results
01612     // in a Scalar.
01613     const Scalar alpha = SCT::one() / theNorm[0];
01614     MVT::MvScale (X, alpha);
01615   }
01616       }
01617       return 0; // The rank of the matrix (actually just one vector) X.
01618     } else {
01619       // NOTE (mfh 09 Nov 2010) I'm assuming that dividing a Scalar by
01620       // a magnitude_type is defined and that it results in a Scalar.
01621       const Scalar alpha = SCT::one() / theNorm[0];
01622       MVT::MvScale (X, alpha);
01623       return 1; // The rank of the matrix (actually just one vector) X.
01624     }
01625   }
01626 
01627 
01628   template<class Scalar, class MV>
01629   void
01630   TsqrOrthoManagerImpl<Scalar, MV>::
01631   rawProject (MV& X, 
01632         Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
01633         Teuchos::ArrayView<Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > > C) const
01634   {
01635 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01636     Teuchos::TimeMonitor timerMonitorNormalize(*timerProject_);
01637 #endif // BELOS_TEUCHOS_TIME_MONITOR
01638 
01639     // "Modified Gram-Schmidt" version of Block Gram-Schmidt.
01640     const int num_Q_blocks = Q.size();
01641     for (int i = 0; i < num_Q_blocks; ++i)
01642       {
01643   // TEUCHOS_TEST_FOR_EXCEPTION(C[i].is_null(), std::logic_error,
01644   //       "TsqrOrthoManagerImpl::rawProject(): C[" 
01645   //       << i << "] is null");
01646   // TEUCHOS_TEST_FOR_EXCEPTION(Q[i].is_null(), std::logic_error,
01647   //       "TsqrOrthoManagerImpl::rawProject(): Q[" 
01648   //       << i << "] is null");
01649   mat_type& Ci = *C[i];
01650   const MV& Qi = *Q[i];
01651   innerProd (Qi, X, Ci);
01652   MVT::MvTimesMatAddMv (-SCT::one(), Qi, Ci, SCT::one(), X);
01653       }
01654   }
01655 
01656 
01657   template<class Scalar, class MV>
01658   void
01659   TsqrOrthoManagerImpl<Scalar, MV>::
01660   rawProject (MV& X, 
01661         const Teuchos::RCP<const MV>& Q, 
01662         const Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> >& C) const
01663   {
01664 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01665     Teuchos::TimeMonitor timerMonitorNormalize(*timerProject_);
01666 #endif // BELOS_TEUCHOS_TIME_MONITOR
01667 
01668     // Block Gram-Schmidt
01669     innerProd (*Q, X, *C);
01670     MVT::MvTimesMatAddMv (-SCT::one(), *Q, *C, SCT::one(), X);
01671   }
01672 
01673   template<class Scalar, class MV>
01674   int
01675   TsqrOrthoManagerImpl<Scalar, MV>::
01676   normalizeImpl (MV& X, 
01677      MV& Q, 
01678      Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > B, 
01679      const bool outOfPlace)
01680   {
01681     using Teuchos::Range1D;
01682     using Teuchos::RCP;
01683     using Teuchos::rcp;
01684     using Teuchos::ScalarTraits;
01685     using Teuchos::tuple;
01686     using std::cerr;
01687     using std::endl;
01688     // Don't set this to true unless you want lots of debugging
01689     // messages written to stderr on every MPI process.
01690     const bool extraDebug = false;
01691 
01692     const int numCols = MVT::GetNumberVecs (X);
01693     if (numCols == 0) {
01694       return 0; // Fast exit for an empty input matrix.
01695     }
01696 
01697     // We allow Q to have more columns than X.  In that case, we only
01698     // touch the first numCols columns of Q.
01699     TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(Q) < numCols, 
01700            std::invalid_argument, 
01701            "TsqrOrthoManagerImpl::normalizeImpl(X,Q,B): Q has "
01702            << MVT::GetNumberVecs(Q) << " columns.  This is too "
01703            "few, since X has " << numCols << " columns.");
01704     // TSQR wants a Q with the same number of columns as X, so have it
01705     // work on a nonconstant view of Q with the same number of columns
01706     // as X.
01707     RCP<MV> Q_view = MVT::CloneViewNonConst (Q, Range1D(0, numCols-1));
01708 
01709     // Make space for the normalization coefficients.  This will
01710     // either be a freshly allocated matrix (if B is null), or a view
01711     // of the appropriately sized upper left submatrix of *B (if B is
01712     // not null).
01713     mat_ptr B_out;
01714     if (B.is_null()) {
01715       B_out = rcp (new mat_type (numCols, numCols));
01716     } else {
01717       // Make sure that B is no smaller than numCols x numCols.
01718       TEUCHOS_TEST_FOR_EXCEPTION(B->numRows() < numCols || B->numCols() < numCols,
01719        std::invalid_argument,
01720        "normalizeOne: Input matrix B must be at "
01721        "least " << numCols << " x " << numCols 
01722        << ", but is instead " << B->numRows()
01723        << " x " << B->numCols() << ".");
01724       // Create a view of the numCols x numCols upper left submatrix
01725       // of *B.  TSQR will write the normalization coefficients there.
01726       B_out = rcp (new mat_type (Teuchos::View, *B, numCols, numCols));
01727     }
01728 
01729     if (extraDebug) {
01730       std::vector<magnitude_type> norms (numCols);
01731       MVT::MvNorm (X, norms);
01732       cerr << "Column norms of X before orthogonalization: ";
01733       typedef typename std::vector<magnitude_type>::const_iterator iter_type;
01734       for (iter_type iter = norms.begin(); iter != norms.end(); ++iter) {
01735   cerr << *iter;
01736   if (iter+1 != norms.end())
01737     cerr << ", ";
01738       }
01739     }
01740 
01741     // Compute rank-revealing decomposition (in this case, TSQR of X
01742     // followed by SVD of the R factor and appropriate updating of the
01743     // resulting Q factor) of X.  X is modified in place and filled
01744     // with garbage, and Q_view contains the resulting explicit Q
01745     // factor.  Later, we will copy this back into X.
01746     //
01747     // The matrix *B_out will only be upper triangular if X is of full
01748     // numerical rank.  Otherwise, the entries below the diagonal may
01749     // be filled in as well.
01750     const int rank = rawNormalize (X, *Q_view, *B_out);
01751     if (B.is_null()) {
01752       // The input matrix B is null, so assign B_out to it.  If B was
01753       // not null on input, then B_out is a view of *B, so we don't
01754       // have to do anything here.  Note that SerialDenseMatrix uses
01755       // raw pointers to store data and represent views, so we have to
01756       // be careful about scope.
01757       B = B_out;
01758     }
01759 
01760     if (extraDebug) {
01761       std::vector<magnitude_type> norms (numCols);
01762       MVT::MvNorm (*Q_view, norms);
01763       cerr << "Column norms of Q_view after orthogonalization: ";
01764       for (typename std::vector<magnitude_type>::const_iterator iter = norms.begin(); 
01765      iter != norms.end(); ++iter) {
01766   cerr << *iter;
01767   if (iter+1 != norms.end())
01768     cerr << ", ";
01769       }
01770     }
01771     TEUCHOS_TEST_FOR_EXCEPTION(rank < 0 || rank > numCols, std::logic_error,
01772            "Belos::TsqrOrthoManagerImpl::normalizeImpl: "
01773            "rawNormalize() returned rank = " << rank << " for a "
01774            "matrix X with " << numCols << " columns.  Please report"
01775            " this bug to the Belos developers.");
01776     if (extraDebug && rank == 0) {
01777       // Sanity check: ensure that the columns of X are sufficiently
01778       // small for X to be reported as rank zero.
01779       const mat_type& B_ref = *B;
01780       std::vector<magnitude_type> norms (B_ref.numCols());
01781       for (typename mat_type::ordinalType j = 0; j < B_ref.numCols(); ++j) {
01782   typedef typename mat_type::scalarType mat_scalar_type;
01783   mat_scalar_type sumOfSquares = ScalarTraits<mat_scalar_type>::zero();
01784   for (typename mat_type::ordinalType i = 0; i <= j; ++i) {
01785     const mat_scalar_type B_ij = 
01786       ScalarTraits<mat_scalar_type>::magnitude (B_ref(i,j));
01787     sumOfSquares += B_ij*B_ij;
01788   }
01789   norms[j] = ScalarTraits<mat_scalar_type>::squareroot (sumOfSquares);
01790       }
01791       using std::cerr;
01792       using std::endl;
01793       cerr << "Norms of columns of B after orthogonalization: ";
01794       for (typename mat_type::ordinalType j = 0; j < B_ref.numCols(); ++j) {
01795   cerr << norms[j];
01796   if (j != B_ref.numCols() - 1)
01797     cerr << ", ";
01798       }
01799       cerr << endl;
01800     }
01801 
01802     // If X is full rank or we don't want to replace its null space
01803     // basis with random vectors, then we're done.
01804     if (rank == numCols || ! randomizeNullSpace_) {
01805       // If we're supposed to be working in place in X, copy the
01806       // results back from Q_view into X.
01807       if (! outOfPlace) {
01808   MVT::Assign (*Q_view, X);
01809       }
01810       return rank;
01811     }
01812 
01813     if (randomizeNullSpace_ && rank < numCols) {
01814       // X wasn't full rank.  Replace the null space basis of X (in
01815       // the last numCols-rank columns of Q_view) with random data,
01816       // project it against the first rank columns of Q_view, and
01817       // normalize.
01818       //
01819       // Number of columns to fill with random data.
01820       const int nullSpaceNumCols = numCols - rank;
01821       // Inclusive range of indices of columns of X to fill with
01822       // random data.
01823       Range1D nullSpaceIndices (rank, numCols-1);
01824 
01825       // rawNormalize() wrote the null space basis vectors into
01826       // Q_view.  We continue to work in place in Q_view by writing
01827       // the random data there and (if there is a nontrival column
01828       // space) projecting in place against the column space basis
01829       // vectors (also in Q_view).
01830       RCP<MV> Q_null = MVT::CloneViewNonConst (*Q_view, nullSpaceIndices);
01831       // Replace the null space basis with random data.
01832       MVT::MvRandom (*Q_null); 
01833 
01834       // Make sure that the "random" data isn't all zeros.  This is
01835       // statistically nearly impossible, but we test for debugging
01836       // purposes.
01837       {
01838   std::vector<magnitude_type> norms (MVT::GetNumberVecs(*Q_null));
01839   MVT::MvNorm (*Q_null, norms);
01840 
01841   bool anyZero = false;
01842   typedef typename std::vector<magnitude_type>::const_iterator iter_type;
01843   for (iter_type it = norms.begin(); it != norms.end(); ++it) {
01844     if (*it == SCTM::zero()) {
01845       anyZero = true;
01846     }
01847   }
01848   if (anyZero) {
01849     std::ostringstream os;
01850     os << "TsqrOrthoManagerImpl::normalizeImpl: "
01851       "We are being asked to randomize the null space, for a matrix "
01852       "with " << numCols << " columns and reported column rank "
01853        << rank << ".  The inclusive range of columns to fill with "
01854       "random data is [" << nullSpaceIndices.lbound() << "," 
01855        << nullSpaceIndices.ubound() << "].  After filling the null "
01856       "space vectors with random numbers, at least one of the vectors"
01857       " has norm zero.  Here are the norms of all the null space "
01858       "vectors: [";
01859     for (iter_type it = norms.begin(); it != norms.end(); ++it) {
01860       os << *it;
01861       if (it+1 != norms.end())
01862         os << ", ";
01863     }
01864     os << "].)  There is a tiny probability that this could happen "
01865       "randomly, but it is likely a bug.  Please report it to the "
01866       "Belos developers, especially if you are able to reproduce the "
01867       "behavior.";
01868     TEUCHOS_TEST_FOR_EXCEPTION(anyZero, TsqrOrthoError, os.str());
01869   }
01870       }
01871 
01872       if (rank > 0) {
01873   // Project the random data against the column space basis of
01874   // X, using a simple block projection ("Block Classical
01875   // Gram-Schmidt").  This is accurate because we've already
01876   // orthogonalized the column space basis of X nearly to
01877   // machine precision via a QR factorization (TSQR) with
01878   // accuracy comparable to Householder QR.
01879   RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D(0, rank-1));
01880 
01881   // Temporary storage for projection coefficients.  We don't
01882   // need to keep them, since they represent the null space
01883   // basis (for which the coefficients are logically zero).
01884   mat_ptr C_null (new mat_type (rank, nullSpaceNumCols));
01885   rawProject (*Q_null, Q_col, C_null);
01886       }
01887       // Normalize the projected random vectors, so that they are
01888       // mutually orthogonal (as well as orthogonal to the column
01889       // space basis of X).  We use X for the output of the
01890       // normalization: for out-of-place normalization (outOfPlace ==
01891       // true), X is overwritten with "invalid values" anyway, and for
01892       // in-place normalization (outOfPlace == false), we want the
01893       // result to be in X anyway.
01894       RCP<MV> X_null = MVT::CloneViewNonConst (X, nullSpaceIndices);
01895       // Normalization coefficients for projected random vectors.
01896       // Will be thrown away.
01897       mat_type B_null (nullSpaceNumCols, nullSpaceNumCols);
01898       // Write the normalized vectors to X_null (in X).
01899       const int nullSpaceBasisRank = rawNormalize (*Q_null, *X_null, B_null);
01900     
01901       // It's possible, but unlikely, that X_null doesn't have full
01902       // rank (after the projection step).  We could recursively fill
01903       // in more random vectors until we finally get a full rank
01904       // matrix, but instead we just throw an exception.
01905       //
01906       // NOTE (mfh 08 Nov 2010) Perhaps we should deal with this case
01907       // more elegantly.  Recursion might be one way to solve it, but
01908       // be sure to check that the recursion will terminate.  We could
01909       // do this by supplying an additional argument to rawNormalize,
01910       // which is the null space basis rank from the previous
01911       // iteration.  The rank has to decrease each time, or the
01912       // recursion may go on forever.
01913       if (nullSpaceBasisRank < nullSpaceNumCols) {
01914   std::vector<magnitude_type> norms (MVT::GetNumberVecs(*X_null));
01915   MVT::MvNorm (*X_null, norms);
01916   std::ostringstream os;
01917   os << "TsqrOrthoManagerImpl::normalizeImpl: "
01918      << "We are being asked to randomize the null space, "
01919      << "for a matrix with " << numCols << " columns and "
01920      << "column rank " << rank << ".  After projecting and "
01921      << "normalizing the generated random vectors, they "
01922      << "only have rank " << nullSpaceBasisRank << ".  They"
01923      << " should have full rank " << nullSpaceNumCols 
01924      << ".  (The inclusive range of columns to fill with "
01925      << "random data is [" << nullSpaceIndices.lbound() 
01926      << "," << nullSpaceIndices.ubound() << "].  The "
01927      << "column norms of the resulting Q factor are: [";
01928   for (typename std::vector<magnitude_type>::size_type k = 0; 
01929        k < norms.size(); ++k) {
01930     os << norms[k];
01931     if (k != norms.size()-1) {
01932       os << ", ";
01933     }
01934   }
01935   os << "].)  There is a tiny probability that this could "
01936      << "happen randomly, but it is likely a bug.  Please "
01937      << "report it to the Belos developers, especially if "
01938      << "you are able to reproduce the behavior.";
01939 
01940   TEUCHOS_TEST_FOR_EXCEPTION(nullSpaceBasisRank < nullSpaceNumCols, 
01941          TsqrOrthoError, os.str());
01942       }
01943       // If we're normalizing out of place, copy the X_null
01944       // vectors back into Q_null; the Q_col vectors are already
01945       // where they are supposed to be in that case.
01946       //
01947       // If we're normalizing in place, leave X_null alone (it's
01948       // already where it needs to be, in X), but copy Q_col back
01949       // into the first rank columns of X.
01950       if (outOfPlace) {
01951   MVT::Assign (*X_null, *Q_null);
01952       } else if (rank > 0) {
01953   // MVT::Assign() doesn't accept empty ranges of columns.
01954   RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D(0, rank-1));
01955   RCP<MV> X_col = MVT::CloneViewNonConst (X, Range1D(0, rank-1));
01956   MVT::Assign (*Q_col, *X_col);
01957       }
01958     }
01959     return rank;
01960   }
01961 
01962 
01963   template<class Scalar, class MV>
01964   void
01965   TsqrOrthoManagerImpl<Scalar, MV>::
01966   checkProjectionDims (int& ncols_X, 
01967            int& num_Q_blocks,
01968            int& ncols_Q_total,
01969            const MV& X, 
01970            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
01971   {
01972     // First assign to temporary values, so the function won't
01973     // commit any externally visible side effects unless it will
01974     // return normally (without throwing an exception).  (I'm being
01975     // cautious; MVT::GetNumberVecs() probably shouldn't have any
01976     // externally visible side effects, unless it is logging to a
01977     // file or something.)
01978     int the_ncols_X, the_num_Q_blocks, the_ncols_Q_total;
01979     the_num_Q_blocks = Q.size();
01980     the_ncols_X = MVT::GetNumberVecs (X);
01981 
01982     // Compute the total number of columns of all the Q[i] blocks.
01983     the_ncols_Q_total = 0;
01984     // You should be angry if your compiler doesn't support type
01985     // inference ("auto").  That's why I need this awful typedef.
01986     using Teuchos::ArrayView;
01987     using Teuchos::RCP;
01988     typedef typename ArrayView<RCP<const MV> >::const_iterator iter_type;
01989     for (iter_type it = Q.begin(); it != Q.end(); ++it) {
01990       const MV& Qi = **it;
01991       the_ncols_Q_total += MVT::GetNumberVecs (Qi);
01992     }
01993 
01994     // Commit temporary values to the output arguments.
01995     ncols_X = the_ncols_X;
01996     num_Q_blocks = the_num_Q_blocks;
01997     ncols_Q_total = the_ncols_Q_total;
01998   }
01999 
02000 } // namespace Belos
02001 
02002 #endif // __BelosTsqrOrthoManagerImpl_hpp
02003 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines