45 #ifndef __BelosTsqrOrthoManagerImpl_hpp 46 #define __BelosTsqrOrthoManagerImpl_hpp 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 129 template<
class Scalar>
131 public std::binary_function<Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>,
132 Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>,
152 operator() (Teuchos::ArrayView<magnitude_type> normsBeforeFirstPass,
153 Teuchos::ArrayView<magnitude_type> normsAfterFirstPass) = 0;
156 template<
class Scalar>
190 template<
class Scalar,
class MV>
192 public Teuchos::ParameterListAcceptorDefaultBase {
198 typedef Teuchos::SerialDenseMatrix<int, Scalar>
mat_type;
202 typedef Teuchos::ScalarTraits<Scalar>
SCT;
203 typedef Teuchos::ScalarTraits<magnitude_type>
SCTM;
249 const std::string& label);
293 #ifdef BELOS_TEUCHOS_TIME_MONITOR 294 clearTimer (label,
"All orthogonalization");
295 clearTimer (label,
"Projection");
296 clearTimer (label,
"Normalization");
298 timerOrtho_ = makeTimer (label,
"All orthogonalization");
299 timerProject_ = makeTimer (label,
"Projection");
300 timerNormalize_ = makeTimer (label,
"Normalization");
301 #endif // BELOS_TEUCHOS_TIME_MONITOR 340 norm (
const MV& X, std::vector<magnitude_type>& normVec)
const;
353 Teuchos::Array<mat_ptr> C,
354 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q);
407 Teuchos::Array<mat_ptr> C,
409 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
438 Teuchos::Array<mat_ptr> C,
440 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
454 const Scalar ONE = SCT::one();
458 for (
int k = 0; k < ncols; ++k) {
461 return XTX.normFrobenius();
471 mat_type X1_T_X2 (ncols_X1, ncols_X2);
473 return X1_T_X2.normFrobenius();
543 #ifdef BELOS_TEUCHOS_TIME_MONITOR 544 Teuchos::RCP<Teuchos::Time> timerOrtho_;
548 Teuchos::RCP<Teuchos::Time> timerProject_;
551 Teuchos::RCP<Teuchos::Time> timerNormalize_;
552 #endif // BELOS_TEUCHOS_TIME_MONITOR 557 #ifdef BELOS_TEUCHOS_TIME_MONITOR 558 static Teuchos::RCP<Teuchos::Time>
567 makeTimer (
const std::string& prefix,
568 const std::string& timerName)
570 const std::string timerLabel =
571 prefix.empty() ? timerName : (prefix +
": " + timerName);
572 return Teuchos::TimeMonitor::getNewCounter (timerLabel);
581 clearTimer (
const std::string& prefix,
582 const std::string& timerName)
584 const std::string timerLabel =
585 prefix.empty() ? timerName : (prefix +
": " + timerName);
586 Teuchos::TimeMonitor::clearCounter (timerLabel);
588 #endif // BELOS_TEUCHOS_TIME_MONITOR 593 const std::vector<magnitude_type>& normsAfterSecondPass,
594 const std::vector<int>& faultIndices);
610 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
624 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
626 const bool attemptToRecycle =
true)
const;
639 const bool outOfPlace,
640 Teuchos::Array<mat_ptr> C,
642 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q);
651 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
652 Teuchos::ArrayView<mat_ptr> C)
const;
657 const Teuchos::RCP<const MV>& Q,
737 template<
class Scalar,
class MV>
742 using Teuchos::ParameterList;
743 using Teuchos::parameterList;
745 using Teuchos::sublist;
748 RCP<const ParameterList> defaultParams = getValidParameters ();
750 RCP<ParameterList> tsqrParams;
752 RCP<ParameterList> theParams;
753 if (params.is_null()) {
754 theParams = parameterList (*defaultParams);
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"));
783 if (! theParams->isSublist (
"TSQR implementation")) {
784 theParams->set (
"TSQR implementation",
785 defaultParams->sublist (
"TSQR implementation"));
787 tsqrParams = sublist (theParams,
"TSQR implementation",
true);
791 tsqrAdaptor_.setParameterList (tsqrParams);
794 setMyParamList (theParams);
797 template<
class Scalar,
class MV>
800 const std::string& label) :
804 randomizeNullSpace_ (true),
805 reorthogonalizeBlocks_ (true),
806 throwOnReorthogFault_ (false),
807 blockReorthogThreshold_ (0),
808 relativeRankTolerance_ (0),
809 forceNonnegativeDiagonal_ (false)
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 820 template<
class Scalar,
class MV>
826 randomizeNullSpace_ (true),
827 reorthogonalizeBlocks_ (true),
828 throwOnReorthogFault_ (false),
829 blockReorthogThreshold_ (0),
830 relativeRankTolerance_ (0),
831 forceNonnegativeDiagonal_ (false)
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 842 template<
class Scalar,
class MV>
845 norm (
const MV& X, std::vector<magnitude_type>& normVec)
const 847 const int numCols = MVT::GetNumberVecs (X);
850 if (normVec.size() <
static_cast<size_t>(numCols))
851 normVec.resize (numCols);
852 MVT::MvNorm (X, normVec);
855 template<
class Scalar,
class MV>
858 Teuchos::Array<mat_ptr> C,
859 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
861 #ifdef BELOS_TEUCHOS_TIME_MONITOR 869 Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
870 #endif // BELOS_TEUCHOS_TIME_MONITOR 872 int ncols_X, num_Q_blocks, ncols_Q_total;
873 checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X, Q);
877 if (ncols_X == 0 || num_Q_blocks == 0 || ncols_Q_total == 0)
881 allocateProjectionCoefficients (C, Q, X,
true);
885 std::vector<magnitude_type> columnNormsBefore (ncols_X,
magnitude_type(0));
886 if (reorthogonalizeBlocks_)
887 MVT::MvNorm (X, columnNormsBefore);
891 rawProject (X, Q, C);
895 if (reorthogonalizeBlocks_) {
896 std::vector<magnitude_type> columnNormsAfter (ncols_X,
magnitude_type(0));
897 MVT::MvNorm (X, columnNormsAfter);
905 bool reorthogonalize =
false;
906 for (
int j = 0; j < ncols_X; ++j) {
907 if (columnNormsAfter[j] < relThres * columnNormsBefore[j]) {
908 reorthogonalize =
true;
912 if (reorthogonalize) {
915 if (! reorthogCallback_.is_null()) {
916 reorthogCallback_->operator() (Teuchos::arrayViewFromVector (columnNormsBefore),
917 Teuchos::arrayViewFromVector (columnNormsAfter));
920 Teuchos::Array<mat_ptr> C2;
921 allocateProjectionCoefficients (C2, Q, X,
false);
925 rawProject (X, Q, C2);
927 for (
int k = 0; k < num_Q_blocks; ++k)
934 template<
class Scalar,
class MV>
938 using Teuchos::Range1D;
941 #ifdef BELOS_TEUCHOS_TIME_MONITOR 942 Teuchos::TimeMonitor timerMonitorNormalize(*timerNormalize_);
947 Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
948 #endif // BELOS_TEUCHOS_TIME_MONITOR 953 const int numCols = MVT::GetNumberVecs (X);
962 return normalizeOne (X, B);
992 MVT::GetGlobalLength(*Q_) != MVT::GetGlobalLength(X) ||
993 numCols > MVT::GetNumberVecs (*Q_)) {
994 Q_ = MVT::Clone (X, numCols);
1003 if (MVT::GetNumberVecs(*Q_) == numCols) {
1004 return normalizeImpl (X, *Q_, B,
false);
1006 RCP<MV> Q_view = MVT::CloneViewNonConst (*Q_, Range1D(0, numCols-1));
1007 return normalizeImpl (X, *Q_view, B,
false);
1011 template<
class Scalar,
class MV>
1015 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
1017 const bool attemptToRecycle)
const 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)
1024 for (
int i = 0; i < num_Q_blocks; ++i)
1026 const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
1030 C[i] = Teuchos::rcp (
new mat_type (ncols_Qi, ncols_X));
1034 if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X)
1035 Ci.shape (ncols_Qi, ncols_X);
1037 Ci.putScalar (SCT::zero());
1043 for (
int i = 0; i < num_Q_blocks; ++i)
1045 const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
1046 C[i] = Teuchos::rcp (
new mat_type (ncols_Qi, ncols_X));
1051 template<
class Scalar,
class MV>
1056 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1057 Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
1058 Teuchos::TimeMonitor timerMonitorNormalize(*timerNormalize_);
1059 #endif // BELOS_TEUCHOS_TIME_MONITOR 1061 const int numVecs = MVT::GetNumberVecs(X);
1064 }
else if (numVecs == 1) {
1066 using Teuchos::Range1D;
1071 const int rank = normalizeOne (X, B);
1073 RCP<MV> Q_0 = MVT::CloneViewNonConst (Q, Range1D(0,0));
1074 MVT::Assign (X, *Q_0);
1080 return normalizeImpl (X, Q, B,
true);
1084 template<
class Scalar,
class MV>
1089 const bool outOfPlace,
1090 Teuchos::Array<mat_ptr> C,
1092 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
1094 using Teuchos::Range1D;
1098 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1101 Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
1102 #endif // BELOS_TEUCHOS_TIME_MONITOR 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.");
1116 int ncols_X, num_Q_blocks, ncols_Q_total;
1117 checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X_in, Q);
1124 if (num_Q_blocks == 0 || ncols_Q_total == 0) {
1126 return normalizeOutOfPlace (X_in, X_out, B);
1128 return normalize (X_in, B);
1135 allocateProjectionCoefficients (C, Q, X_in,
true);
1140 std::vector<magnitude_type> normsBeforeFirstPass (ncols_X, SCTM::zero());
1141 if (reorthogonalizeBlocks_) {
1142 MVT::MvNorm (X_in, normsBeforeFirstPass);
1146 rawProject (X_in, Q, C);
1159 B_out = rcp (
new mat_type (ncols_X, ncols_X));
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() <<
".");
1171 B_out = rcp (
new mat_type (Teuchos::View, *B, ncols_X, ncols_X));
1177 const int firstPassRank = outOfPlace ?
1178 normalizeOutOfPlace (X_in, X_out, B_out) :
1179 normalize (X_in, B_out);
1188 int rank = firstPassRank;
1204 if (firstPassRank < ncols_X && randomizeNullSpace_) {
1205 const int numNullSpaceCols = ncols_X - firstPassRank;
1206 const Range1D nullSpaceIndices (firstPassRank, ncols_X - 1);
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));
1215 RCP<mat_type> B_null (
new mat_type (numNullSpaceCols, numNullSpaceCols));
1217 int randomVectorsRank;
1221 RCP<MV> X_out_null = MVT::CloneViewNonConst (X_out, nullSpaceIndices);
1226 RCP<MV> X_in_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
1227 MVT::Assign (*X_out_null, *X_in_null);
1230 rawProject (*X_in_null, Q, C_null);
1231 randomVectorsRank = normalizeOutOfPlace (*X_in_null, *X_out_null, B_null);
1235 RCP<MV> X_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
1238 rawProject (*X_null, Q, C_null);
1239 randomVectorsRank = normalize (*X_null, B_null);
1245 TEUCHOS_TEST_FOR_EXCEPTION(randomVectorsRank != numNullSpaceCols,
1247 "Belos::TsqrOrthoManagerImpl::projectAndNormalize" 1248 "OutOfPlace(): After projecting and normalizing the " 1249 "random vectors (used to replace the null space " 1250 "basis vectors from normalizing X), they have rank " 1251 << randomVectorsRank <<
", but should have full " 1252 "rank " << numNullSpaceCols <<
".");
1257 if (reorthogonalizeBlocks_) {
1260 std::vector<magnitude_type>
1261 normsAfterFirstPass (firstPassRank, SCTM::zero());
1262 std::vector<magnitude_type>
1263 normsAfterSecondPass (firstPassRank, SCTM::zero());
1277 Teuchos::BLAS<int, Scalar> blas;
1278 for (
int j = 0; j < firstPassRank; ++j) {
1279 const Scalar*
const B_j = &(*B_out)(0,j);
1282 normsAfterFirstPass[j] = blas.NRM2 (firstPassRank, B_j, 1);
1286 bool reorthogonalize =
false;
1287 for (
int j = 0; j < firstPassRank; ++j) {
1295 blockReorthogThreshold() * normsBeforeFirstPass[j];
1296 if (normsAfterFirstPass[j] < curThreshold) {
1297 reorthogonalize =
true;
1304 if (! reorthogCallback_.is_null()) {
1305 using Teuchos::arrayViewFromVector;
1306 (*reorthogCallback_) (arrayViewFromVector (normsBeforeFirstPass),
1307 arrayViewFromVector (normsAfterFirstPass));
1321 bool reorthogFault =
false;
1323 std::vector<int> faultIndices;
1324 if (reorthogonalize) {
1325 using Teuchos::Copy;
1326 using Teuchos::NO_TRANS;
1332 MVT::Assign (X_out, X_in);
1338 Teuchos::Array<mat_ptr> C2;
1339 allocateProjectionCoefficients (C2, Q, X_in,
false);
1344 rawProject (X_in, Q, C2);
1347 RCP<mat_type> B2 (
new mat_type (ncols_X, ncols_X));
1350 const int secondPassRank = outOfPlace ?
1351 normalizeOutOfPlace (X_in, X_out, B2) :
1352 normalize (X_in, B2);
1353 rank = secondPassRank;
1358 mat_type B_copy (Copy, *B_out, B_out->numRows(), B_out->numCols());
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");
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());
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");
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);
1388 reorthogFault =
false;
1389 for (
int j = 0; j < rank; ++j) {
1391 blockReorthogThreshold() * normsAfterFirstPass[j];
1392 if (normsAfterSecondPass[j] < relativeLowerBound) {
1393 reorthogFault =
true;
1394 faultIndices.push_back (j);
1399 if (reorthogFault) {
1400 if (throwOnReorthogFault_) {
1401 raiseReorthogFault (normsAfterFirstPass,
1402 normsAfterSecondPass,
1410 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
1411 "TsqrOrthoManagerImpl has not yet implemented" 1412 " recovery from an orthogonalization fault.");
1420 template<
class Scalar,
class MV>
1424 const std::vector<magnitude_type>& normsAfterSecondPass,
1425 const std::vector<int>& faultIndices)
1428 typedef std::vector<int>::size_type size_type;
1429 std::ostringstream os;
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];
1436 normsAfterSecondPass[index] / normsAfterFirstPass[index];
1437 os << index <<
"\t" << decreaseFactor << endl;
1442 template<
class Scalar,
class MV>
1443 Teuchos::RCP<const Teuchos::ParameterList>
1446 using Teuchos::ParameterList;
1447 using Teuchos::parameterList;
1450 if (defaultParams_.is_null()) {
1451 RCP<ParameterList> params = parameterList (
"TsqrOrthoManagerImpl");
1455 params->set (
"TSQR implementation", *(tsqrAdaptor_.getValidParameters()),
1456 "TSQR implementation parameters.");
1460 const bool defaultRandomizeNullSpace =
true;
1461 params->set (
"randomizeNullSpace", defaultRandomizeNullSpace,
1462 "Whether to fill in null space vectors with random data.");
1464 const bool defaultReorthogonalizeBlocks =
true;
1465 params->set (
"reorthogonalizeBlocks", defaultReorthogonalizeBlocks,
1466 "Whether to do block reorthogonalization as necessary.");
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.");
1480 Teuchos::as<magnitude_type>(10) * SCTM::eps();
1485 params->set (
"relativeRankTolerance", defaultRelativeRankTolerance,
1486 "Relative tolerance to determine the numerical rank of a " 1487 "block when normalizing.");
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).");
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.");
1502 defaultParams_ = params;
1504 return defaultParams_;
1507 template<
class Scalar,
class MV>
1508 Teuchos::RCP<const Teuchos::ParameterList>
1511 using Teuchos::ParameterList;
1515 RCP<const ParameterList> defaultParams = getValidParameters();
1517 RCP<ParameterList> params = rcp (
new ParameterList (*defaultParams));
1526 const bool randomizeNullSpace =
false;
1527 params->set (
"randomizeNullSpace", randomizeNullSpace);
1528 const bool reorthogonalizeBlocks =
false;
1529 params->set (
"reorthogonalizeBlocks", reorthogonalizeBlocks);
1534 template<
class Scalar,
class MV>
1539 Teuchos::SerialDenseMatrix<int, Scalar>& B)
1546 tsqrAdaptor_.factorExplicit (X, Q, B, forceNonnegativeDiagonal_);
1549 rank = tsqrAdaptor_.revealRank (Q, B, relativeRankTolerance_);
1550 }
catch (std::exception& e) {
1556 template<
class Scalar,
class MV>
1560 Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > B)
const 1567 B_out = Teuchos::rcp (
new mat_type (1, 1));
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 <<
".");
1576 B_out = Teuchos::rcp (
new mat_type (Teuchos::View, *B, 1, 1));
1580 std::vector<magnitude_type> theNorm (1, SCTM::zero());
1581 MVT::MvNorm (X, theNorm);
1582 (*B_out)(0,0) = theNorm[0];
1596 if (theNorm[0] == SCTM::zero()) {
1599 if (randomizeNullSpace_) {
1601 MVT::MvNorm (X, theNorm);
1602 if (theNorm[0] == SCTM::zero()) {
1608 "vector has norm zero!");
1613 const Scalar alpha = SCT::one() / theNorm[0];
1614 MVT::MvScale (X, alpha);
1621 const Scalar alpha = SCT::one() / theNorm[0];
1622 MVT::MvScale (X, alpha);
1628 template<
class Scalar,
class MV>
1632 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
1633 Teuchos::ArrayView<Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > > C)
const 1635 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1636 Teuchos::TimeMonitor timerMonitorNormalize(*timerProject_);
1637 #endif // BELOS_TEUCHOS_TIME_MONITOR 1640 const int num_Q_blocks = Q.size();
1641 for (
int i = 0; i < num_Q_blocks; ++i)
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);
1657 template<
class Scalar,
class MV>
1661 const Teuchos::RCP<const MV>& Q,
1662 const Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> >& C)
const 1664 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1665 Teuchos::TimeMonitor timerMonitorNormalize(*timerProject_);
1666 #endif // BELOS_TEUCHOS_TIME_MONITOR 1669 innerProd (*Q, X, *C);
1670 MVT::MvTimesMatAddMv (-SCT::one(), *Q, *C, SCT::one(), X);
1673 template<
class Scalar,
class MV>
1678 Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > B,
1679 const bool outOfPlace)
1681 using Teuchos::Range1D;
1684 using Teuchos::ScalarTraits;
1685 using Teuchos::tuple;
1687 const int numCols = MVT::GetNumberVecs (X);
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.");
1702 RCP<MV> Q_view = MVT::CloneViewNonConst (Q, Range1D (0, numCols-1));
1710 B_out = rcp (
new mat_type (numCols, 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() <<
".");
1720 B_out = rcp (
new mat_type (Teuchos::View, *B, numCols, numCols));
1732 const int rank = rawNormalize (X, *Q_view, *B_out);
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.");
1750 if (rank == numCols || ! randomizeNullSpace_) {
1754 MVT::Assign (*Q_view, X);
1759 if (randomizeNullSpace_ && rank < numCols) {
1766 const int nullSpaceNumCols = numCols - rank;
1769 Range1D nullSpaceIndices (rank, numCols-1);
1776 RCP<MV> Q_null = MVT::CloneViewNonConst (*Q_view, nullSpaceIndices);
1778 MVT::MvRandom (*Q_null);
1784 std::vector<magnitude_type> norms (MVT::GetNumberVecs (*Q_null));
1785 MVT::MvNorm (*Q_null, norms);
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()) {
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 " 1805 for (iter_type it = norms.begin(); it != norms.end(); ++it) {
1807 if (it+1 != norms.end())
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 " 1825 RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D (0, rank-1));
1831 rawProject (*Q_null, Q_col, C_null);
1840 RCP<MV> X_null = MVT::CloneViewNonConst (X, nullSpaceIndices);
1843 mat_type B_null (nullSpaceNumCols, nullSpaceNumCols);
1845 const int nullSpaceBasisRank = rawNormalize (*Q_null, *X_null, B_null);
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) {
1877 if (k != norms.size()-1) {
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.";
1886 TEUCHOS_TEST_FOR_EXCEPTION(nullSpaceBasisRank < nullSpaceNumCols,
1897 MVT::Assign (*X_null, *Q_null);
1898 }
else if (rank > 0) {
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);
1909 template<
class Scalar,
class MV>
1916 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 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);
1929 the_ncols_Q_total = 0;
1932 using Teuchos::ArrayView;
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);
1941 ncols_X = the_ncols_X;
1942 num_Q_blocks = the_num_Q_blocks;
1943 ncols_Q_total = the_ncols_Q_total;
1948 #endif // __BelosTsqrOrthoManagerImpl_hpp magnitude_type orthonormError(const MV &X) const
Return .
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Default valid parameter list.
MVT::tsqr_adaptor_type tsqr_adaptor_type
Teuchos::ScalarTraits< magnitude_type > SCTM
magnitude_type relativeRankTolerance_
Relative tolerance for measuring the numerical rank of a matrix.
int projectAndNormalize(MV &X, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Project X against Q and normalize X.
void norm(const MV &X, std::vector< magnitude_type > &normVec) const
Compute the 2-norm of each column j of X.
Interface of callback invoked by TsqrOrthoManager on reorthogonalization.
magnitude_type blockReorthogThreshold_
Relative reorthogonalization threshold in Block Gram-Schmidt.
void project(MV &X, Teuchos::Array< mat_ptr > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Compute and .
TsqrOrthoManagerImpl(const Teuchos::RCP< Teuchos::ParameterList > ¶ms, const std::string &label)
Constructor (that sets user-specified parameters).
void rawProject(MV &X, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q, Teuchos::ArrayView< mat_ptr > C) const
One projection pass of X against the Q[i] blocks.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
Type of the projection and normalization coefficients.
Teuchos::RCP< MV > Q_
Scratch space for TSQR.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Set parameters from the given parameter list.
int projectAndNormalizeImpl(MV &X_in, MV &X_out, const bool outOfPlace, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Implementation of projection and normalization.
magnitude_type eps_
Machine precision for Scalar.
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters()
Get "fast" parameters for TsqrOrthoManagerImpl.
TSQR-based OrthoManager subclass implementation.
Declaration of basic traits for the multivector type.
void setReorthogonalizationCallback(const Teuchos::RCP< ReorthogonalizationCallback< Scalar > > &callback)
Set callback to be invoked on reorthogonalization.
Teuchos::RCP< ReorthogonalizationCallback< Scalar > > reorthogCallback_
Callback invoked if reorthogonalization is necessary.
Error in TsqrOrthoManager or TsqrOrthoManagerImpl.
bool forceNonnegativeDiagonal_
Force R factor of normalization to have a nonnegative diagonal.
int normalize(MV &X, mat_ptr B)
Orthogonalize the columns of X in place.
bool throwOnReorthogFault_
Whether to throw an exception on a orthogonalization fault.
void setLabel(const std::string &label)
Set the label for timers.
Teuchos::RCP< mat_type > mat_ptr
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
void innerProd(const MV &X, const MV &Y, mat_type &Z) const
Euclidean inner product.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
The type of a norm result.
tsqr_adaptor_type tsqrAdaptor_
Interface to TSQR implementation.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< const Teuchos::ParameterList > defaultParams_
Default configuration parameters.
Scalar scalar_type
The template parameter of this class; the type of an inner product result.
void checkProjectionDims(int &ncols_X, int &num_Q_blocks, int &ncols_Q_total, const MV &X, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Return through output arguments some relevant dimension information about X and Q.
Teuchos::ScalarTraits< Scalar > SCT
virtual void operator()(Teuchos::ArrayView< magnitude_type > normsBeforeFirstPass, Teuchos::ArrayView< magnitude_type > normsAfterFirstPass)=0
Callback invoked by TsqrOrthoManager on reorthogonalization.
magnitude_type blockReorthogThreshold() const
Relative tolerance for triggering a block reorthogonalization.
Exception thrown to signal error in an orthogonalization manager method.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
TsqrOrthoError(const std::string &what_arg)
void allocateProjectionCoefficients(Teuchos::Array< mat_ptr > &C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q, const MV &X, const bool attemptToRecycle=true) const
Allocate projection coefficients.
int normalizeImpl(MV &X, MV &Q, mat_ptr B, const bool outOfPlace)
Normalize X into Q*B, with out-of-place option.
Teuchos::RCP< Teuchos::ParameterList > params_
Configuration parameters.
virtual ~ReorthogonalizationCallback()
Destructor (virtual for memory safety of derived classes)
Templated virtual class for providing orthogonalization/orthonormalization methods.
magnitude_type relativeRankTolerance() const
Relative tolerance for determining (via the SVD) whether a block is of full numerical rank...
int rawNormalize(MV &X, MV &Q, mat_type &B)
One out-of-place normalization pass.
bool reorthogonalizeBlocks_
Whether to reorthogonalize blocks at all.
TsqrOrthoFault(const std::string &what_arg)
void raiseReorthogFault(const std::vector< magnitude_type > &normsAfterFirstPass, const std::vector< magnitude_type > &normsAfterSecondPass, const std::vector< int > &faultIndices)
Throw an exception indicating a reorthgonalization fault.
std::string label_
Label for timers (if timers are used).
int projectAndNormalizeOutOfPlace(MV &X_in, MV &X_out, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Project and normalize X_in into X_out; overwrite X_in.
int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B)
Normalize X into Q*B, overwriting X.
const std::string & getLabel() const
Get the label for timers (if timers are enabled).
Belos header file which uses auto-configuration information to include necessary C++ headers...
bool randomizeNullSpace_
Whether to fill null space vectors with random data.
int normalizeOne(MV &X, mat_ptr B) const
Normalize a "multivector" of only one column.
magnitude_type orthogError(const MV &X1, const MV &X2) const
Return the Frobenius norm of the inner product of X1 with itself.
MultiVecTraits< Scalar, MV > MVT