Anasazi  Version of the Day
AnasaziTsqrOrthoManagerImpl.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2010) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
32 #ifndef __AnasaziTsqrOrthoManagerImpl_hpp
33 #define __AnasaziTsqrOrthoManagerImpl_hpp
34 
35 #include "AnasaziConfigDefs.hpp"
37 #include "AnasaziOrthoManager.hpp" // OrthoError, etc.
38 
39 #include "Teuchos_as.hpp"
40 #include "Teuchos_LAPACK.hpp"
41 #include "Teuchos_ParameterList.hpp"
42 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
43 #include <algorithm>
44 
45 
46 namespace Anasazi {
47 
51  class TsqrOrthoError : public OrthoError {
52  public:
53  TsqrOrthoError (const std::string& what_arg) :
54  OrthoError (what_arg) {}
55  };
56 
76  class TsqrOrthoFault : public OrthoError {
77  public:
78  TsqrOrthoFault (const std::string& what_arg) :
79  OrthoError (what_arg) {}
80  };
81 
114  template<class Scalar, class MV>
116  public Teuchos::ParameterListAcceptorDefaultBase {
117  public:
118  typedef Scalar scalar_type;
119  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
120  typedef MV multivector_type;
123  typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
124  typedef Teuchos::RCP<mat_type> mat_ptr;
125 
126  private:
127  typedef Teuchos::ScalarTraits<Scalar> SCT;
128  typedef Teuchos::ScalarTraits<magnitude_type> SCTM;
130  typedef typename MVT::tsqr_adaptor_type tsqr_adaptor_type;
131 
132  public:
140  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters () const;
141 
143  void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params);
144 
155  Teuchos::RCP<const Teuchos::ParameterList> getFastParameters ();
156 
173  TsqrOrthoManagerImpl (const Teuchos::RCP<Teuchos::ParameterList>& params,
174  const std::string& label);
175 
180  TsqrOrthoManagerImpl (const std::string& label);
181 
189  void setLabel (const std::string& label) {
190  if (label != label_) {
191  label_ = label;
192  }
193  }
194 
196  const std::string& getLabel () const { return label_; }
197 
206  void
207  innerProd (const MV& X, const MV& Y, mat_type& Z) const
208  {
209  MVT::MvTransMv (SCT::one(), X, Y, Z);
210  }
211 
229  void
230  norm (const MV& X, std::vector<magnitude_type>& normvec) const;
231 
241  void
242  project (MV& X,
243  Teuchos::Array<mat_ptr> C,
244  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q);
245 
259  int normalize (MV& X, mat_ptr B);
260 
279  int
280  normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B);
281 
295  int
297  Teuchos::Array<mat_ptr> C,
298  mat_ptr B,
299  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
300  {
301  // "false" means we work on X in place. The second argument is
302  // not read or written in that case.
303  return projectAndNormalizeImpl (X, X, false, C, B, Q);
304  }
305 
325  int
327  MV& X_out,
328  Teuchos::Array<mat_ptr> C,
329  mat_ptr B,
330  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
331  {
332  // "true" means we work on X_in out of place, writing the
333  // results into X_out.
334  return projectAndNormalizeImpl (X_in, X_out, true, C, B, Q);
335  }
336 
341  magnitude_type
342  orthonormError (const MV &X) const
343  {
344  const Scalar ONE = SCT::one();
345  const int ncols = MVT::GetNumberVecs(X);
346  mat_type XTX (ncols, ncols);
347  innerProd (X, X, XTX);
348  for (int k = 0; k < ncols; ++k) {
349  XTX(k,k) -= ONE;
350  }
351  return XTX.normFrobenius();
352  }
353 
355  magnitude_type
356  orthogError (const MV &X1,
357  const MV &X2) const
358  {
359  const int ncols_X1 = MVT::GetNumberVecs (X1);
360  const int ncols_X2 = MVT::GetNumberVecs (X2);
361  mat_type X1_T_X2 (ncols_X1, ncols_X2);
362  innerProd (X1, X2, X1_T_X2);
363  return X1_T_X2.normFrobenius();
364  }
365 
369  magnitude_type blockReorthogThreshold() const { return blockReorthogThreshold_; }
370 
373  magnitude_type relativeRankTolerance() const { return relativeRankTolerance_; }
374 
375  private:
377  Teuchos::RCP<Teuchos::ParameterList> params_;
378 
380  mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
381 
383  std::string label_;
384 
386  tsqr_adaptor_type tsqrAdaptor_;
387 
397  Teuchos::RCP<MV> Q_;
398 
400  magnitude_type eps_;
401 
402 
406  bool randomizeNullSpace_;
407 
413  bool reorthogonalizeBlocks_;
414 
418  bool throwOnReorthogFault_;
419 
421  magnitude_type blockReorthogThreshold_;
422 
424  magnitude_type relativeRankTolerance_;
425 
432  bool forceNonnegativeDiagonal_;
433 
435  void
436  raiseReorthogFault (const std::vector<magnitude_type>& normsAfterFirstPass,
437  const std::vector<magnitude_type>& normsAfterSecondPass,
438  const std::vector<int>& faultIndices);
439 
449  void
450  checkProjectionDims (int& ncols_X,
451  int& num_Q_blocks,
452  int& ncols_Q_total,
453  const MV& X,
454  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
455 
466  void
467  allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C,
468  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
469  const MV& X,
470  const bool attemptToRecycle = true) const;
471 
480  int
481  projectAndNormalizeImpl (MV& X_in,
482  MV& X_out,
483  const bool outOfPlace,
484  Teuchos::Array<mat_ptr> C,
485  mat_ptr B,
486  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q);
487 
493  void
494  rawProject (MV& X,
495  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
496  Teuchos::ArrayView<mat_ptr> C) const;
497 
499  void
500  rawProject (MV& X,
501  const Teuchos::RCP<const MV>& Q,
502  const mat_ptr& C) const;
503 
530  int rawNormalize (MV& X, MV& Q, mat_type& B);
531 
549  int normalizeOne (MV& X, mat_ptr B) const;
550 
578  int normalizeImpl (MV& X, MV& Q, mat_ptr B, const bool outOfPlace);
579  };
580 
581  template<class Scalar, class MV>
582  void
584  setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params)
585  {
586  using Teuchos::ParameterList;
587  using Teuchos::parameterList;
588  using Teuchos::RCP;
589  using Teuchos::sublist;
590  typedef magnitude_type M; // abbreviation.
591 
592  RCP<const ParameterList> defaultParams = getValidParameters ();
593  // Sublist of TSQR implementation parameters; to get below.
594  RCP<ParameterList> tsqrParams;
595 
596  RCP<ParameterList> theParams;
597  if (params.is_null()) {
598  theParams = parameterList (*defaultParams);
599  } else {
600  theParams = params;
601 
602  // Don't call validateParametersAndSetDefaults(); we prefer to
603  // ignore parameters that we don't recognize, at least for now.
604  // However, we do fill in missing parameters with defaults.
605 
606  randomizeNullSpace_ =
607  theParams->get<bool> ("randomizeNullSpace",
608  defaultParams->get<bool> ("randomizeNullSpace"));
609  reorthogonalizeBlocks_ =
610  theParams->get<bool> ("reorthogonalizeBlocks",
611  defaultParams->get<bool> ("reorthogonalizeBlocks"));
612  throwOnReorthogFault_ =
613  theParams->get<bool> ("throwOnReorthogFault",
614  defaultParams->get<bool> ("throwOnReorthogFault"));
615  blockReorthogThreshold_ =
616  theParams->get<M> ("blockReorthogThreshold",
617  defaultParams->get<M> ("blockReorthogThreshold"));
618  relativeRankTolerance_ =
619  theParams->get<M> ("relativeRankTolerance",
620  defaultParams->get<M> ("relativeRankTolerance"));
621  forceNonnegativeDiagonal_ =
622  theParams->get<bool> ("forceNonnegativeDiagonal",
623  defaultParams->get<bool> ("forceNonnegativeDiagonal"));
624 
625  // Get the sublist of TSQR implementation parameters. Use the
626  // default sublist if one isn't provided.
627  if (! theParams->isSublist ("TSQR implementation")) {
628  theParams->set ("TSQR implementation",
629  defaultParams->sublist ("TSQR implementation"));
630  }
631  tsqrParams = sublist (theParams, "TSQR implementation", true);
632  }
633 
634  // Send the TSQR implementation parameters to the TSQR adaptor.
635  tsqrAdaptor_.setParameterList (tsqrParams);
636 
637  // Save the input parameter list.
638  setMyParamList (theParams);
639  }
640 
641  template<class Scalar, class MV>
643  TsqrOrthoManagerImpl (const Teuchos::RCP<Teuchos::ParameterList>& params,
644  const std::string& label) :
645  label_ (label),
646  Q_ (Teuchos::null), // Initialized on demand
647  eps_ (SCTM::eps()), // Machine precision
648  randomizeNullSpace_ (true),
649  reorthogonalizeBlocks_ (true),
650  throwOnReorthogFault_ (false),
651  blockReorthogThreshold_ (0),
652  relativeRankTolerance_ (0),
653  forceNonnegativeDiagonal_ (false)
654  {
655  setParameterList (params); // This also sets tsqrAdaptor_'s parameters.
656  }
657 
658  template<class Scalar, class MV>
660  TsqrOrthoManagerImpl (const std::string& label) :
661  label_ (label),
662  Q_ (Teuchos::null), // Initialized on demand
663  eps_ (SCTM::eps()), // Machine precision
664  randomizeNullSpace_ (true),
665  reorthogonalizeBlocks_ (true),
666  throwOnReorthogFault_ (false),
667  blockReorthogThreshold_ (0),
668  relativeRankTolerance_ (0),
669  forceNonnegativeDiagonal_ (false)
670  {
671  setParameterList (Teuchos::null); // Set default parameters.
672  }
673 
674  template<class Scalar, class MV>
675  void
677  norm (const MV& X, std::vector<magnitude_type>& normVec) const
678  {
679  const int numCols = MVT::GetNumberVecs (X);
680  // std::vector<T>::size_type is unsigned; int is signed. Mixed
681  // unsigned/signed comparisons trigger compiler warnings.
682  if (normVec.size() < static_cast<size_t>(numCols))
683  normVec.resize (numCols); // Resize normvec if necessary.
684  MVT::MvNorm (X, normVec);
685  }
686 
687  template<class Scalar, class MV>
688  void
690  Teuchos::Array<mat_ptr> C,
691  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
692  {
693  int ncols_X, num_Q_blocks, ncols_Q_total;
694  checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X, Q);
695  // Test for quick exit: any dimension of X is zero, or there are
696  // zero Q blocks, or the total number of columns of the Q blocks
697  // is zero.
698  if (ncols_X == 0 || num_Q_blocks == 0 || ncols_Q_total == 0)
699  return;
700 
701  // Make space for first-pass projection coefficients
702  allocateProjectionCoefficients (C, Q, X, true);
703 
704  // We only use columnNormsBefore and compute pre-projection column
705  // norms if doing block reorthogonalization.
706  std::vector<magnitude_type> columnNormsBefore (ncols_X, magnitude_type(0));
707  if (reorthogonalizeBlocks_)
708  MVT::MvNorm (X, columnNormsBefore);
709 
710  // Project (first block orthogonalization step):
711  // C := Q^* X, X := X - Q C.
712  rawProject (X, Q, C);
713 
714  // If we are doing block reorthogonalization, reorthogonalize X if
715  // necessary.
716  if (reorthogonalizeBlocks_)
717  {
718  std::vector<magnitude_type> columnNormsAfter (ncols_X, magnitude_type(0));
719  MVT::MvNorm (X, columnNormsAfter);
720 
721  // Relative block reorthogonalization threshold
722  const magnitude_type relThres = blockReorthogThreshold();
723  // Reorthogonalize X if any of its column norms decreased by a
724  // factor more than the block reorthogonalization threshold.
725  // Don't bother trying to subset the columns; that will make
726  // the columns noncontiguous and thus hinder BLAS 3
727  // optimizations.
728  bool reorthogonalize = false;
729  for (int j = 0; j < ncols_X; ++j)
730  if (columnNormsAfter[j] < relThres * columnNormsBefore[j])
731  {
732  reorthogonalize = true;
733  break;
734  }
735  if (reorthogonalize)
736  {
737  // Second-pass projection coefficients
738  Teuchos::Array<mat_ptr> C2;
739  allocateProjectionCoefficients (C2, Q, X, false);
740 
741  // Perform the second projection pass:
742  // C2 = Q' X, X = X - Q*C2
743  rawProject (X, Q, C2);
744  // Update the projection coefficients
745  for (int k = 0; k < num_Q_blocks; ++k)
746  *C[k] += *C2[k];
747  }
748  }
749  }
750 
751 
752 
753  template<class Scalar, class MV>
754  int
756  {
757  using Teuchos::Range1D;
758  using Teuchos::RCP;
759 
760  // MVT returns int for this, even though the "local ordinal
761  // type" of the MV may be some other type (for example,
762  // Tpetra::MultiVector<double, int32_t, int64_t, ...>).
763  const int numCols = MVT::GetNumberVecs (X);
764 
765  // This special case (for X having only one column) makes
766  // TsqrOrthoManagerImpl equivalent to Modified Gram-Schmidt
767  // orthogonalization with conditional full reorthogonalization,
768  // if all multivector inputs have only one column. It also
769  // avoids allocating Q_ scratch space and copying data when we
770  // don't need to invoke TSQR at all.
771  if (numCols == 1) {
772  return normalizeOne (X, B);
773  }
774 
775  // We use Q_ as scratch space for the normalization, since TSQR
776  // requires a scratch multivector (it can't factor in place). Q_
777  // should come from a vector space compatible with X's vector
778  // space, and Q_ should have at least as many columns as X.
779  // Otherwise, we have to reallocate. We also have to allocate
780  // (not "re-") Q_ if we haven't allocated it before. (We can't
781  // allocate Q_ until we have some X, so we need a multivector as
782  // the "prototype.")
783  //
784  // NOTE (mfh 11 Jan 2011) We only increase the number of columsn
785  // in Q_, never decrease. This is OK for typical uses of TSQR,
786  // but you might prefer different behavior in some cases.
787  //
788  // NOTE (mfh 10 Mar 2011) We should only reuse the scratch space
789  // Q_ if X and Q_ have compatible data distributions. However,
790  // Belos' current MultiVecTraits interface does not let us check
791  // for this. Thus, we can only check whether X and Q_ have the
792  // same number of rows. This will behave correctly for the common
793  // case in Belos that all multivectors with the same number of
794  // rows have the same data distribution.
795  //
796  // The specific MV implementation may do more checks than this on
797  // its own and throw an exception if X and Q_ are not compatible,
798  // but it may not. If you find that recycling the Q_ space causes
799  // troubles, you may consider modifying the code below to
800  // reallocate Q_ for every X that comes in.
801  if (Q_.is_null() ||
802  MVT::GetGlobalLength(*Q_) != MVT::GetGlobalLength(X) ||
803  numCols > MVT::GetNumberVecs (*Q_)) {
804  Q_ = MVT::Clone (X, numCols);
805  }
806 
807  // normalizeImpl() wants the second MV argument to have the same
808  // number of columns as X. To ensure this, we pass it a view of
809  // Q_ if Q_ has more columns than X. (This is possible if we
810  // previously called normalize() with a different multivector,
811  // since we never reallocate Q_ if it has more columns than
812  // necessary.)
813  if (MVT::GetNumberVecs(*Q_) == numCols) {
814  return normalizeImpl (X, *Q_, B, false);
815  } else {
816  RCP<MV> Q_view = MVT::CloneViewNonConst (*Q_, Range1D(0, numCols-1));
817  return normalizeImpl (X, *Q_view, B, false);
818  }
819  }
820 
821  template<class Scalar, class MV>
822  void
824  allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C,
825  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
826  const MV& X,
827  const bool attemptToRecycle) const
828  {
829  const int num_Q_blocks = Q.size();
830  const int ncols_X = MVT::GetNumberVecs (X);
831  C.resize (num_Q_blocks);
832  if (attemptToRecycle)
833  {
834  for (int i = 0; i < num_Q_blocks; ++i)
835  {
836  const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
837  // Create a new C[i] if necessary, otherwise resize if
838  // necessary, otherwise fill with zeros.
839  if (C[i].is_null())
840  C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
841  else
842  {
843  mat_type& Ci = *C[i];
844  if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X)
845  Ci.shape (ncols_Qi, ncols_X);
846  else
847  Ci.putScalar (SCT::zero());
848  }
849  }
850  }
851  else
852  {
853  for (int i = 0; i < num_Q_blocks; ++i)
854  {
855  const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
856  C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
857  }
858  }
859  }
860 
861  template<class Scalar, class MV>
862  int
864  normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B)
865  {
866  const int numVecs = MVT::GetNumberVecs(X);
867  if (numVecs == 0) {
868  return 0; // Nothing to do.
869  } else if (numVecs == 1) {
870  // Special case for a single column; scale and copy over.
871  using Teuchos::Range1D;
872  using Teuchos::RCP;
873  using Teuchos::rcp;
874 
875  // Normalize X in place (faster than TSQR for one column).
876  const int rank = normalizeOne (X, B);
877  // Copy results to first column of Q.
878  RCP<MV> Q_0 = MVT::CloneViewNonConst (Q, Range1D(0,0));
879  MVT::Assign (X, *Q_0);
880  return rank;
881  } else {
882  // The "true" argument to normalizeImpl() means the output
883  // vectors go into Q, and the contents of X are overwritten with
884  // invalid values.
885  return normalizeImpl (X, Q, B, true);
886  }
887  }
888 
889  template<class Scalar, class MV>
890  int
892  projectAndNormalizeImpl (MV& X_in,
893  MV& X_out, // Only written if outOfPlace==false.
894  const bool outOfPlace,
895  Teuchos::Array<mat_ptr> C,
896  mat_ptr B,
897  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
898  {
899  using Teuchos::Range1D;
900  using Teuchos::RCP;
901  using Teuchos::rcp;
902 
903  if (outOfPlace) {
904  // Make sure that X_out has at least as many columns as X_in.
905  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(X_out) < MVT::GetNumberVecs(X_in),
906  std::invalid_argument,
907  "Belos::TsqrOrthoManagerImpl::"
908  "projectAndNormalizeOutOfPlace(...):"
909  "X_out has " << MVT::GetNumberVecs(X_out)
910  << " columns, but X_in has "
911  << MVT::GetNumberVecs(X_in) << " columns.");
912  }
913  // Fetch dimensions of X_in and Q, and allocate space for first-
914  // and second-pass projection coefficients (C resp. C2).
915  int ncols_X, num_Q_blocks, ncols_Q_total;
916  checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X_in, Q);
917 
918  // Test for quick exit: if any dimension of X is zero.
919  if (ncols_X == 0) {
920  return 0;
921  }
922  // If there are zero Q blocks or zero Q columns, just normalize!
923  if (num_Q_blocks == 0 || ncols_Q_total == 0) {
924  if (outOfPlace) {
925  return normalizeOutOfPlace (X_in, X_out, B);
926  } else {
927  return normalize (X_in, B);
928  }
929  }
930 
931  // The typical case is that the entries of C have been allocated
932  // before, so we attempt to recycle the allocations. The call
933  // below will reallocate if it cannot recycle.
934  allocateProjectionCoefficients (C, Q, X_in, true);
935 
936  // If we are doing block reorthogonalization, then compute the
937  // column norms of X before projecting for the first time. This
938  // will help us decide whether we need to reorthogonalize X.
939  std::vector<magnitude_type> normsBeforeFirstPass (ncols_X, SCTM::zero());
940  if (reorthogonalizeBlocks_) {
941  MVT::MvNorm (X_in, normsBeforeFirstPass);
942  }
943 
944  // First (Modified) Block Gram-Schmidt pass, in place in X_in.
945  rawProject (X_in, Q, C);
946 
947  // Make space for the normalization coefficients. This will
948  // either be a freshly allocated matrix (if B is null), or a view
949  // of the appropriately sized upper left submatrix of *B (if B is
950  // not null).
951  //
952  // Note that if we let the normalize() routine allocate (in the
953  // case that B is null), that storage will go away at the end of
954  // normalize(). (This is because it passes the RCP by value, not
955  // by reference.)
956  mat_ptr B_out;
957  if (B.is_null()) {
958  B_out = rcp (new mat_type (ncols_X, ncols_X));
959  } else {
960  // Make sure that B is no smaller than numCols x numCols.
961  TEUCHOS_TEST_FOR_EXCEPTION(B->numRows() < ncols_X || B->numCols() < ncols_X,
962  std::invalid_argument,
963  "normalizeOne: Input matrix B must be at "
964  "least " << ncols_X << " x " << ncols_X
965  << ", but is instead " << B->numRows()
966  << " x " << B->numCols() << ".");
967  // Create a view of the ncols_X by ncols_X upper left
968  // submatrix of *B. TSQR will write the normalization
969  // coefficients there.
970  B_out = rcp (new mat_type (Teuchos::View, *B, ncols_X, ncols_X));
971  }
972 
973  // Rank of X(_in) after first projection pass. If outOfPlace,
974  // this overwrites X_in with invalid values, and the results go in
975  // X_out. Otherwise, it's done in place in X_in.
976  const int firstPassRank = outOfPlace ?
977  normalizeOutOfPlace (X_in, X_out, B_out) :
978  normalize (X_in, B_out);
979  if (B.is_null()) {
980  // The input matrix B is null, so assign B_out to it. If B was
981  // not null on input, then B_out is a view of *B, so we don't
982  // have to do anything here. Note that SerialDenseMatrix uses
983  // raw pointers to store data and represent views, so we have to
984  // be careful about scope.
985  B = B_out;
986  }
987  int rank = firstPassRank; // Current rank of X.
988 
989  // If X was not full rank after projection and randomizeNullSpace_
990  // is true, then normalize(OutOfPlace)() replaced the null space
991  // basis of X with random vectors, and orthogonalized them against
992  // the column space basis of X. However, we still need to
993  // orthogonalize the random vectors against the Q[i], after which
994  // we will need to renormalize them.
995  //
996  // If outOfPlace, then we need to work in X_out (where
997  // normalizeOutOfPlace() wrote the normalized vectors).
998  // Otherwise, we need to work in X_in.
999  //
1000  // Note: We don't need to keep the new projection coefficients,
1001  // since they are multiplied by the "small" part of B
1002  // corresponding to the null space of the original X.
1003  if (firstPassRank < ncols_X && randomizeNullSpace_) {
1004  const int numNullSpaceCols = ncols_X - firstPassRank;
1005  const Range1D nullSpaceIndices (firstPassRank, ncols_X - 1);
1006 
1007  // Space for projection coefficients (will be thrown away)
1008  Teuchos::Array<mat_ptr> C_null (num_Q_blocks);
1009  for (int k = 0; k < num_Q_blocks; ++k) {
1010  const int numColsQk = MVT::GetNumberVecs(*Q[k]);
1011  C_null[k] = rcp (new mat_type (numColsQk, numNullSpaceCols));
1012  }
1013  // Space for normalization coefficients (will be thrown away).
1014  RCP<mat_type> B_null (new mat_type (numNullSpaceCols, numNullSpaceCols));
1015 
1016  int randomVectorsRank;
1017  if (outOfPlace) {
1018  // View of the null space basis columns of X.
1019  // normalizeOutOfPlace() wrote them into X_out.
1020  RCP<MV> X_out_null = MVT::CloneViewNonConst (X_out, nullSpaceIndices);
1021  // Use X_in for scratch space. Copy X_out_null into the
1022  // last few columns of X_in (X_in_null) and do projections
1023  // in there. (This saves us a copy wen we renormalize
1024  // (out of place) back into X_out.)
1025  RCP<MV> X_in_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
1026  MVT::Assign (*X_out_null, *X_in_null);
1027  // Project the new random vectors against the Q blocks, and
1028  // renormalize the result into X_out_null.
1029  rawProject (*X_in_null, Q, C_null);
1030  randomVectorsRank = normalizeOutOfPlace (*X_in_null, *X_out_null, B_null);
1031  } else {
1032  // View of the null space columns of X.
1033  // They live in X_in.
1034  RCP<MV> X_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
1035  // Project the new random vectors against the Q blocks,
1036  // and renormalize the result (in place).
1037  rawProject (*X_null, Q, C_null);
1038  randomVectorsRank = normalize (*X_null, B_null);
1039  }
1040  // While unusual, it is still possible for the random data not
1041  // to be full rank after projection and normalization. In that
1042  // case, we could try another set of random data and recurse as
1043  // necessary, but instead for now we just raise an exception.
1044  TEUCHOS_TEST_FOR_EXCEPTION(randomVectorsRank != numNullSpaceCols,
1045  TsqrOrthoError,
1046  "Belos::TsqrOrthoManagerImpl::projectAndNormalize"
1047  "OutOfPlace(): After projecting and normalizing the "
1048  "random vectors (used to replace the null space "
1049  "basis vectors from normalizing X), they have rank "
1050  << randomVectorsRank << ", but should have full "
1051  "rank " << numNullSpaceCols << ".");
1052  }
1053 
1054  // Whether or not X_in was full rank after projection, we still
1055  // might want to reorthogonalize against Q.
1056  if (reorthogonalizeBlocks_) {
1057  // We are only interested in the column space basis of X
1058  // resp. X_out.
1059  std::vector<magnitude_type>
1060  normsAfterFirstPass (firstPassRank, SCTM::zero());
1061  std::vector<magnitude_type>
1062  normsAfterSecondPass (firstPassRank, SCTM::zero());
1063 
1064  // Compute post-first-pass (pre-normalization) norms. We could
1065  // have done that using MVT::MvNorm() on X_in after projecting,
1066  // but before the first normalization. However, that operation
1067  // may be expensive. It is also unnecessary: after calling
1068  // normalize(), the 2-norm of B(:,j) is the 2-norm of X_in(:,j)
1069  // before normalization, in exact arithmetic.
1070  //
1071  // NOTE (mfh 06 Nov 2010) This is one way that combining
1072  // projection and normalization into a single kernel --
1073  // projectAndNormalize() -- pays off. In project(), we have to
1074  // compute column norms of X before and after projection. Here,
1075  // we get them for free from the normalization coefficients.
1076  Teuchos::BLAS<int, Scalar> blas;
1077  for (int j = 0; j < firstPassRank; ++j) {
1078  const Scalar* const B_j = &(*B_out)(0,j);
1079  // Teuchos::BLAS::NRM2 returns a magnitude_type result on
1080  // Scalar inputs.
1081  normsAfterFirstPass[j] = blas.NRM2 (firstPassRank, B_j, 1);
1082  }
1083  // Test whether any of the norms dropped below the
1084  // reorthogonalization threshold.
1085  bool reorthogonalize = false;
1086  for (int j = 0; j < firstPassRank; ++j) {
1087  // If any column's norm decreased too much, mark this block
1088  // for reorthogonalization. Note that this test will _not_
1089  // activate reorthogonalization if a column's norm before the
1090  // first project-and-normalize step was zero. It _will_
1091  // activate reorthogonalization if the column's norm before
1092  // was not zero, but is zero now.
1093  const magnitude_type curThreshold =
1094  blockReorthogThreshold() * normsBeforeFirstPass[j];
1095  if (normsAfterFirstPass[j] < curThreshold) {
1096  reorthogonalize = true;
1097  break;
1098  }
1099  }
1100  // Perform another Block Gram-Schmidt pass if necessary. "Twice
1101  // is enough" (Kahan's theorem) for a Krylov method, unless
1102  // (using Stewart's term) there is an "orthogonalization fault"
1103  // (indicated by reorthogFault).
1104  //
1105  // NOTE (mfh 07 Nov 2010) For now, we include the entire block
1106  // of X, including possible random data (that was already
1107  // projected and normalized above). It might make more sense
1108  // just to process the first firstPassRank columns of X.
1109  // However, the resulting reorthogonalization should still be
1110  // correct regardless.
1111  bool reorthogFault = false;
1112  // Indices of X at which there was an orthogonalization fault.
1113  std::vector<int> faultIndices;
1114  if (reorthogonalize) {
1115  using Teuchos::Copy;
1116  using Teuchos::NO_TRANS;
1117 
1118  // If we're using out-of-place normalization, copy X_out
1119  // (results of first project and normalize pass) back into
1120  // X_in, for the second project and normalize pass.
1121  if (outOfPlace)
1122  MVT::Assign (X_out, X_in);
1123 
1124  // C2 is only used internally, so we know that we are
1125  // allocating fresh and not recycling allocations. Stating
1126  // this lets us save time checking dimensions.
1127  Teuchos::Array<mat_ptr> C2;
1128  allocateProjectionCoefficients (C2, Q, X_in, false);
1129 
1130  // Block Gram-Schmidt (again). Delay updating the block
1131  // coefficients until we have the new normalization
1132  // coefficients, which we need in order to do the update.
1133  rawProject (X_in, Q, C2);
1134 
1135  // Coefficients for (re)normalization of X_in.
1136  RCP<mat_type> B2 (new mat_type (ncols_X, ncols_X));
1137 
1138  // Normalize X_in (into X_out, if working out of place).
1139  const int secondPassRank = outOfPlace ?
1140  normalizeOutOfPlace (X_in, X_out, B2) :
1141  normalize (X_in, B2);
1142  rank = secondPassRank; // Current rank of X
1143 
1144  // Update normalization coefficients. We begin with copying
1145  // B_out, since the BLAS' _GEMM routine doesn't let us alias
1146  // its input and output arguments.
1147  mat_type B_copy (Copy, *B_out, B_out->numRows(), B_out->numCols());
1148  // B_out := B2 * B_out (where input B_out is in B_copy).
1149  const int err = B_out->multiply (NO_TRANS, NO_TRANS, SCT::one(),
1150  *B2, B_copy, SCT::zero());
1151  TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::logic_error,
1152  "Teuchos::SerialDenseMatrix::multiply "
1153  "returned err = " << err << " != 0");
1154  // Update the block coefficients from the projection step. We
1155  // use B_copy for this (a copy of B_out, the first-pass
1156  // normalization coefficients).
1157  for (int k = 0; k < num_Q_blocks; ++k) {
1158  mat_type C_k_copy (Copy, *C[k], C[k]->numRows(), C[k]->numCols());
1159 
1160  // C[k] := C2[k]*B_copy + C[k].
1161  const int err2 = C[k]->multiply (NO_TRANS, NO_TRANS, SCT::one(),
1162  *C2[k], B_copy, SCT::one());
1163  TEUCHOS_TEST_FOR_EXCEPTION(err2 != 0, std::logic_error,
1164  "Teuchos::SerialDenseMatrix::multiply "
1165  "returned err = " << err << " != 0");
1166  }
1167  // Compute post-second-pass (pre-normalization) norms, using
1168  // B2 (the coefficients from the second normalization step) in
1169  // the same way as with B_out before.
1170  for (int j = 0; j < rank; ++j) {
1171  const Scalar* const B2_j = &(*B2)(0,j);
1172  normsAfterSecondPass[j] = blas.NRM2 (rank, B2_j, 1);
1173  }
1174  // Test whether any of the norms dropped below the
1175  // reorthogonalization threshold. If so, it's an
1176  // orthogonalization fault, which requires expensive recovery.
1177  reorthogFault = false;
1178  for (int j = 0; j < rank; ++j) {
1179  const magnitude_type relativeLowerBound =
1180  blockReorthogThreshold() * normsAfterFirstPass[j];
1181  if (normsAfterSecondPass[j] < relativeLowerBound) {
1182  reorthogFault = true;
1183  faultIndices.push_back (j);
1184  }
1185  }
1186  } // if (reorthogonalize) // reorthogonalization pass
1187 
1188  if (reorthogFault) {
1189  if (throwOnReorthogFault_) {
1190  raiseReorthogFault (normsAfterFirstPass,
1191  normsAfterSecondPass,
1192  faultIndices);
1193  } else {
1194  // NOTE (mfh 19 Jan 2011) We could handle the fault here by
1195  // slowly reorthogonalizing, one vector at a time, the
1196  // offending vectors of X. However, we choose not to
1197  // implement this for now. If it becomes a problem, let us
1198  // know and we will prioritize implementing this.
1199  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
1200  "TsqrOrthoManagerImpl has not yet implemented"
1201  " recovery from an orthogonalization fault.");
1202  }
1203  }
1204  } // if (reorthogonalizeBlocks_)
1205  return rank;
1206  }
1207 
1208 
1209  template<class Scalar, class MV>
1210  void
1211  TsqrOrthoManagerImpl<Scalar, MV>::
1212  raiseReorthogFault (const std::vector<magnitude_type>& normsAfterFirstPass,
1213  const std::vector<magnitude_type>& normsAfterSecondPass,
1214  const std::vector<int>& faultIndices)
1215  {
1216  using std::endl;
1217  typedef std::vector<int>::size_type size_type;
1218  std::ostringstream os;
1219 
1220  os << "Orthogonalization fault at the following column(s) of X:" << endl;
1221  os << "Column\tNorm decrease factor" << endl;
1222  for (size_type k = 0; k < faultIndices.size(); ++k)
1223  {
1224  const int index = faultIndices[k];
1225  const magnitude_type decreaseFactor =
1226  normsAfterSecondPass[index] / normsAfterFirstPass[index];
1227  os << index << "\t" << decreaseFactor << endl;
1228  }
1229  throw TsqrOrthoFault (os.str());
1230  }
1231 
1232  template<class Scalar, class MV>
1233  Teuchos::RCP<const Teuchos::ParameterList>
1235  {
1236  using Teuchos::ParameterList;
1237  using Teuchos::parameterList;
1238  using Teuchos::RCP;
1239 
1240  if (defaultParams_.is_null()) {
1241  RCP<ParameterList> params = parameterList ("TsqrOrthoManagerImpl");
1242  //
1243  // TSQR parameters (set as a sublist).
1244  //
1245  params->set ("TSQR implementation", *(tsqrAdaptor_.getValidParameters()),
1246  "TSQR implementation parameters.");
1247  //
1248  // Orthogonalization parameters
1249  //
1250  const bool defaultRandomizeNullSpace = true;
1251  params->set ("randomizeNullSpace", defaultRandomizeNullSpace,
1252  "Whether to fill in null space vectors with random data.");
1253 
1254  const bool defaultReorthogonalizeBlocks = true;
1255  params->set ("reorthogonalizeBlocks", defaultReorthogonalizeBlocks,
1256  "Whether to do block reorthogonalization as necessary.");
1257 
1258  // This parameter corresponds to the "blk_tol_" parameter in
1259  // Belos' DGKSOrthoManager. We choose the same default value.
1260  const magnitude_type defaultBlockReorthogThreshold =
1261  magnitude_type(10) * SCTM::squareroot (SCTM::eps());
1262  params->set ("blockReorthogThreshold", defaultBlockReorthogThreshold,
1263  "If reorthogonalizeBlocks==true, and if the norm of "
1264  "any column within a block decreases by this much or "
1265  "more after orthogonalization, we reorthogonalize.");
1266 
1267  // This parameter corresponds to the "sing_tol_" parameter in
1268  // Belos' DGKSOrthoManager. We choose the same default value.
1269  const magnitude_type defaultRelativeRankTolerance =
1270  Teuchos::as<magnitude_type>(10) * SCTM::eps();
1271 
1272  // If the relative rank tolerance is zero, then we will always
1273  // declare blocks to be numerically full rank, as long as no
1274  // singular values are zero.
1275  params->set ("relativeRankTolerance", defaultRelativeRankTolerance,
1276  "Relative tolerance to determine the numerical rank of a "
1277  "block when normalizing.");
1278 
1279  // See Stewart's 2008 paper on block Gram-Schmidt for a definition
1280  // of "orthogonalization fault."
1281  const bool defaultThrowOnReorthogFault = true;
1282  params->set ("throwOnReorthogFault", defaultThrowOnReorthogFault,
1283  "Whether to throw an exception if an orthogonalization "
1284  "fault occurs. This only matters if reorthogonalization "
1285  "is enabled (reorthogonalizeBlocks==true).");
1286 
1287  const bool defaultForceNonnegativeDiagonal = false;
1288  params->set ("forceNonnegativeDiagonal", defaultForceNonnegativeDiagonal,
1289  "Whether to force the R factor produced by the normalization "
1290  "step to have a nonnegative diagonal.");
1291 
1292  defaultParams_ = params;
1293  }
1294  return defaultParams_;
1295  }
1296 
1297  template<class Scalar, class MV>
1298  Teuchos::RCP<const Teuchos::ParameterList>
1300  {
1301  using Teuchos::ParameterList;
1302  using Teuchos::RCP;
1303  using Teuchos::rcp;
1304 
1305  RCP<const ParameterList> defaultParams = getValidParameters();
1306  // Start with a clone of the default parameters.
1307  RCP<ParameterList> params = rcp (new ParameterList (*defaultParams));
1308 
1309  // Disable reorthogonalization and randomization of the null
1310  // space basis. Reorthogonalization tolerances don't matter,
1311  // since we aren't reorthogonalizing blocks in the fast
1312  // settings. We can leave the default values. Also,
1313  // (re)orthogonalization faults may only occur with
1314  // reorthogonalization, so we don't have to worry about the
1315  // "throwOnReorthogFault" setting.
1316  const bool randomizeNullSpace = false;
1317  params->set ("randomizeNullSpace", randomizeNullSpace);
1318  const bool reorthogonalizeBlocks = false;
1319  params->set ("reorthogonalizeBlocks", reorthogonalizeBlocks);
1320 
1321  return params;
1322  }
1323 
1324  template<class Scalar, class MV>
1325  int
1327  rawNormalize (MV& X,
1328  MV& Q,
1329  Teuchos::SerialDenseMatrix<int, Scalar>& B)
1330  {
1331  int rank;
1332  try {
1333  // This call only computes the QR factorization X = Q B.
1334  // It doesn't compute the rank of X. That comes from
1335  // revealRank() below.
1336  tsqrAdaptor_.factorExplicit (X, Q, B, forceNonnegativeDiagonal_);
1337  // This call will only modify *B if *B on input is not of full
1338  // numerical rank.
1339  rank = tsqrAdaptor_.revealRank (Q, B, relativeRankTolerance_);
1340  } catch (std::exception& e) {
1341  throw TsqrOrthoError (e.what()); // Toss the exception up the chain.
1342  }
1343  return rank;
1344  }
1345 
1346  template<class Scalar, class MV>
1347  int
1348  TsqrOrthoManagerImpl<Scalar, MV>::
1349  normalizeOne (MV& X,
1350  Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > B) const
1351  {
1352  // Make space for the normalization coefficient. This will either
1353  // be a freshly allocated matrix (if B is null), or a view of the
1354  // 1x1 upper left submatrix of *B (if B is not null).
1355  mat_ptr B_out;
1356  if (B.is_null()) {
1357  B_out = rcp (new mat_type (1, 1));
1358  } else {
1359  const int numRows = B->numRows();
1360  const int numCols = B->numCols();
1361  TEUCHOS_TEST_FOR_EXCEPTION(numRows < 1 || numCols < 1,
1362  std::invalid_argument,
1363  "normalizeOne: Input matrix B must be at "
1364  "least 1 x 1, but is instead " << numRows
1365  << " x " << numCols << ".");
1366  // Create a view of the 1x1 upper left submatrix of *B.
1367  B_out = rcp (new mat_type (Teuchos::View, *B, 1, 1));
1368  }
1369 
1370  // Compute the norm of X, and write the result to B_out.
1371  std::vector<magnitude_type> theNorm (1, SCTM::zero());
1372  MVT::MvNorm (X, theNorm);
1373  (*B_out)(0,0) = theNorm[0];
1374 
1375  if (B.is_null()) {
1376  // The input matrix B is null, so assign B_out to it. If B was
1377  // not null on input, then B_out is a view of *B, so we don't
1378  // have to do anything here. Note that SerialDenseMatrix uses
1379  // raw pointers to store data and represent views, so we have to
1380  // be careful about scope.
1381  B = B_out;
1382  }
1383 
1384  // Scale X by its norm, if its norm is zero. Otherwise, do the
1385  // right thing based on whether the user wants us to fill the null
1386  // space with random vectors.
1387  if (theNorm[0] == SCTM::zero()) {
1388  // Make a view of the first column of Q, fill it with random
1389  // data, and normalize it. Throw away the resulting norm.
1390  if (randomizeNullSpace_) {
1391  MVT::MvRandom(X);
1392  MVT::MvNorm (X, theNorm);
1393  if (theNorm[0] == SCTM::zero()) {
1394  // It is possible that a random vector could have all zero
1395  // entries, but unlikely. We could try again, but it's also
1396  // possible that multiple tries could result in zero
1397  // vectors. We choose instead to give up.
1398  throw TsqrOrthoError("normalizeOne: a supposedly random "
1399  "vector has norm zero!");
1400  } else {
1401  // NOTE (mfh 09 Nov 2010) I'm assuming that dividing a
1402  // Scalar by a magnitude_type is defined and that it results
1403  // in a Scalar.
1404  const Scalar alpha = SCT::one() / theNorm[0];
1405  MVT::MvScale (X, alpha);
1406  }
1407  }
1408  return 0; // The rank of the matrix (actually just one vector) X.
1409  } else {
1410  // NOTE (mfh 09 Nov 2010) I'm assuming that dividing a Scalar by
1411  // a magnitude_type is defined and that it results in a Scalar.
1412  const Scalar alpha = SCT::one() / theNorm[0];
1413  MVT::MvScale (X, alpha);
1414  return 1; // The rank of the matrix (actually just one vector) X.
1415  }
1416  }
1417 
1418 
1419  template<class Scalar, class MV>
1420  void
1421  TsqrOrthoManagerImpl<Scalar, MV>::
1422  rawProject (MV& X,
1423  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
1424  Teuchos::ArrayView<Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > > C) const
1425  {
1426  // "Modified Gram-Schmidt" version of Block Gram-Schmidt.
1427  const int num_Q_blocks = Q.size();
1428  for (int i = 0; i < num_Q_blocks; ++i)
1429  {
1430  // TEUCHOS_TEST_FOR_EXCEPTION(C[i].is_null(), std::logic_error,
1431  // "TsqrOrthoManagerImpl::rawProject(): C["
1432  // << i << "] is null");
1433  // TEUCHOS_TEST_FOR_EXCEPTION(Q[i].is_null(), std::logic_error,
1434  // "TsqrOrthoManagerImpl::rawProject(): Q["
1435  // << i << "] is null");
1436  mat_type& Ci = *C[i];
1437  const MV& Qi = *Q[i];
1438  innerProd (Qi, X, Ci);
1439  MVT::MvTimesMatAddMv (-SCT::one(), Qi, Ci, SCT::one(), X);
1440  }
1441  }
1442 
1443 
1444  template<class Scalar, class MV>
1445  void
1446  TsqrOrthoManagerImpl<Scalar, MV>::
1447  rawProject (MV& X,
1448  const Teuchos::RCP<const MV>& Q,
1449  const Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> >& C) const
1450  {
1451  // Block Gram-Schmidt
1452  innerProd (*Q, X, *C);
1453  MVT::MvTimesMatAddMv (-SCT::one(), *Q, *C, SCT::one(), X);
1454  }
1455 
1456  template<class Scalar, class MV>
1457  int
1458  TsqrOrthoManagerImpl<Scalar, MV>::
1459  normalizeImpl (MV& X,
1460  MV& Q,
1461  Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > B,
1462  const bool outOfPlace)
1463  {
1464  using Teuchos::Range1D;
1465  using Teuchos::RCP;
1466  using Teuchos::rcp;
1467  using Teuchos::ScalarTraits;
1468  using Teuchos::tuple;
1469  using std::cerr;
1470  using std::endl;
1471  // Don't set this to true unless you want lots of debugging
1472  // messages written to stderr on every MPI process.
1473  const bool extraDebug = false;
1474 
1475  const int numCols = MVT::GetNumberVecs (X);
1476  if (numCols == 0) {
1477  return 0; // Fast exit for an empty input matrix.
1478  }
1479 
1480  // We allow Q to have more columns than X. In that case, we only
1481  // touch the first numCols columns of Q.
1482  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(Q) < numCols,
1483  std::invalid_argument,
1484  "TsqrOrthoManagerImpl::normalizeImpl(X,Q,B): Q has "
1485  << MVT::GetNumberVecs(Q) << " columns. This is too "
1486  "few, since X has " << numCols << " columns.");
1487  // TSQR wants a Q with the same number of columns as X, so have it
1488  // work on a nonconstant view of Q with the same number of columns
1489  // as X.
1490  RCP<MV> Q_view = MVT::CloneViewNonConst (Q, Range1D(0, numCols-1));
1491 
1492  // Make space for the normalization coefficients. This will
1493  // either be a freshly allocated matrix (if B is null), or a view
1494  // of the appropriately sized upper left submatrix of *B (if B is
1495  // not null).
1496  mat_ptr B_out;
1497  if (B.is_null()) {
1498  B_out = rcp (new mat_type (numCols, numCols));
1499  } else {
1500  // Make sure that B is no smaller than numCols x numCols.
1501  TEUCHOS_TEST_FOR_EXCEPTION(B->numRows() < numCols || B->numCols() < numCols,
1502  std::invalid_argument,
1503  "normalizeOne: Input matrix B must be at "
1504  "least " << numCols << " x " << numCols
1505  << ", but is instead " << B->numRows()
1506  << " x " << B->numCols() << ".");
1507  // Create a view of the numCols x numCols upper left submatrix
1508  // of *B. TSQR will write the normalization coefficients there.
1509  B_out = rcp (new mat_type (Teuchos::View, *B, numCols, numCols));
1510  }
1511 
1512  if (extraDebug) {
1513  std::vector<magnitude_type> norms (numCols);
1514  MVT::MvNorm (X, norms);
1515  cerr << "Column norms of X before orthogonalization: ";
1516  typedef typename std::vector<magnitude_type>::const_iterator iter_type;
1517  for (iter_type iter = norms.begin(); iter != norms.end(); ++iter) {
1518  cerr << *iter;
1519  if (iter+1 != norms.end())
1520  cerr << ", ";
1521  }
1522  }
1523 
1524  // Compute rank-revealing decomposition (in this case, TSQR of X
1525  // followed by SVD of the R factor and appropriate updating of the
1526  // resulting Q factor) of X. X is modified in place and filled
1527  // with garbage, and Q_view contains the resulting explicit Q
1528  // factor. Later, we will copy this back into X.
1529  //
1530  // The matrix *B_out will only be upper triangular if X is of full
1531  // numerical rank. Otherwise, the entries below the diagonal may
1532  // be filled in as well.
1533  const int rank = rawNormalize (X, *Q_view, *B_out);
1534  if (B.is_null()) {
1535  // The input matrix B is null, so assign B_out to it. If B was
1536  // not null on input, then B_out is a view of *B, so we don't
1537  // have to do anything here. Note that SerialDenseMatrix uses
1538  // raw pointers to store data and represent views, so we have to
1539  // be careful about scope.
1540  B = B_out;
1541  }
1542 
1543  if (extraDebug) {
1544  std::vector<magnitude_type> norms (numCols);
1545  MVT::MvNorm (*Q_view, norms);
1546  cerr << "Column norms of Q_view after orthogonalization: ";
1547  for (typename std::vector<magnitude_type>::const_iterator iter = norms.begin();
1548  iter != norms.end(); ++iter) {
1549  cerr << *iter;
1550  if (iter+1 != norms.end())
1551  cerr << ", ";
1552  }
1553  }
1554  TEUCHOS_TEST_FOR_EXCEPTION(rank < 0 || rank > numCols, std::logic_error,
1555  "Belos::TsqrOrthoManagerImpl::normalizeImpl: "
1556  "rawNormalize() returned rank = " << rank << " for a "
1557  "matrix X with " << numCols << " columns. Please report"
1558  " this bug to the Belos developers.");
1559  if (extraDebug && rank == 0) {
1560  // Sanity check: ensure that the columns of X are sufficiently
1561  // small for X to be reported as rank zero.
1562  const mat_type& B_ref = *B;
1563  std::vector<magnitude_type> norms (B_ref.numCols());
1564  for (typename mat_type::ordinalType j = 0; j < B_ref.numCols(); ++j) {
1565  typedef typename mat_type::scalarType mat_scalar_type;
1566  mat_scalar_type sumOfSquares = ScalarTraits<mat_scalar_type>::zero();
1567  for (typename mat_type::ordinalType i = 0; i <= j; ++i) {
1568  const mat_scalar_type B_ij =
1569  ScalarTraits<mat_scalar_type>::magnitude (B_ref(i,j));
1570  sumOfSquares += B_ij*B_ij;
1571  }
1572  norms[j] = ScalarTraits<mat_scalar_type>::squareroot (sumOfSquares);
1573  }
1574  using std::cerr;
1575  using std::endl;
1576  cerr << "Norms of columns of B after orthogonalization: ";
1577  for (typename mat_type::ordinalType j = 0; j < B_ref.numCols(); ++j) {
1578  cerr << norms[j];
1579  if (j != B_ref.numCols() - 1)
1580  cerr << ", ";
1581  }
1582  cerr << endl;
1583  }
1584 
1585  // If X is full rank or we don't want to replace its null space
1586  // basis with random vectors, then we're done.
1587  if (rank == numCols || ! randomizeNullSpace_) {
1588  // If we're supposed to be working in place in X, copy the
1589  // results back from Q_view into X.
1590  if (! outOfPlace) {
1591  MVT::Assign (*Q_view, X);
1592  }
1593  return rank;
1594  }
1595 
1596  if (randomizeNullSpace_ && rank < numCols) {
1597  // X wasn't full rank. Replace the null space basis of X (in
1598  // the last numCols-rank columns of Q_view) with random data,
1599  // project it against the first rank columns of Q_view, and
1600  // normalize.
1601  //
1602  // Number of columns to fill with random data.
1603  const int nullSpaceNumCols = numCols - rank;
1604  // Inclusive range of indices of columns of X to fill with
1605  // random data.
1606  Range1D nullSpaceIndices (rank, numCols-1);
1607 
1608  // rawNormalize() wrote the null space basis vectors into
1609  // Q_view. We continue to work in place in Q_view by writing
1610  // the random data there and (if there is a nontrival column
1611  // space) projecting in place against the column space basis
1612  // vectors (also in Q_view).
1613  RCP<MV> Q_null = MVT::CloneViewNonConst (*Q_view, nullSpaceIndices);
1614  // Replace the null space basis with random data.
1615  MVT::MvRandom (*Q_null);
1616 
1617  // Make sure that the "random" data isn't all zeros. This is
1618  // statistically nearly impossible, but we test for debugging
1619  // purposes.
1620  {
1621  std::vector<magnitude_type> norms (MVT::GetNumberVecs(*Q_null));
1622  MVT::MvNorm (*Q_null, norms);
1623 
1624  bool anyZero = false;
1625  typedef typename std::vector<magnitude_type>::const_iterator iter_type;
1626  for (iter_type it = norms.begin(); it != norms.end(); ++it) {
1627  if (*it == SCTM::zero()) {
1628  anyZero = true;
1629  }
1630  }
1631  if (anyZero) {
1632  std::ostringstream os;
1633  os << "TsqrOrthoManagerImpl::normalizeImpl: "
1634  "We are being asked to randomize the null space, for a matrix "
1635  "with " << numCols << " columns and reported column rank "
1636  << rank << ". The inclusive range of columns to fill with "
1637  "random data is [" << nullSpaceIndices.lbound() << ","
1638  << nullSpaceIndices.ubound() << "]. After filling the null "
1639  "space vectors with random numbers, at least one of the vectors"
1640  " has norm zero. Here are the norms of all the null space "
1641  "vectors: [";
1642  for (iter_type it = norms.begin(); it != norms.end(); ++it) {
1643  os << *it;
1644  if (it+1 != norms.end())
1645  os << ", ";
1646  }
1647  os << "].) There is a tiny probability that this could happen "
1648  "randomly, but it is likely a bug. Please report it to the "
1649  "Belos developers, especially if you are able to reproduce the "
1650  "behavior.";
1651  TEUCHOS_TEST_FOR_EXCEPTION(anyZero, TsqrOrthoError, os.str());
1652  }
1653  }
1654 
1655  if (rank > 0) {
1656  // Project the random data against the column space basis of
1657  // X, using a simple block projection ("Block Classical
1658  // Gram-Schmidt"). This is accurate because we've already
1659  // orthogonalized the column space basis of X nearly to
1660  // machine precision via a QR factorization (TSQR) with
1661  // accuracy comparable to Householder QR.
1662  RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D(0, rank-1));
1663 
1664  // Temporary storage for projection coefficients. We don't
1665  // need to keep them, since they represent the null space
1666  // basis (for which the coefficients are logically zero).
1667  mat_ptr C_null (new mat_type (rank, nullSpaceNumCols));
1668  rawProject (*Q_null, Q_col, C_null);
1669  }
1670  // Normalize the projected random vectors, so that they are
1671  // mutually orthogonal (as well as orthogonal to the column
1672  // space basis of X). We use X for the output of the
1673  // normalization: for out-of-place normalization (outOfPlace ==
1674  // true), X is overwritten with "invalid values" anyway, and for
1675  // in-place normalization (outOfPlace == false), we want the
1676  // result to be in X anyway.
1677  RCP<MV> X_null = MVT::CloneViewNonConst (X, nullSpaceIndices);
1678  // Normalization coefficients for projected random vectors.
1679  // Will be thrown away.
1680  mat_type B_null (nullSpaceNumCols, nullSpaceNumCols);
1681  // Write the normalized vectors to X_null (in X).
1682  const int nullSpaceBasisRank = rawNormalize (*Q_null, *X_null, B_null);
1683 
1684  // It's possible, but unlikely, that X_null doesn't have full
1685  // rank (after the projection step). We could recursively fill
1686  // in more random vectors until we finally get a full rank
1687  // matrix, but instead we just throw an exception.
1688  //
1689  // NOTE (mfh 08 Nov 2010) Perhaps we should deal with this case
1690  // more elegantly. Recursion might be one way to solve it, but
1691  // be sure to check that the recursion will terminate. We could
1692  // do this by supplying an additional argument to rawNormalize,
1693  // which is the null space basis rank from the previous
1694  // iteration. The rank has to decrease each time, or the
1695  // recursion may go on forever.
1696  if (nullSpaceBasisRank < nullSpaceNumCols) {
1697  std::vector<magnitude_type> norms (MVT::GetNumberVecs(*X_null));
1698  MVT::MvNorm (*X_null, norms);
1699  std::ostringstream os;
1700  os << "TsqrOrthoManagerImpl::normalizeImpl: "
1701  << "We are being asked to randomize the null space, "
1702  << "for a matrix with " << numCols << " columns and "
1703  << "column rank " << rank << ". After projecting and "
1704  << "normalizing the generated random vectors, they "
1705  << "only have rank " << nullSpaceBasisRank << ". They"
1706  << " should have full rank " << nullSpaceNumCols
1707  << ". (The inclusive range of columns to fill with "
1708  << "random data is [" << nullSpaceIndices.lbound()
1709  << "," << nullSpaceIndices.ubound() << "]. The "
1710  << "column norms of the resulting Q factor are: [";
1711  for (typename std::vector<magnitude_type>::size_type k = 0;
1712  k < norms.size(); ++k) {
1713  os << norms[k];
1714  if (k != norms.size()-1) {
1715  os << ", ";
1716  }
1717  }
1718  os << "].) There is a tiny probability that this could "
1719  << "happen randomly, but it is likely a bug. Please "
1720  << "report it to the Belos developers, especially if "
1721  << "you are able to reproduce the behavior.";
1722 
1723  TEUCHOS_TEST_FOR_EXCEPTION(nullSpaceBasisRank < nullSpaceNumCols,
1724  TsqrOrthoError, os.str());
1725  }
1726  // If we're normalizing out of place, copy the X_null
1727  // vectors back into Q_null; the Q_col vectors are already
1728  // where they are supposed to be in that case.
1729  //
1730  // If we're normalizing in place, leave X_null alone (it's
1731  // already where it needs to be, in X), but copy Q_col back
1732  // into the first rank columns of X.
1733  if (outOfPlace) {
1734  MVT::Assign (*X_null, *Q_null);
1735  } else if (rank > 0) {
1736  // MVT::Assign() doesn't accept empty ranges of columns.
1737  RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D(0, rank-1));
1738  RCP<MV> X_col = MVT::CloneViewNonConst (X, Range1D(0, rank-1));
1739  MVT::Assign (*Q_col, *X_col);
1740  }
1741  }
1742  return rank;
1743  }
1744 
1745 
1746  template<class Scalar, class MV>
1747  void
1748  TsqrOrthoManagerImpl<Scalar, MV>::
1749  checkProjectionDims (int& ncols_X,
1750  int& num_Q_blocks,
1751  int& ncols_Q_total,
1752  const MV& X,
1753  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1754  {
1755  // First assign to temporary values, so the function won't
1756  // commit any externally visible side effects unless it will
1757  // return normally (without throwing an exception). (I'm being
1758  // cautious; MVT::GetNumberVecs() probably shouldn't have any
1759  // externally visible side effects, unless it is logging to a
1760  // file or something.)
1761  int the_ncols_X, the_num_Q_blocks, the_ncols_Q_total;
1762  the_num_Q_blocks = Q.size();
1763  the_ncols_X = MVT::GetNumberVecs (X);
1764 
1765  // Compute the total number of columns of all the Q[i] blocks.
1766  the_ncols_Q_total = 0;
1767  // You should be angry if your compiler doesn't support type
1768  // inference ("auto"). That's why I need this awful typedef.
1769  using Teuchos::ArrayView;
1770  using Teuchos::RCP;
1771  typedef typename ArrayView<RCP<const MV> >::const_iterator iter_type;
1772  for (iter_type it = Q.begin(); it != Q.end(); ++it) {
1773  const MV& Qi = **it;
1774  the_ncols_Q_total += MVT::GetNumberVecs (Qi);
1775  }
1776 
1777  // Commit temporary values to the output arguments.
1778  ncols_X = the_ncols_X;
1779  num_Q_blocks = the_num_Q_blocks;
1780  ncols_Q_total = the_ncols_Q_total;
1781  }
1782 
1783 } // namespace Anasazi
1784 
1785 #endif // __AnasaziTsqrOrthoManagerImpl_hpp
void norm(const MV &X, std::vector< magnitude_type > &normvec) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Default valid parameter list.
TSQR-based OrthoManager subclass implementation.
magnitude_type orthonormError(const MV &X) const
Return .
magnitude_type orthogError(const MV &X1, const MV &X2) const
Return the Frobenius norm of the inner product of X1 with itself.
Declaration of basic traits for the multivector type.
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.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &params)
Set parameters from the given parameter list.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
TsqrOrthoManagerImpl(const Teuchos::RCP< Teuchos::ParameterList > &params, const std::string &label)
Constructor (that sets user-specified parameters).
TsqrOrthoManager(Impl) error.
const std::string & getLabel() const
Get the label for timers (if timers are enabled).
void setLabel(const std::string &label)
Set the label for timers.
Templated virtual class for providing orthogonalization/orthonormalization methods.
Traits class which defines basic operations on multivectors.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters()
Get "fast" parameters for TsqrOrthoManagerImpl.
void project(MV &X, Teuchos::Array< mat_ptr > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Compute and .
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Exception thrown to signal error in an orthogonalization manager method.
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.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
int normalize(MV &X, mat_ptr B)
Orthogonalize the columns of X in place.