Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosTsqrOrthoManagerImpl.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
45 #ifndef __BelosTsqrOrthoManagerImpl_hpp
46 #define __BelosTsqrOrthoManagerImpl_hpp
47 
48 #include "BelosConfigDefs.hpp" // HAVE_BELOS_TSQR
49 #include "BelosMultiVecTraits.hpp"
50 #include "BelosOrthoManager.hpp" // OrthoError, etc.
51 
52 #include "Teuchos_as.hpp"
53 #include "Teuchos_LAPACK.hpp"
55 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
56 #ifdef BELOS_TEUCHOS_TIME_MONITOR
57 # include "Teuchos_TimeMonitor.hpp"
58 #endif // BELOS_TEUCHOS_TIME_MONITOR
59 #include <algorithm>
60 #include <functional> // std::binary_function
61 
62 namespace Belos {
63 
67  class TsqrOrthoError : public OrthoError {
68  public:
69  TsqrOrthoError (const std::string& what_arg) :
70  OrthoError (what_arg) {}
71  };
72 
92  class TsqrOrthoFault : public OrthoError {
93  public:
94  TsqrOrthoFault (const std::string& what_arg) :
95  OrthoError (what_arg) {}
96  };
97 
130  template<class Scalar>
132  public std::binary_function<Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>,
133  Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>,
134  void>
135  {
136  public:
138  typedef Scalar scalar_type;
144 
146  virtual ~ReorthogonalizationCallback ();
147 
152  virtual void
153  operator() (Teuchos::ArrayView<magnitude_type> normsBeforeFirstPass,
154  Teuchos::ArrayView<magnitude_type> normsAfterFirstPass) = 0;
155  };
156 
157  template<class Scalar>
159 
191  template<class Scalar, class MV>
194  public:
195  typedef Scalar scalar_type;
197  typedef MV multivector_type;
201 
202  private:
206  typedef typename MVT::tsqr_adaptor_type tsqr_adaptor_type;
207 
208  public:
217 
220 
232 
250  const std::string& label);
251 
256  TsqrOrthoManagerImpl (const std::string& label);
257 
277  void
279  {
280  reorthogCallback_ = callback;
281  }
282 
290  void setLabel (const std::string& label) {
291  if (label != label_) {
292  label_ = label;
293 
294 #ifdef BELOS_TEUCHOS_TIME_MONITOR
295  clearTimer (label, "All orthogonalization");
296  clearTimer (label, "Projection");
297  clearTimer (label, "Normalization");
298 
299  timerOrtho_ = makeTimer (label, "All orthogonalization");
300  timerProject_ = makeTimer (label, "Projection");
301  timerNormalize_ = makeTimer (label, "Normalization");
302 #endif // BELOS_TEUCHOS_TIME_MONITOR
303  }
304  }
305 
307  const std::string& getLabel () const { return label_; }
308 
317  void
318  innerProd (const MV& X, const MV& Y, mat_type& Z) const
319  {
320  MVT::MvTransMv (SCT::one(), X, Y, Z);
321  }
322 
340  void
341  norm (const MV& X, std::vector<magnitude_type>& normVec) const;
342 
352  void
353  project (MV& X,
356 
370  int normalize (MV& X, mat_ptr B);
371 
390  int
391  normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B);
392 
406  int
409  mat_ptr B,
411  {
412  // "false" means we work on X in place. The second argument is
413  // not read or written in that case.
414  return projectAndNormalizeImpl (X, X, false, C, B, Q);
415  }
416 
436  int
438  MV& X_out,
440  mat_ptr B,
442  {
443  // "true" means we work on X_in out of place, writing the
444  // results into X_out.
445  return projectAndNormalizeImpl (X_in, X_out, true, C, B, Q);
446  }
447 
453  orthonormError (const MV &X) const
454  {
455  const Scalar ONE = SCT::one();
456  const int ncols = MVT::GetNumberVecs(X);
457  mat_type XTX (ncols, ncols);
458  innerProd (X, X, XTX);
459  for (int k = 0; k < ncols; ++k) {
460  XTX(k,k) -= ONE;
461  }
462  return XTX.normFrobenius();
463  }
464 
467  orthogError (const MV &X1,
468  const MV &X2) const
469  {
470  const int ncols_X1 = MVT::GetNumberVecs (X1);
471  const int ncols_X2 = MVT::GetNumberVecs (X2);
472  mat_type X1_T_X2 (ncols_X1, ncols_X2);
473  innerProd (X1, X2, X1_T_X2);
474  return X1_T_X2.normFrobenius();
475  }
476 
480  magnitude_type blockReorthogThreshold() const { return blockReorthogThreshold_; }
481 
484  magnitude_type relativeRankTolerance() const { return relativeRankTolerance_; }
485 
486  private:
489 
491  mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
492 
494  std::string label_;
495 
497  tsqr_adaptor_type tsqrAdaptor_;
498 
508  Teuchos::RCP<MV> Q_;
509 
511  magnitude_type eps_;
512 
516  bool randomizeNullSpace_;
517 
523  bool reorthogonalizeBlocks_;
524 
528  bool throwOnReorthogFault_;
529 
531  magnitude_type blockReorthogThreshold_;
532 
534  magnitude_type relativeRankTolerance_;
535 
542  bool forceNonnegativeDiagonal_;
543 
544 #ifdef BELOS_TEUCHOS_TIME_MONITOR
545  Teuchos::RCP<Teuchos::Time> timerOrtho_;
547 
549  Teuchos::RCP<Teuchos::Time> timerProject_;
550 
552  Teuchos::RCP<Teuchos::Time> timerNormalize_;
553 #endif // BELOS_TEUCHOS_TIME_MONITOR
554 
557 
558 #ifdef BELOS_TEUCHOS_TIME_MONITOR
568  makeTimer (const std::string& prefix,
569  const std::string& timerName)
570  {
571  const std::string timerLabel =
572  prefix.empty() ? timerName : (prefix + ": " + timerName);
573  return Teuchos::TimeMonitor::getNewCounter (timerLabel);
574  }
575 
581  void
582  clearTimer (const std::string& prefix,
583  const std::string& timerName)
584  {
585  const std::string timerLabel =
586  prefix.empty() ? timerName : (prefix + ": " + timerName);
588  }
589 #endif // BELOS_TEUCHOS_TIME_MONITOR
590 
592  void
593  raiseReorthogFault (const std::vector<magnitude_type>& normsAfterFirstPass,
594  const std::vector<magnitude_type>& normsAfterSecondPass,
595  const std::vector<int>& faultIndices);
596 
606  void
607  checkProjectionDims (int& ncols_X,
608  int& num_Q_blocks,
609  int& ncols_Q_total,
610  const MV& X,
612 
623  void
624  allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C,
626  const MV& X,
627  const bool attemptToRecycle = true) const;
628 
637  int
638  projectAndNormalizeImpl (MV& X_in,
639  MV& X_out,
640  const bool outOfPlace,
642  mat_ptr B,
644 
650  void
651  rawProject (MV& X,
654 
656  void
657  rawProject (MV& X,
658  const Teuchos::RCP<const MV>& Q,
659  const mat_ptr& C) const;
660 
687  int rawNormalize (MV& X, MV& Q, mat_type& B);
688 
706  int normalizeOne (MV& X, mat_ptr B) const;
707 
735  int normalizeImpl (MV& X, MV& Q, mat_ptr B, const bool outOfPlace);
736  };
737 
738  template<class Scalar, class MV>
739  void
742  {
744  using Teuchos::parameterList;
745  using Teuchos::RCP;
746  using Teuchos::sublist;
747  typedef magnitude_type M; // abbreviation.
748 
749  RCP<const ParameterList> defaultParams = getValidParameters ();
750  // Sublist of TSQR implementation parameters; to get below.
751  RCP<ParameterList> tsqrParams;
752 
753  RCP<ParameterList> theParams;
754  if (params.is_null()) {
755  theParams = parameterList (*defaultParams);
756  } else {
757  theParams = params;
758 
759  // Don't call validateParametersAndSetDefaults(); we prefer to
760  // ignore parameters that we don't recognize, at least for now.
761  // However, we do fill in missing parameters with defaults.
762 
763  randomizeNullSpace_ =
764  theParams->get<bool> ("randomizeNullSpace",
765  defaultParams->get<bool> ("randomizeNullSpace"));
766  reorthogonalizeBlocks_ =
767  theParams->get<bool> ("reorthogonalizeBlocks",
768  defaultParams->get<bool> ("reorthogonalizeBlocks"));
769  throwOnReorthogFault_ =
770  theParams->get<bool> ("throwOnReorthogFault",
771  defaultParams->get<bool> ("throwOnReorthogFault"));
772  blockReorthogThreshold_ =
773  theParams->get<M> ("blockReorthogThreshold",
774  defaultParams->get<M> ("blockReorthogThreshold"));
775  relativeRankTolerance_ =
776  theParams->get<M> ("relativeRankTolerance",
777  defaultParams->get<M> ("relativeRankTolerance"));
778  forceNonnegativeDiagonal_ =
779  theParams->get<bool> ("forceNonnegativeDiagonal",
780  defaultParams->get<bool> ("forceNonnegativeDiagonal"));
781 
782  // Get the sublist of TSQR implementation parameters. Use the
783  // default sublist if one isn't provided.
784  if (! theParams->isSublist ("TSQR implementation")) {
785  theParams->set ("TSQR implementation",
786  defaultParams->sublist ("TSQR implementation"));
787  }
788  tsqrParams = sublist (theParams, "TSQR implementation", true);
789  }
790 
791  // Send the TSQR implementation parameters to the TSQR adaptor.
792  tsqrAdaptor_.setParameterList (tsqrParams);
793 
794  // Save the input parameter list.
795  setMyParamList (theParams);
796  }
797 
798  template<class Scalar, class MV>
801  const std::string& label) :
802  label_ (label),
803  Q_ (Teuchos::null), // Initialized on demand
804  eps_ (SCTM::eps()), // Machine precision
805  randomizeNullSpace_ (true),
806  reorthogonalizeBlocks_ (true),
807  throwOnReorthogFault_ (false),
808  blockReorthogThreshold_ (0),
809  relativeRankTolerance_ (0),
810  forceNonnegativeDiagonal_ (false)
811  {
812  setParameterList (params); // This also sets tsqrAdaptor_'s parameters.
813 
814 #ifdef BELOS_TEUCHOS_TIME_MONITOR
815  timerOrtho_ = makeTimer (label, "All orthogonalization");
816  timerProject_ = makeTimer (label, "Projection");
817  timerNormalize_ = makeTimer (label, "Normalization");
818 #endif // BELOS_TEUCHOS_TIME_MONITOR
819  }
820 
821  template<class Scalar, class MV>
823  TsqrOrthoManagerImpl (const std::string& label) :
824  label_ (label),
825  Q_ (Teuchos::null), // Initialized on demand
826  eps_ (SCTM::eps()), // Machine precision
827  randomizeNullSpace_ (true),
828  reorthogonalizeBlocks_ (true),
829  throwOnReorthogFault_ (false),
830  blockReorthogThreshold_ (0),
831  relativeRankTolerance_ (0),
832  forceNonnegativeDiagonal_ (false)
833  {
834  setParameterList (Teuchos::null); // Set default parameters.
835 
836 #ifdef BELOS_TEUCHOS_TIME_MONITOR
837  timerOrtho_ = makeTimer (label, "All orthogonalization");
838  timerProject_ = makeTimer (label, "Projection");
839  timerNormalize_ = makeTimer (label, "Normalization");
840 #endif // BELOS_TEUCHOS_TIME_MONITOR
841  }
842 
843  template<class Scalar, class MV>
844  void
846  norm (const MV& X, std::vector<magnitude_type>& normVec) const
847  {
848  const int numCols = MVT::GetNumberVecs (X);
849  // std::vector<T>::size_type is unsigned; int is signed. Mixed
850  // unsigned/signed comparisons trigger compiler warnings.
851  if (normVec.size() < static_cast<size_t>(numCols))
852  normVec.resize (numCols); // Resize normvec if necessary.
853  MVT::MvNorm (X, normVec);
854  }
855 
856  template<class Scalar, class MV>
857  void
861  {
862 #ifdef BELOS_TEUCHOS_TIME_MONITOR
863  // "Projection" only happens in rawProject(), so we only time
864  // projection inside rawProject(). However, we count the time
865  // spend in project() as part of the whole orthogonalization.
866  //
867  // If project() is called from projectAndNormalize(), the
868  // TimeMonitor won't start timerOrtho_, because it is already
869  // running in projectAndNormalize().
870  Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
871 #endif // BELOS_TEUCHOS_TIME_MONITOR
872 
873  int ncols_X, num_Q_blocks, ncols_Q_total;
874  checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X, Q);
875  // Test for quick exit: any dimension of X is zero, or there are
876  // zero Q blocks, or the total number of columns of the Q blocks
877  // is zero.
878  if (ncols_X == 0 || num_Q_blocks == 0 || ncols_Q_total == 0)
879  return;
880 
881  // Make space for first-pass projection coefficients
882  allocateProjectionCoefficients (C, Q, X, true);
883 
884  // We only use columnNormsBefore and compute pre-projection column
885  // norms if doing block reorthogonalization.
886  std::vector<magnitude_type> columnNormsBefore (ncols_X, magnitude_type(0));
887  if (reorthogonalizeBlocks_)
888  MVT::MvNorm (X, columnNormsBefore);
889 
890  // Project (first block orthogonalization step):
891  // C := Q^* X, X := X - Q C.
892  rawProject (X, Q, C);
893 
894  // If we are doing block reorthogonalization, reorthogonalize X if
895  // necessary.
896  if (reorthogonalizeBlocks_) {
897  std::vector<magnitude_type> columnNormsAfter (ncols_X, magnitude_type(0));
898  MVT::MvNorm (X, columnNormsAfter);
899 
900  // Relative block reorthogonalization threshold.
901  const magnitude_type relThres = blockReorthogThreshold();
902  // Reorthogonalize X if any of its column norms decreased by a
903  // factor more than the block reorthogonalization threshold.
904  // Don't bother trying to subset the columns; that will make the
905  // columns noncontiguous and thus hinder BLAS 3 optimizations.
906  bool reorthogonalize = false;
907  for (int j = 0; j < ncols_X; ++j) {
908  if (columnNormsAfter[j] < relThres * columnNormsBefore[j]) {
909  reorthogonalize = true;
910  break;
911  }
912  }
913  if (reorthogonalize) {
914  // Notify the caller via callback about the need for
915  // reorthogonalization.
916  if (! reorthogCallback_.is_null()) {
917  reorthogCallback_->operator() (Teuchos::arrayViewFromVector (columnNormsBefore),
918  Teuchos::arrayViewFromVector (columnNormsAfter));
919  }
920  // Second-pass projection coefficients
922  allocateProjectionCoefficients (C2, Q, X, false);
923 
924  // Perform the second projection pass:
925  // C2 = Q' X, X = X - Q*C2
926  rawProject (X, Q, C2);
927  // Update the projection coefficients
928  for (int k = 0; k < num_Q_blocks; ++k)
929  *C[k] += *C2[k];
930  }
931  }
932  }
933 
934 
935  template<class Scalar, class MV>
936  int
938  {
939  using Teuchos::Range1D;
940  using Teuchos::RCP;
941 
942 #ifdef BELOS_TEUCHOS_TIME_MONITOR
943  Teuchos::TimeMonitor timerMonitorNormalize(*timerNormalize_);
944  // If normalize() is called internally -- i.e., called from
945  // projectAndNormalize() -- the timer will not be started or
946  // stopped, because it is already running. TimeMonitor handles
947  // recursive invocation by doing nothing.
948  Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
949 #endif // BELOS_TEUCHOS_TIME_MONITOR
950 
951  // MVT returns int for this, even though the "local ordinal
952  // type" of the MV may be some other type (for example,
953  // Tpetra::MultiVector<double, int32_t, int64_t, ...>).
954  const int numCols = MVT::GetNumberVecs (X);
955 
956  // This special case (for X having only one column) makes
957  // TsqrOrthoManagerImpl equivalent to Modified Gram-Schmidt
958  // orthogonalization with conditional full reorthogonalization,
959  // if all multivector inputs have only one column. It also
960  // avoids allocating Q_ scratch space and copying data when we
961  // don't need to invoke TSQR at all.
962  if (numCols == 1) {
963  return normalizeOne (X, B);
964  }
965 
966  // We use Q_ as scratch space for the normalization, since TSQR
967  // requires a scratch multivector (it can't factor in place). Q_
968  // should come from a vector space compatible with X's vector
969  // space, and Q_ should have at least as many columns as X.
970  // Otherwise, we have to reallocate. We also have to allocate
971  // (not "re-") Q_ if we haven't allocated it before. (We can't
972  // allocate Q_ until we have some X, so we need a multivector as
973  // the "prototype.")
974  //
975  // NOTE (mfh 11 Jan 2011) We only increase the number of columsn
976  // in Q_, never decrease. This is OK for typical uses of TSQR,
977  // but you might prefer different behavior in some cases.
978  //
979  // NOTE (mfh 10 Mar 2011) We should only reuse the scratch space
980  // Q_ if X and Q_ have compatible data distributions. However,
981  // Belos' current MultiVecTraits interface does not let us check
982  // for this. Thus, we can only check whether X and Q_ have the
983  // same number of rows. This will behave correctly for the common
984  // case in Belos that all multivectors with the same number of
985  // rows have the same data distribution.
986  //
987  // The specific MV implementation may do more checks than this on
988  // its own and throw an exception if X and Q_ are not compatible,
989  // but it may not. If you find that recycling the Q_ space causes
990  // troubles, you may consider modifying the code below to
991  // reallocate Q_ for every X that comes in.
992  if (Q_.is_null() ||
993  MVT::GetGlobalLength(*Q_) != MVT::GetGlobalLength(X) ||
994  numCols > MVT::GetNumberVecs (*Q_)) {
995  Q_ = MVT::Clone (X, numCols);
996  }
997 
998  // normalizeImpl() wants the second MV argument to have the same
999  // number of columns as X. To ensure this, we pass it a view of
1000  // Q_ if Q_ has more columns than X. (This is possible if we
1001  // previously called normalize() with a different multivector,
1002  // since we never reallocate Q_ if it has more columns than
1003  // necessary.)
1004  if (MVT::GetNumberVecs(*Q_) == numCols) {
1005  return normalizeImpl (X, *Q_, B, false);
1006  } else {
1007  RCP<MV> Q_view = MVT::CloneViewNonConst (*Q_, Range1D(0, numCols-1));
1008  return normalizeImpl (X, *Q_view, B, false);
1009  }
1010  }
1011 
1012  template<class Scalar, class MV>
1013  void
1017  const MV& X,
1018  const bool attemptToRecycle) const
1019  {
1020  const int num_Q_blocks = Q.size();
1021  const int ncols_X = MVT::GetNumberVecs (X);
1022  C.resize (num_Q_blocks);
1023  if (attemptToRecycle)
1024  {
1025  for (int i = 0; i < num_Q_blocks; ++i)
1026  {
1027  const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
1028  // Create a new C[i] if necessary, otherwise resize if
1029  // necessary, otherwise fill with zeros.
1030  if (C[i].is_null())
1031  C[i] = Teuchos::rcp (new mat_type (ncols_Qi, ncols_X));
1032  else
1033  {
1034  mat_type& Ci = *C[i];
1035  if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X)
1036  Ci.shape (ncols_Qi, ncols_X);
1037  else
1038  Ci.putScalar (SCT::zero());
1039  }
1040  }
1041  }
1042  else
1043  {
1044  for (int i = 0; i < num_Q_blocks; ++i)
1045  {
1046  const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
1047  C[i] = Teuchos::rcp (new mat_type (ncols_Qi, ncols_X));
1048  }
1049  }
1050  }
1051 
1052  template<class Scalar, class MV>
1053  int
1056  {
1057 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1058  Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
1059  Teuchos::TimeMonitor timerMonitorNormalize(*timerNormalize_);
1060 #endif // BELOS_TEUCHOS_TIME_MONITOR
1061 
1062  const int numVecs = MVT::GetNumberVecs(X);
1063  if (numVecs == 0) {
1064  return 0; // Nothing to do.
1065  } else if (numVecs == 1) {
1066  // Special case for a single column; scale and copy over.
1067  using Teuchos::Range1D;
1068  using Teuchos::RCP;
1069  using Teuchos::rcp;
1070 
1071  // Normalize X in place (faster than TSQR for one column).
1072  const int rank = normalizeOne (X, B);
1073  // Copy results to first column of Q.
1074  RCP<MV> Q_0 = MVT::CloneViewNonConst (Q, Range1D(0,0));
1075  MVT::Assign (X, *Q_0);
1076  return rank;
1077  } else {
1078  // The "true" argument to normalizeImpl() means the output
1079  // vectors go into Q, and the contents of X are overwritten with
1080  // invalid values.
1081  return normalizeImpl (X, Q, B, true);
1082  }
1083  }
1084 
1085  template<class Scalar, class MV>
1086  int
1088  projectAndNormalizeImpl (MV& X_in,
1089  MV& X_out, // Only written if outOfPlace==false.
1090  const bool outOfPlace,
1092  mat_ptr B,
1094  {
1095  using Teuchos::Range1D;
1096  using Teuchos::RCP;
1097  using Teuchos::rcp;
1098 
1099 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1100  // Projection is only timed in rawProject(), and normalization is
1101  // only timed in normalize() and normalizeOutOfPlace().
1102  Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
1103 #endif // BELOS_TEUCHOS_TIME_MONITOR
1104 
1105  if (outOfPlace) {
1106  // Make sure that X_out has at least as many columns as X_in.
1107  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(X_out) < MVT::GetNumberVecs(X_in),
1108  std::invalid_argument,
1109  "Belos::TsqrOrthoManagerImpl::"
1110  "projectAndNormalizeImpl(..., outOfPlace=true, ...):"
1111  "X_out has " << MVT::GetNumberVecs(X_out)
1112  << " columns, but X_in has "
1113  << MVT::GetNumberVecs(X_in) << " columns.");
1114  }
1115  // Fetch dimensions of X_in and Q, and allocate space for first-
1116  // and second-pass projection coefficients (C resp. C2).
1117  int ncols_X, num_Q_blocks, ncols_Q_total;
1118  checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X_in, Q);
1119 
1120  // Test for quick exit: if any dimension of X is zero.
1121  if (ncols_X == 0) {
1122  return 0;
1123  }
1124  // If there are zero Q blocks or zero Q columns, just normalize!
1125  if (num_Q_blocks == 0 || ncols_Q_total == 0) {
1126  if (outOfPlace) {
1127  return normalizeOutOfPlace (X_in, X_out, B);
1128  } else {
1129  return normalize (X_in, B);
1130  }
1131  }
1132 
1133  // The typical case is that the entries of C have been allocated
1134  // before, so we attempt to recycle the allocations. The call
1135  // below will reallocate if it cannot recycle.
1136  allocateProjectionCoefficients (C, Q, X_in, true);
1137 
1138  // If we are doing block reorthogonalization, then compute the
1139  // column norms of X before projecting for the first time. This
1140  // will help us decide whether we need to reorthogonalize X.
1141  std::vector<magnitude_type> normsBeforeFirstPass (ncols_X, SCTM::zero());
1142  if (reorthogonalizeBlocks_) {
1143  MVT::MvNorm (X_in, normsBeforeFirstPass);
1144  }
1145 
1146  // First (Modified) Block Gram-Schmidt pass, in place in X_in.
1147  rawProject (X_in, Q, C);
1148 
1149  // Make space for the normalization coefficients. This will
1150  // either be a freshly allocated matrix (if B is null), or a view
1151  // of the appropriately sized upper left submatrix of *B (if B is
1152  // not null).
1153  //
1154  // Note that if we let the normalize() routine allocate (in the
1155  // case that B is null), that storage will go away at the end of
1156  // normalize(). (This is because it passes the RCP by value, not
1157  // by reference.)
1158  mat_ptr B_out;
1159  if (B.is_null()) {
1160  B_out = rcp (new mat_type (ncols_X, ncols_X));
1161  } else {
1162  // Make sure that B is no smaller than numCols x numCols.
1163  TEUCHOS_TEST_FOR_EXCEPTION(B->numRows() < ncols_X || B->numCols() < ncols_X,
1164  std::invalid_argument,
1165  "normalizeOne: Input matrix B must be at "
1166  "least " << ncols_X << " x " << ncols_X
1167  << ", but is instead " << B->numRows()
1168  << " x " << B->numCols() << ".");
1169  // Create a view of the ncols_X by ncols_X upper left
1170  // submatrix of *B. TSQR will write the normalization
1171  // coefficients there.
1172  B_out = rcp (new mat_type (Teuchos::View, *B, ncols_X, ncols_X));
1173  }
1174 
1175  // Rank of X(_in) after first projection pass. If outOfPlace,
1176  // this overwrites X_in with invalid values, and the results go in
1177  // X_out. Otherwise, it's done in place in X_in.
1178  const int firstPassRank = outOfPlace ?
1179  normalizeOutOfPlace (X_in, X_out, B_out) :
1180  normalize (X_in, B_out);
1181  if (B.is_null()) {
1182  // The input matrix B is null, so assign B_out to it. If B was
1183  // not null on input, then B_out is a view of *B, so we don't
1184  // have to do anything here. Note that SerialDenseMatrix uses
1185  // raw pointers to store data and represent views, so we have to
1186  // be careful about scope.
1187  B = B_out;
1188  }
1189  int rank = firstPassRank; // Current rank of X.
1190 
1191  // If X was not full rank after projection and randomizeNullSpace_
1192  // is true, then normalize(OutOfPlace)() replaced the null space
1193  // basis of X with random vectors, and orthogonalized them against
1194  // the column space basis of X. However, we still need to
1195  // orthogonalize the random vectors against the Q[i], after which
1196  // we will need to renormalize them.
1197  //
1198  // If outOfPlace, then we need to work in X_out (where
1199  // normalizeOutOfPlace() wrote the normalized vectors).
1200  // Otherwise, we need to work in X_in.
1201  //
1202  // Note: We don't need to keep the new projection coefficients,
1203  // since they are multiplied by the "small" part of B
1204  // corresponding to the null space of the original X.
1205  if (firstPassRank < ncols_X && randomizeNullSpace_) {
1206  const int numNullSpaceCols = ncols_X - firstPassRank;
1207  const Range1D nullSpaceIndices (firstPassRank, ncols_X - 1);
1208 
1209  // Space for projection coefficients (will be thrown away)
1210  Teuchos::Array<mat_ptr> C_null (num_Q_blocks);
1211  for (int k = 0; k < num_Q_blocks; ++k) {
1212  const int numColsQk = MVT::GetNumberVecs(*Q[k]);
1213  C_null[k] = rcp (new mat_type (numColsQk, numNullSpaceCols));
1214  }
1215  // Space for normalization coefficients (will be thrown away).
1216  RCP<mat_type> B_null (new mat_type (numNullSpaceCols, numNullSpaceCols));
1217 
1218  int randomVectorsRank;
1219  if (outOfPlace) {
1220  // View of the null space basis columns of X.
1221  // normalizeOutOfPlace() wrote them into X_out.
1222  RCP<MV> X_out_null = MVT::CloneViewNonConst (X_out, nullSpaceIndices);
1223  // Use X_in for scratch space. Copy X_out_null into the
1224  // last few columns of X_in (X_in_null) and do projections
1225  // in there. (This saves us a copy wen we renormalize
1226  // (out of place) back into X_out.)
1227  RCP<MV> X_in_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
1228  MVT::Assign (*X_out_null, *X_in_null);
1229  // Project the new random vectors against the Q blocks, and
1230  // renormalize the result into X_out_null.
1231  rawProject (*X_in_null, Q, C_null);
1232  randomVectorsRank = normalizeOutOfPlace (*X_in_null, *X_out_null, B_null);
1233  } else {
1234  // View of the null space columns of X.
1235  // They live in X_in.
1236  RCP<MV> X_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
1237  // Project the new random vectors against the Q blocks,
1238  // and renormalize the result (in place).
1239  rawProject (*X_null, Q, C_null);
1240  randomVectorsRank = normalize (*X_null, B_null);
1241  }
1242  // While unusual, it is still possible for the random data not
1243  // to be full rank after projection and normalization. In that
1244  // case, we could try another set of random data and recurse as
1245  // necessary, but instead for now we just raise an exception.
1246  TEUCHOS_TEST_FOR_EXCEPTION(randomVectorsRank != numNullSpaceCols,
1247  TsqrOrthoError,
1248  "Belos::TsqrOrthoManagerImpl::projectAndNormalize"
1249  "OutOfPlace(): After projecting and normalizing the "
1250  "random vectors (used to replace the null space "
1251  "basis vectors from normalizing X), they have rank "
1252  << randomVectorsRank << ", but should have full "
1253  "rank " << numNullSpaceCols << ".");
1254  }
1255 
1256  // Whether or not X_in was full rank after projection, we still
1257  // might want to reorthogonalize against Q.
1258  if (reorthogonalizeBlocks_) {
1259  // We are only interested in the column space basis of X
1260  // resp. X_out.
1261  std::vector<magnitude_type>
1262  normsAfterFirstPass (firstPassRank, SCTM::zero());
1263  std::vector<magnitude_type>
1264  normsAfterSecondPass (firstPassRank, SCTM::zero());
1265 
1266  // Compute post-first-pass (pre-normalization) norms. We could
1267  // have done that using MVT::MvNorm() on X_in after projecting,
1268  // but before the first normalization. However, that operation
1269  // may be expensive. It is also unnecessary: after calling
1270  // normalize(), the 2-norm of B(:,j) is the 2-norm of X_in(:,j)
1271  // before normalization, in exact arithmetic.
1272  //
1273  // NOTE (mfh 06 Nov 2010) This is one way that combining
1274  // projection and normalization into a single kernel --
1275  // projectAndNormalize() -- pays off. In project(), we have to
1276  // compute column norms of X before and after projection. Here,
1277  // we get them for free from the normalization coefficients.
1279  for (int j = 0; j < firstPassRank; ++j) {
1280  const Scalar* const B_j = &(*B_out)(0,j);
1281  // Teuchos::BLAS::NRM2 returns a magnitude_type result on
1282  // Scalar inputs.
1283  normsAfterFirstPass[j] = blas.NRM2 (firstPassRank, B_j, 1);
1284  }
1285  // Test whether any of the norms dropped below the
1286  // reorthogonalization threshold.
1287  bool reorthogonalize = false;
1288  for (int j = 0; j < firstPassRank; ++j) {
1289  // If any column's norm decreased too much, mark this block
1290  // for reorthogonalization. Note that this test will _not_
1291  // activate reorthogonalization if a column's norm before the
1292  // first project-and-normalize step was zero. It _will_
1293  // activate reorthogonalization if the column's norm before
1294  // was not zero, but is zero now.
1295  const magnitude_type curThreshold =
1296  blockReorthogThreshold() * normsBeforeFirstPass[j];
1297  if (normsAfterFirstPass[j] < curThreshold) {
1298  reorthogonalize = true;
1299  break;
1300  }
1301  }
1302 
1303  // Notify the caller via callback about the need for
1304  // reorthogonalization.
1305  if (! reorthogCallback_.is_null()) {
1306  using Teuchos::arrayViewFromVector;
1307  (*reorthogCallback_) (arrayViewFromVector (normsBeforeFirstPass),
1308  arrayViewFromVector (normsAfterFirstPass));
1309  }
1310 
1311  // Perform another Block Gram-Schmidt pass if necessary. "Twice
1312  // is enough" (Kahan's theorem) for a Krylov method, unless
1313  // (using Stewart's term) there is an "orthogonalization fault"
1314  // (indicated by reorthogFault).
1315  //
1316  // NOTE (mfh 07 Nov 2010) For now, we include the entire block
1317  // of X, including possible random data (that was already
1318  // projected and normalized above). It might make more sense
1319  // just to process the first firstPassRank columns of X.
1320  // However, the resulting reorthogonalization should still be
1321  // correct regardless.
1322  bool reorthogFault = false;
1323  // Indices of X at which there was an orthogonalization fault.
1324  std::vector<int> faultIndices;
1325  if (reorthogonalize) {
1326  using Teuchos::Copy;
1327  using Teuchos::NO_TRANS;
1328 
1329  // If we're using out-of-place normalization, copy X_out
1330  // (results of first project and normalize pass) back into
1331  // X_in, for the second project and normalize pass.
1332  if (outOfPlace) {
1333  MVT::Assign (X_out, X_in);
1334  }
1335 
1336  // C2 is only used internally, so we know that we are
1337  // allocating fresh and not recycling allocations. Stating
1338  // this lets us save time checking dimensions.
1340  allocateProjectionCoefficients (C2, Q, X_in, false);
1341 
1342  // Block Gram-Schmidt (again). Delay updating the block
1343  // coefficients until we have the new normalization
1344  // coefficients, which we need in order to do the update.
1345  rawProject (X_in, Q, C2);
1346 
1347  // Coefficients for (re)normalization of X_in.
1348  RCP<mat_type> B2 (new mat_type (ncols_X, ncols_X));
1349 
1350  // Normalize X_in (into X_out, if working out of place).
1351  const int secondPassRank = outOfPlace ?
1352  normalizeOutOfPlace (X_in, X_out, B2) :
1353  normalize (X_in, B2);
1354  rank = secondPassRank; // Current rank of X
1355 
1356  // Update normalization coefficients. We begin with copying
1357  // B_out, since the BLAS' _GEMM routine doesn't let us alias
1358  // its input and output arguments.
1359  mat_type B_copy (Copy, *B_out, B_out->numRows(), B_out->numCols());
1360  // B_out := B2 * B_out (where input B_out is in B_copy).
1361  const int err = B_out->multiply (NO_TRANS, NO_TRANS, SCT::one(),
1362  *B2, B_copy, SCT::zero());
1363  TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::logic_error,
1364  "Teuchos::SerialDenseMatrix::multiply "
1365  "returned err = " << err << " != 0");
1366  // Update the block coefficients from the projection step. We
1367  // use B_copy for this (a copy of B_out, the first-pass
1368  // normalization coefficients).
1369  for (int k = 0; k < num_Q_blocks; ++k) {
1370  mat_type C_k_copy (Copy, *C[k], C[k]->numRows(), C[k]->numCols());
1371 
1372  // C[k] := C2[k]*B_copy + C[k].
1373  const int err1 = C[k]->multiply (NO_TRANS, NO_TRANS, SCT::one(),
1374  *C2[k], B_copy, SCT::one());
1375  TEUCHOS_TEST_FOR_EXCEPTION(err1 != 0, std::logic_error,
1376  "Teuchos::SerialDenseMatrix::multiply "
1377  "returned err = " << err1 << " != 0");
1378  }
1379  // Compute post-second-pass (pre-normalization) norms, using
1380  // B2 (the coefficients from the second normalization step) in
1381  // the same way as with B_out before.
1382  for (int j = 0; j < rank; ++j) {
1383  const Scalar* const B2_j = &(*B2)(0,j);
1384  normsAfterSecondPass[j] = blas.NRM2 (rank, B2_j, 1);
1385  }
1386  // Test whether any of the norms dropped below the
1387  // reorthogonalization threshold. If so, it's an
1388  // orthogonalization fault, which requires expensive recovery.
1389  reorthogFault = false;
1390  for (int j = 0; j < rank; ++j) {
1391  const magnitude_type relativeLowerBound =
1392  blockReorthogThreshold() * normsAfterFirstPass[j];
1393  if (normsAfterSecondPass[j] < relativeLowerBound) {
1394  reorthogFault = true;
1395  faultIndices.push_back (j);
1396  }
1397  }
1398  } // if (reorthogonalize) // reorthogonalization pass
1399 
1400  if (reorthogFault) {
1401  if (throwOnReorthogFault_) {
1402  raiseReorthogFault (normsAfterFirstPass,
1403  normsAfterSecondPass,
1404  faultIndices);
1405  } else {
1406  // NOTE (mfh 19 Jan 2011) We could handle the fault here by
1407  // slowly reorthogonalizing, one vector at a time, the
1408  // offending vectors of X. However, we choose not to
1409  // implement this for now. If it becomes a problem, let us
1410  // know and we will prioritize implementing this.
1411  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
1412  "TsqrOrthoManagerImpl has not yet implemented"
1413  " recovery from an orthogonalization fault.");
1414  }
1415  }
1416  } // if (reorthogonalizeBlocks_)
1417  return rank;
1418  }
1419 
1420 
1421  template<class Scalar, class MV>
1422  void
1423  TsqrOrthoManagerImpl<Scalar, MV>::
1424  raiseReorthogFault (const std::vector<magnitude_type>& normsAfterFirstPass,
1425  const std::vector<magnitude_type>& normsAfterSecondPass,
1426  const std::vector<int>& faultIndices)
1427  {
1428  using std::endl;
1429  typedef std::vector<int>::size_type size_type;
1430  std::ostringstream os;
1431 
1432  os << "Orthogonalization fault at the following column(s) of X:" << endl;
1433  os << "Column\tNorm decrease factor" << endl;
1434  for (size_type k = 0; k < faultIndices.size(); ++k) {
1435  const int index = faultIndices[k];
1436  const magnitude_type decreaseFactor =
1437  normsAfterSecondPass[index] / normsAfterFirstPass[index];
1438  os << index << "\t" << decreaseFactor << endl;
1439  }
1440  throw TsqrOrthoFault (os.str());
1441  }
1442 
1443  template<class Scalar, class MV>
1446  {
1447  using Teuchos::ParameterList;
1448  using Teuchos::parameterList;
1449  using Teuchos::RCP;
1450 
1451  if (defaultParams_.is_null()) {
1452  RCP<ParameterList> params = parameterList ("TsqrOrthoManagerImpl");
1453  //
1454  // TSQR parameters (set as a sublist).
1455  //
1456  params->set ("TSQR implementation", *(tsqrAdaptor_.getValidParameters()),
1457  "TSQR implementation parameters.");
1458  //
1459  // Orthogonalization parameters
1460  //
1461  const bool defaultRandomizeNullSpace = true;
1462  params->set ("randomizeNullSpace", defaultRandomizeNullSpace,
1463  "Whether to fill in null space vectors with random data.");
1464 
1465  const bool defaultReorthogonalizeBlocks = true;
1466  params->set ("reorthogonalizeBlocks", defaultReorthogonalizeBlocks,
1467  "Whether to do block reorthogonalization as necessary.");
1468 
1469  // This parameter corresponds to the "blk_tol_" parameter in
1470  // Belos' DGKSOrthoManager. We choose the same default value.
1471  const magnitude_type defaultBlockReorthogThreshold =
1472  magnitude_type(10) * SCTM::squareroot (SCTM::eps());
1473  params->set ("blockReorthogThreshold", defaultBlockReorthogThreshold,
1474  "If reorthogonalizeBlocks==true, and if the norm of "
1475  "any column within a block decreases by this much or "
1476  "more after orthogonalization, we reorthogonalize.");
1477 
1478  // This parameter corresponds to the "sing_tol_" parameter in
1479  // Belos' DGKSOrthoManager. We choose the same default value.
1480  const magnitude_type defaultRelativeRankTolerance =
1481  Teuchos::as<magnitude_type>(10) * SCTM::eps();
1482 
1483  // If the relative rank tolerance is zero, then we will always
1484  // declare blocks to be numerically full rank, as long as no
1485  // singular values are zero.
1486  params->set ("relativeRankTolerance", defaultRelativeRankTolerance,
1487  "Relative tolerance to determine the numerical rank of a "
1488  "block when normalizing.");
1489 
1490  // See Stewart's 2008 paper on block Gram-Schmidt for a definition
1491  // of "orthogonalization fault."
1492  const bool defaultThrowOnReorthogFault = true;
1493  params->set ("throwOnReorthogFault", defaultThrowOnReorthogFault,
1494  "Whether to throw an exception if an orthogonalization "
1495  "fault occurs. This only matters if reorthogonalization "
1496  "is enabled (reorthogonalizeBlocks==true).");
1497 
1498  const bool defaultForceNonnegativeDiagonal = false;
1499  params->set ("forceNonnegativeDiagonal", defaultForceNonnegativeDiagonal,
1500  "Whether to force the R factor produced by the normalization "
1501  "step to have a nonnegative diagonal.");
1502 
1503  defaultParams_ = params;
1504  }
1505  return defaultParams_;
1506  }
1507 
1508  template<class Scalar, class MV>
1511  {
1512  using Teuchos::ParameterList;
1513  using Teuchos::RCP;
1514  using Teuchos::rcp;
1515 
1516  RCP<const ParameterList> defaultParams = getValidParameters();
1517  // Start with a clone of the default parameters.
1518  RCP<ParameterList> params = rcp (new ParameterList (*defaultParams));
1519 
1520  // Disable reorthogonalization and randomization of the null
1521  // space basis. Reorthogonalization tolerances don't matter,
1522  // since we aren't reorthogonalizing blocks in the fast
1523  // settings. We can leave the default values. Also,
1524  // (re)orthogonalization faults may only occur with
1525  // reorthogonalization, so we don't have to worry about the
1526  // "throwOnReorthogFault" setting.
1527  const bool randomizeNullSpace = false;
1528  params->set ("randomizeNullSpace", randomizeNullSpace);
1529  const bool reorthogonalizeBlocks = false;
1530  params->set ("reorthogonalizeBlocks", reorthogonalizeBlocks);
1531 
1532  return params;
1533  }
1534 
1535  template<class Scalar, class MV>
1536  int
1538  rawNormalize (MV& X,
1539  MV& Q,
1541  {
1542  int rank;
1543  try {
1544  // This call only computes the QR factorization X = Q B.
1545  // It doesn't compute the rank of X. That comes from
1546  // revealRank() below.
1547  tsqrAdaptor_.factorExplicit (X, Q, B, forceNonnegativeDiagonal_);
1548  // This call will only modify *B if *B on input is not of full
1549  // numerical rank.
1550  rank = tsqrAdaptor_.revealRank (Q, B, relativeRankTolerance_);
1551  } catch (std::exception& e) {
1552  throw TsqrOrthoError (e.what()); // Toss the exception up the chain.
1553  }
1554  return rank;
1555  }
1556 
1557  template<class Scalar, class MV>
1558  int
1559  TsqrOrthoManagerImpl<Scalar, MV>::
1560  normalizeOne (MV& X,
1562  {
1563  // Make space for the normalization coefficient. This will either
1564  // be a freshly allocated matrix (if B is null), or a view of the
1565  // 1x1 upper left submatrix of *B (if B is not null).
1566  mat_ptr B_out;
1567  if (B.is_null()) {
1568  B_out = Teuchos::rcp (new mat_type (1, 1));
1569  } else {
1570  const int theNumRows = B->numRows ();
1571  const int theNumCols = B->numCols ();
1573  theNumRows < 1 || theNumCols < 1, std::invalid_argument,
1574  "normalizeOne: Input matrix B must be at least 1 x 1, but "
1575  "is instead " << theNumRows << " x " << theNumCols << ".");
1576  // Create a view of the 1x1 upper left submatrix of *B.
1577  B_out = Teuchos::rcp (new mat_type (Teuchos::View, *B, 1, 1));
1578  }
1579 
1580  // Compute the norm of X, and write the result to B_out.
1581  std::vector<magnitude_type> theNorm (1, SCTM::zero());
1582  MVT::MvNorm (X, theNorm);
1583  (*B_out)(0,0) = theNorm[0];
1584 
1585  if (B.is_null()) {
1586  // The input matrix B is null, so assign B_out to it. If B was
1587  // not null on input, then B_out is a view of *B, so we don't
1588  // have to do anything here. Note that SerialDenseMatrix uses
1589  // raw pointers to store data and represent views, so we have to
1590  // be careful about scope.
1591  B = B_out;
1592  }
1593 
1594  // Scale X by its norm, if its norm is zero. Otherwise, do the
1595  // right thing based on whether the user wants us to fill the null
1596  // space with random vectors.
1597  if (theNorm[0] == SCTM::zero()) {
1598  // Make a view of the first column of Q, fill it with random
1599  // data, and normalize it. Throw away the resulting norm.
1600  if (randomizeNullSpace_) {
1601  MVT::MvRandom(X);
1602  MVT::MvNorm (X, theNorm);
1603  if (theNorm[0] == SCTM::zero()) {
1604  // It is possible that a random vector could have all zero
1605  // entries, but unlikely. We could try again, but it's also
1606  // possible that multiple tries could result in zero
1607  // vectors. We choose instead to give up.
1608  throw TsqrOrthoError("normalizeOne: a supposedly random "
1609  "vector has norm zero!");
1610  } else {
1611  // NOTE (mfh 09 Nov 2010) I'm assuming that dividing a
1612  // Scalar by a magnitude_type is defined and that it results
1613  // in a Scalar.
1614  const Scalar alpha = SCT::one() / theNorm[0];
1615  MVT::MvScale (X, alpha);
1616  }
1617  }
1618  return 0; // The rank of the matrix (actually just one vector) X.
1619  } else {
1620  // NOTE (mfh 09 Nov 2010) I'm assuming that dividing a Scalar by
1621  // a magnitude_type is defined and that it results in a Scalar.
1622  const Scalar alpha = SCT::one() / theNorm[0];
1623  MVT::MvScale (X, alpha);
1624  return 1; // The rank of the matrix (actually just one vector) X.
1625  }
1626  }
1627 
1628 
1629  template<class Scalar, class MV>
1630  void
1631  TsqrOrthoManagerImpl<Scalar, MV>::
1632  rawProject (MV& X,
1635  {
1636 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1637  Teuchos::TimeMonitor timerMonitorNormalize(*timerProject_);
1638 #endif // BELOS_TEUCHOS_TIME_MONITOR
1639 
1640  // "Modified Gram-Schmidt" version of Block Gram-Schmidt.
1641  const int num_Q_blocks = Q.size();
1642  for (int i = 0; i < num_Q_blocks; ++i)
1643  {
1644  // TEUCHOS_TEST_FOR_EXCEPTION(C[i].is_null(), std::logic_error,
1645  // "TsqrOrthoManagerImpl::rawProject(): C["
1646  // << i << "] is null");
1647  // TEUCHOS_TEST_FOR_EXCEPTION(Q[i].is_null(), std::logic_error,
1648  // "TsqrOrthoManagerImpl::rawProject(): Q["
1649  // << i << "] is null");
1650  mat_type& Ci = *C[i];
1651  const MV& Qi = *Q[i];
1652  innerProd (Qi, X, Ci);
1653  MVT::MvTimesMatAddMv (-SCT::one(), Qi, Ci, SCT::one(), X);
1654  }
1655  }
1656 
1657 
1658  template<class Scalar, class MV>
1659  void
1660  TsqrOrthoManagerImpl<Scalar, MV>::
1661  rawProject (MV& X,
1662  const Teuchos::RCP<const MV>& Q,
1664  {
1665 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1666  Teuchos::TimeMonitor timerMonitorNormalize(*timerProject_);
1667 #endif // BELOS_TEUCHOS_TIME_MONITOR
1668 
1669  // Block Gram-Schmidt
1670  innerProd (*Q, X, *C);
1671  MVT::MvTimesMatAddMv (-SCT::one(), *Q, *C, SCT::one(), X);
1672  }
1673 
1674  template<class Scalar, class MV>
1675  int
1676  TsqrOrthoManagerImpl<Scalar, MV>::
1677  normalizeImpl (MV& X,
1678  MV& Q,
1680  const bool outOfPlace)
1681  {
1682  using Teuchos::Range1D;
1683  using Teuchos::RCP;
1684  using Teuchos::rcp;
1685  using Teuchos::ScalarTraits;
1686  using Teuchos::tuple;
1687 
1688  const int numCols = MVT::GetNumberVecs (X);
1689  if (numCols == 0) {
1690  return 0; // Fast exit for an empty input matrix.
1691  }
1692 
1693  // We allow Q to have more columns than X. In that case, we only
1694  // touch the first numCols columns of Q.
1696  MVT::GetNumberVecs (Q) < numCols, std::invalid_argument,
1697  "TsqrOrthoManagerImpl::normalizeImpl: Q has "
1698  << MVT::GetNumberVecs(Q) << " columns. This is too "
1699  "few, since X has " << numCols << " columns.");
1700  // TSQR wants a Q with the same number of columns as X, so have it
1701  // work on a nonconstant view of Q with the same number of columns
1702  // as X.
1703  RCP<MV> Q_view = MVT::CloneViewNonConst (Q, Range1D (0, numCols-1));
1704 
1705  // Make space for the normalization coefficients. This will
1706  // either be a freshly allocated matrix (if B is null), or a view
1707  // of the appropriately sized upper left submatrix of *B (if B is
1708  // not null).
1709  mat_ptr B_out;
1710  if (B.is_null ()) {
1711  B_out = rcp (new mat_type (numCols, numCols));
1712  } else {
1713  // Make sure that B is no smaller than numCols x numCols.
1715  B->numRows() < numCols || B->numCols() < numCols, std::invalid_argument,
1716  "TsqrOrthoManagerImpl::normalizeImpl: Input matrix B must be at least "
1717  << numCols << " x " << numCols << ", but is instead " << B->numRows ()
1718  << " x " << B->numCols() << ".");
1719  // Create a view of the numCols x numCols upper left submatrix
1720  // of *B. TSQR will write the normalization coefficients there.
1721  B_out = rcp (new mat_type (Teuchos::View, *B, numCols, numCols));
1722  }
1723 
1724  // Compute rank-revealing decomposition (in this case, TSQR of X
1725  // followed by SVD of the R factor and appropriate updating of the
1726  // resulting Q factor) of X. X is modified in place and filled
1727  // with garbage, and Q_view contains the resulting explicit Q
1728  // factor. Later, we will copy this back into X.
1729  //
1730  // The matrix *B_out will only be upper triangular if X is of full
1731  // numerical rank. Otherwise, the entries below the diagonal may
1732  // be filled in as well.
1733  const int rank = rawNormalize (X, *Q_view, *B_out);
1734  if (B.is_null ()) {
1735  // The input matrix B is null, so assign B_out to it. If B was
1736  // not null on input, then B_out is a view of *B, so we don't
1737  // have to do anything here. Note that SerialDenseMatrix uses
1738  // raw pointers to store data and represent views, so we have to
1739  // be careful about scope.
1740  B = B_out;
1741  }
1742 
1744  rank < 0 || rank > numCols, std::logic_error,
1745  "Belos::TsqrOrthoManagerImpl::normalizeImpl: rawNormalize returned rank "
1746  " = " << rank << " for a matrix X with " << numCols << " columns. "
1747  "Please report this bug to the Belos developers.");
1748 
1749  // If X is full rank or we don't want to replace its null space
1750  // basis with random vectors, then we're done.
1751  if (rank == numCols || ! randomizeNullSpace_) {
1752  // If we're supposed to be working in place in X, copy the
1753  // results back from Q_view into X.
1754  if (! outOfPlace) {
1755  MVT::Assign (*Q_view, X);
1756  }
1757  return rank;
1758  }
1759 
1760  if (randomizeNullSpace_ && rank < numCols) {
1761  // X wasn't full rank. Replace the null space basis of X (in
1762  // the last numCols-rank columns of Q_view) with random data,
1763  // project it against the first rank columns of Q_view, and
1764  // normalize.
1765  //
1766  // Number of columns to fill with random data.
1767  const int nullSpaceNumCols = numCols - rank;
1768  // Inclusive range of indices of columns of X to fill with
1769  // random data.
1770  Range1D nullSpaceIndices (rank, numCols-1);
1771 
1772  // rawNormalize wrote the null space basis vectors into Q_view.
1773  // We continue to work in place in Q_view by writing the random
1774  // data there and (if there is a nontrival column space)
1775  // projecting in place against the column space basis vectors
1776  // (also in Q_view).
1777  RCP<MV> Q_null = MVT::CloneViewNonConst (*Q_view, nullSpaceIndices);
1778  // Replace the null space basis with random data.
1779  MVT::MvRandom (*Q_null);
1780 
1781  // Make sure that the "random" data isn't all zeros. This is
1782  // statistically nearly impossible, but we test for debugging
1783  // purposes.
1784  {
1785  std::vector<magnitude_type> norms (MVT::GetNumberVecs (*Q_null));
1786  MVT::MvNorm (*Q_null, norms);
1787 
1788  bool anyZero = false;
1789  typedef typename std::vector<magnitude_type>::const_iterator iter_type;
1790  for (iter_type it = norms.begin(); it != norms.end(); ++it) {
1791  if (*it == SCTM::zero()) {
1792  anyZero = true;
1793  }
1794  }
1795  if (anyZero) {
1796  std::ostringstream os;
1797  os << "TsqrOrthoManagerImpl::normalizeImpl: "
1798  "We are being asked to randomize the null space, for a matrix "
1799  "with " << numCols << " columns and reported column rank "
1800  << rank << ". The inclusive range of columns to fill with "
1801  "random data is [" << nullSpaceIndices.lbound() << ","
1802  << nullSpaceIndices.ubound() << "]. After filling the null "
1803  "space vectors with random numbers, at least one of the vectors"
1804  " has norm zero. Here are the norms of all the null space "
1805  "vectors: [";
1806  for (iter_type it = norms.begin(); it != norms.end(); ++it) {
1807  os << *it;
1808  if (it+1 != norms.end())
1809  os << ", ";
1810  }
1811  os << "].) There is a tiny probability that this could happen "
1812  "randomly, but it is likely a bug. Please report it to the "
1813  "Belos developers, especially if you are able to reproduce the "
1814  "behavior.";
1815  TEUCHOS_TEST_FOR_EXCEPTION(anyZero, TsqrOrthoError, os.str());
1816  }
1817  }
1818 
1819  if (rank > 0) {
1820  // Project the random data against the column space basis of
1821  // X, using a simple block projection ("Block Classical
1822  // Gram-Schmidt"). This is accurate because we've already
1823  // orthogonalized the column space basis of X nearly to
1824  // machine precision via a QR factorization (TSQR) with
1825  // accuracy comparable to Householder QR.
1826  RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D (0, rank-1));
1827 
1828  // Temporary storage for projection coefficients. We don't
1829  // need to keep them, since they represent the null space
1830  // basis (for which the coefficients are logically zero).
1831  mat_ptr C_null (new mat_type (rank, nullSpaceNumCols));
1832  rawProject (*Q_null, Q_col, C_null);
1833  }
1834  // Normalize the projected random vectors, so that they are
1835  // mutually orthogonal (as well as orthogonal to the column
1836  // space basis of X). We use X for the output of the
1837  // normalization: for out-of-place normalization (outOfPlace ==
1838  // true), X is overwritten with "invalid values" anyway, and for
1839  // in-place normalization (outOfPlace == false), we want the
1840  // result to be in X anyway.
1841  RCP<MV> X_null = MVT::CloneViewNonConst (X, nullSpaceIndices);
1842  // Normalization coefficients for projected random vectors.
1843  // Will be thrown away.
1844  mat_type B_null (nullSpaceNumCols, nullSpaceNumCols);
1845  // Write the normalized vectors to X_null (in X).
1846  const int nullSpaceBasisRank = rawNormalize (*Q_null, *X_null, B_null);
1847 
1848  // It's possible, but unlikely, that X_null doesn't have full
1849  // rank (after the projection step). We could recursively fill
1850  // in more random vectors until we finally get a full rank
1851  // matrix, but instead we just throw an exception.
1852  //
1853  // NOTE (mfh 08 Nov 2010) Perhaps we should deal with this case
1854  // more elegantly. Recursion might be one way to solve it, but
1855  // be sure to check that the recursion will terminate. We could
1856  // do this by supplying an additional argument to rawNormalize,
1857  // which is the null space basis rank from the previous
1858  // iteration. The rank has to decrease each time, or the
1859  // recursion may go on forever.
1860  if (nullSpaceBasisRank < nullSpaceNumCols) {
1861  std::vector<magnitude_type> norms (MVT::GetNumberVecs(*X_null));
1862  MVT::MvNorm (*X_null, norms);
1863  std::ostringstream os;
1864  os << "TsqrOrthoManagerImpl::normalizeImpl: "
1865  << "We are being asked to randomize the null space, "
1866  << "for a matrix with " << numCols << " columns and "
1867  << "column rank " << rank << ". After projecting and "
1868  << "normalizing the generated random vectors, they "
1869  << "only have rank " << nullSpaceBasisRank << ". They"
1870  << " should have full rank " << nullSpaceNumCols
1871  << ". (The inclusive range of columns to fill with "
1872  << "random data is [" << nullSpaceIndices.lbound()
1873  << "," << nullSpaceIndices.ubound() << "]. The "
1874  << "column norms of the resulting Q factor are: [";
1875  for (typename std::vector<magnitude_type>::size_type k = 0;
1876  k < norms.size(); ++k) {
1877  os << norms[k];
1878  if (k != norms.size()-1) {
1879  os << ", ";
1880  }
1881  }
1882  os << "].) There is a tiny probability that this could "
1883  << "happen randomly, but it is likely a bug. Please "
1884  << "report it to the Belos developers, especially if "
1885  << "you are able to reproduce the behavior.";
1886 
1887  TEUCHOS_TEST_FOR_EXCEPTION(nullSpaceBasisRank < nullSpaceNumCols,
1888  TsqrOrthoError, os.str ());
1889  }
1890  // If we're normalizing out of place, copy the X_null
1891  // vectors back into Q_null; the Q_col vectors are already
1892  // where they are supposed to be in that case.
1893  //
1894  // If we're normalizing in place, leave X_null alone (it's
1895  // already where it needs to be, in X), but copy Q_col back
1896  // into the first rank columns of X.
1897  if (outOfPlace) {
1898  MVT::Assign (*X_null, *Q_null);
1899  } else if (rank > 0) {
1900  // MVT::Assign() doesn't accept empty ranges of columns.
1901  RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D (0, rank-1));
1902  RCP<MV> X_col = MVT::CloneViewNonConst (X, Range1D (0, rank-1));
1903  MVT::Assign (*Q_col, *X_col);
1904  }
1905  }
1906  return rank;
1907  }
1908 
1909 
1910  template<class Scalar, class MV>
1911  void
1912  TsqrOrthoManagerImpl<Scalar, MV>::
1913  checkProjectionDims (int& ncols_X,
1914  int& num_Q_blocks,
1915  int& ncols_Q_total,
1916  const MV& X,
1918  {
1919  // First assign to temporary values, so the function won't
1920  // commit any externally visible side effects unless it will
1921  // return normally (without throwing an exception). (I'm being
1922  // cautious; MVT::GetNumberVecs() probably shouldn't have any
1923  // externally visible side effects, unless it is logging to a
1924  // file or something.)
1925  int the_ncols_X, the_num_Q_blocks, the_ncols_Q_total;
1926  the_num_Q_blocks = Q.size();
1927  the_ncols_X = MVT::GetNumberVecs (X);
1928 
1929  // Compute the total number of columns of all the Q[i] blocks.
1930  the_ncols_Q_total = 0;
1931  // You should be angry if your compiler doesn't support type
1932  // inference ("auto"). That's why I need this awful typedef.
1933  using Teuchos::ArrayView;
1934  using Teuchos::RCP;
1935  typedef typename ArrayView<RCP<const MV> >::const_iterator iter_type;
1936  for (iter_type it = Q.begin(); it != Q.end(); ++it) {
1937  const MV& Qi = **it;
1938  the_ncols_Q_total += MVT::GetNumberVecs (Qi);
1939  }
1940 
1941  // Commit temporary values to the output arguments.
1942  ncols_X = the_ncols_X;
1943  num_Q_blocks = the_num_Q_blocks;
1944  ncols_Q_total = the_ncols_Q_total;
1945  }
1946 
1947 } // namespace Belos
1948 
1949 #endif // __BelosTsqrOrthoManagerImpl_hpp
1950 
int projectAndNormalize(MV &X, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Project X against Q and normalize X.
Interface of callback invoked by TsqrOrthoManager on reorthogonalization.
bool is_null(const boost::shared_ptr< T > &p)
void project(MV &X, Teuchos::Array< mat_ptr > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Compute and .
void norm(const MV &X, std::vector< magnitude_type > &normVec) const
Compute the 2-norm of each column j of X.
TsqrOrthoManagerImpl(const Teuchos::RCP< Teuchos::ParameterList > &params, const std::string &label)
Constructor (that sets user-specified parameters).
static void clearCounter(const std::string &name)
magnitude_type orthogError(const MV &X1, const MV &X2) const
Return the Frobenius norm of the inner product of X1 with itself.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
Type of the projection and normalization coefficients.
T & get(const std::string &name, T def_value)
static RCP< Time > getNewCounter(const std::string &name)
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &params)
Set parameters from the given parameter list.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters()
Get &quot;fast&quot; parameters for TsqrOrthoManagerImpl.
TSQR-based OrthoManager subclass implementation.
Declaration of basic traits for the multivector type.
void setReorthogonalizationCallback(const Teuchos::RCP< ReorthogonalizationCallback< Scalar > > &callback)
Set callback to be invoked on reorthogonalization.
Error in TsqrOrthoManager or TsqrOrthoManagerImpl.
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
int normalize(MV &X, mat_ptr B)
Orthogonalize the columns of X in place.
void setLabel(const std::string &label)
Set the label for timers.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
The type of a norm result.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Default valid parameter list.
Traits class which defines basic operations on multivectors.
Scalar scalar_type
The template parameter of this class; the type of an inner product result.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual void operator()(Teuchos::ArrayView< magnitude_type > normsBeforeFirstPass, Teuchos::ArrayView< magnitude_type > normsAfterFirstPass)=0
Callback invoked by TsqrOrthoManager on reorthogonalization.
const std::string & getLabel() const
Get the label for timers (if timers are enabled).
magnitude_type orthonormError(const MV &X) const
Return .
Exception thrown to signal error in an orthogonalization manager method.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
magnitude_type relativeRankTolerance() const
Relative tolerance for determining (via the SVD) whether a block is of full numerical rank...
void resize(size_type new_size, const value_type &x=value_type())
TsqrOrthoError(const std::string &what_arg)
OrdinalType numCols() const
virtual ~ReorthogonalizationCallback()
Destructor (virtual for memory safety of derived classes)
Templated virtual class for providing orthogonalization/orthonormalization methods.
TsqrOrthoFault(const std::string &what_arg)
magnitude_type blockReorthogThreshold() const
Relative tolerance for triggering a block reorthogonalization.
int projectAndNormalizeOutOfPlace(MV &X_in, MV &X_out, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Project and normalize X_in into X_out; overwrite X_in.
int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B)
Normalize X into Q*B, overwriting X.
void innerProd(const MV &X, const MV &Y, mat_type &Z) const
Euclidean inner product.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Orthogonalization fault.
OrdinalType numRows() const
bool is_null() const

Generated on Fri Aug 14 2020 10:48:34 for Belos by doxygen 1.8.5