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