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

Generated for Belos by doxygen 1.10.0