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