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