45 #ifndef __BelosTsqrOrthoManagerImpl_hpp
46 #define __BelosTsqrOrthoManagerImpl_hpp
54 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
55 #ifdef BELOS_TEUCHOS_TIME_MONITOR
57 #endif // BELOS_TEUCHOS_TIME_MONITOR
129 template<
class Scalar>
153 template<
class Scalar>
187 template<
class Scalar,
class MV>
202 typedef typename MVT::tsqr_adaptor_type tsqr_adaptor_type;
246 const std::string& label);
276 reorthogCallback_ = callback;
287 if (label != label_) {
290 #ifdef BELOS_TEUCHOS_TIME_MONITOR
291 clearTimer (label,
"All orthogonalization");
292 clearTimer (label,
"Projection");
293 clearTimer (label,
"Normalization");
295 timerOrtho_ = makeTimer (label,
"All orthogonalization");
296 timerProject_ = makeTimer (label,
"Projection");
297 timerNormalize_ = makeTimer (label,
"Normalization");
298 #endif // BELOS_TEUCHOS_TIME_MONITOR
303 const std::string&
getLabel ()
const {
return label_; }
337 norm (
const MV& X, std::vector<magnitude_type>& normVec)
const;
410 return projectAndNormalizeImpl (X, X,
false, C, B, Q);
441 return projectAndNormalizeImpl (X_in, X_out,
true, C, B, Q);
455 for (
int k = 0; k < ncols; ++k) {
468 mat_type X1_T_X2 (ncols_X1, ncols_X2);
493 tsqr_adaptor_type tsqrAdaptor_;
512 bool randomizeNullSpace_;
519 bool reorthogonalizeBlocks_;
524 bool throwOnReorthogFault_;
538 bool forceNonnegativeDiagonal_;
540 #ifdef BELOS_TEUCHOS_TIME_MONITOR
549 #endif // BELOS_TEUCHOS_TIME_MONITOR
554 #ifdef BELOS_TEUCHOS_TIME_MONITOR
564 makeTimer (
const std::string& prefix,
565 const std::string& timerName)
567 const std::string timerLabel =
568 prefix.empty() ? timerName : (prefix +
": " + timerName);
578 clearTimer (
const std::string& prefix,
579 const std::string& timerName)
581 const std::string timerLabel =
582 prefix.empty() ? timerName : (prefix +
": " + timerName);
585 #endif // BELOS_TEUCHOS_TIME_MONITOR
589 raiseReorthogFault (
const std::vector<magnitude_type>& normsAfterFirstPass,
590 const std::vector<magnitude_type>& normsAfterSecondPass,
591 const std::vector<int>& faultIndices);
603 checkProjectionDims (
int& ncols_X,
623 const bool attemptToRecycle =
true)
const;
634 projectAndNormalizeImpl (MV& X_in,
636 const bool outOfPlace,
683 int rawNormalize (MV& X, MV& Q,
mat_type& B);
702 int normalizeOne (MV& X,
mat_ptr B)
const;
731 int normalizeImpl (MV& X, MV& Q,
mat_ptr B,
const bool outOfPlace);
734 template<
class Scalar,
class MV>
740 using Teuchos::parameterList;
742 using Teuchos::sublist;
745 RCP<const ParameterList> defaultParams = getValidParameters ();
747 RCP<ParameterList> tsqrParams;
749 RCP<ParameterList> theParams;
751 theParams = parameterList (*defaultParams);
759 randomizeNullSpace_ =
760 theParams->
get<
bool> (
"randomizeNullSpace",
761 defaultParams->get<
bool> (
"randomizeNullSpace"));
762 reorthogonalizeBlocks_ =
763 theParams->get<
bool> (
"reorthogonalizeBlocks",
764 defaultParams->get<
bool> (
"reorthogonalizeBlocks"));
765 throwOnReorthogFault_ =
766 theParams->get<
bool> (
"throwOnReorthogFault",
767 defaultParams->get<
bool> (
"throwOnReorthogFault"));
768 blockReorthogThreshold_ =
769 theParams->get<M> (
"blockReorthogThreshold",
770 defaultParams->get<M> (
"blockReorthogThreshold"));
771 relativeRankTolerance_ =
772 theParams->get<M> (
"relativeRankTolerance",
773 defaultParams->get<M> (
"relativeRankTolerance"));
774 forceNonnegativeDiagonal_ =
775 theParams->get<
bool> (
"forceNonnegativeDiagonal",
776 defaultParams->get<
bool> (
"forceNonnegativeDiagonal"));
780 if (! theParams->isSublist (
"TSQR implementation")) {
781 theParams->set (
"TSQR implementation",
782 defaultParams->sublist (
"TSQR implementation"));
784 tsqrParams = sublist (theParams,
"TSQR implementation",
true);
788 tsqrAdaptor_.setParameterList (tsqrParams);
791 setMyParamList (theParams);
794 template<
class Scalar,
class MV>
797 const std::string& label) :
801 randomizeNullSpace_ (true),
802 reorthogonalizeBlocks_ (true),
803 throwOnReorthogFault_ (false),
804 blockReorthogThreshold_ (0),
805 relativeRankTolerance_ (0),
806 forceNonnegativeDiagonal_ (false)
810 #ifdef BELOS_TEUCHOS_TIME_MONITOR
811 timerOrtho_ = makeTimer (label,
"All orthogonalization");
812 timerProject_ = makeTimer (label,
"Projection");
813 timerNormalize_ = makeTimer (label,
"Normalization");
814 #endif // BELOS_TEUCHOS_TIME_MONITOR
817 template<
class Scalar,
class MV>
823 randomizeNullSpace_ (true),
824 reorthogonalizeBlocks_ (true),
825 throwOnReorthogFault_ (false),
826 blockReorthogThreshold_ (0),
827 relativeRankTolerance_ (0),
828 forceNonnegativeDiagonal_ (false)
832 #ifdef BELOS_TEUCHOS_TIME_MONITOR
833 timerOrtho_ = makeTimer (label,
"All orthogonalization");
834 timerProject_ = makeTimer (label,
"Projection");
835 timerNormalize_ = makeTimer (label,
"Normalization");
836 #endif // BELOS_TEUCHOS_TIME_MONITOR
839 template<
class Scalar,
class MV>
842 norm (
const MV& X, std::vector<magnitude_type>& normVec)
const
844 const int numCols = MVT::GetNumberVecs (X);
847 if (normVec.size() <
static_cast<size_t>(numCols))
848 normVec.resize (numCols);
849 MVT::MvNorm (X, normVec);
852 template<
class Scalar,
class MV>
858 #ifdef BELOS_TEUCHOS_TIME_MONITOR
867 #endif // BELOS_TEUCHOS_TIME_MONITOR
869 int ncols_X, num_Q_blocks, ncols_Q_total;
870 checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X, Q);
874 if (ncols_X == 0 || num_Q_blocks == 0 || ncols_Q_total == 0)
878 allocateProjectionCoefficients (C, Q, X,
true);
882 std::vector<magnitude_type> columnNormsBefore (ncols_X,
magnitude_type(0));
883 if (reorthogonalizeBlocks_)
884 MVT::MvNorm (X, columnNormsBefore);
888 rawProject (X, Q, C);
892 if (reorthogonalizeBlocks_) {
893 std::vector<magnitude_type> columnNormsAfter (ncols_X,
magnitude_type(0));
894 MVT::MvNorm (X, columnNormsAfter);
902 bool reorthogonalize =
false;
903 for (
int j = 0; j < ncols_X; ++j) {
904 if (columnNormsAfter[j] < relThres * columnNormsBefore[j]) {
905 reorthogonalize =
true;
909 if (reorthogonalize) {
912 if (! reorthogCallback_.is_null()) {
913 reorthogCallback_->operator() (Teuchos::arrayViewFromVector (columnNormsBefore),
914 Teuchos::arrayViewFromVector (columnNormsAfter));
918 allocateProjectionCoefficients (C2, Q, X,
false);
922 rawProject (X, Q, C2);
924 for (
int k = 0; k < num_Q_blocks; ++k)
931 template<
class Scalar,
class MV>
938 #ifdef BELOS_TEUCHOS_TIME_MONITOR
945 #endif // BELOS_TEUCHOS_TIME_MONITOR
950 const int numCols = MVT::GetNumberVecs (X);
959 return normalizeOne (X, B);
989 MVT::GetGlobalLength(*Q_) != MVT::GetGlobalLength(X) ||
990 numCols > MVT::GetNumberVecs (*Q_)) {
991 Q_ = MVT::Clone (X, numCols);
1000 if (MVT::GetNumberVecs(*Q_) == numCols) {
1001 return normalizeImpl (X, *Q_, B,
false);
1003 RCP<MV> Q_view = MVT::CloneViewNonConst (*Q_, Range1D(0, numCols-1));
1004 return normalizeImpl (X, *Q_view, B,
false);
1008 template<
class Scalar,
class MV>
1014 const bool attemptToRecycle)
const
1016 const int num_Q_blocks = Q.size();
1017 const int ncols_X = MVT::GetNumberVecs (X);
1019 if (attemptToRecycle)
1021 for (
int i = 0; i < num_Q_blocks; ++i)
1023 const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
1027 C[i] =
Teuchos::rcp (
new mat_type (ncols_Qi, ncols_X));
1030 mat_type& Ci = *C[i];
1031 if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X)
1032 Ci.shape (ncols_Qi, ncols_X);
1034 Ci.putScalar (SCT::zero());
1040 for (
int i = 0; i < num_Q_blocks; ++i)
1042 const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
1043 C[i] =
Teuchos::rcp (
new mat_type (ncols_Qi, ncols_X));
1048 template<
class Scalar,
class MV>
1053 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1056 #endif // BELOS_TEUCHOS_TIME_MONITOR
1058 const int numVecs = MVT::GetNumberVecs(X);
1061 }
else if (numVecs == 1) {
1068 const int rank = normalizeOne (X, B);
1070 RCP<MV> Q_0 = MVT::CloneViewNonConst (Q, Range1D(0,0));
1071 MVT::Assign (X, *Q_0);
1077 return normalizeImpl (X, Q, B,
true);
1081 template<
class Scalar,
class MV>
1086 const bool outOfPlace,
1095 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1099 #endif // BELOS_TEUCHOS_TIME_MONITOR
1104 std::invalid_argument,
1105 "Belos::TsqrOrthoManagerImpl::"
1106 "projectAndNormalizeImpl(..., outOfPlace=true, ...):"
1107 "X_out has " << MVT::GetNumberVecs(X_out)
1108 <<
" columns, but X_in has "
1109 << MVT::GetNumberVecs(X_in) <<
" columns.");
1113 int ncols_X, num_Q_blocks, ncols_Q_total;
1114 checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X_in, Q);
1121 if (num_Q_blocks == 0 || ncols_Q_total == 0) {
1123 return normalizeOutOfPlace (X_in, X_out, B);
1125 return normalize (X_in, B);
1132 allocateProjectionCoefficients (C, Q, X_in,
true);
1137 std::vector<magnitude_type> normsBeforeFirstPass (ncols_X, SCTM::zero());
1138 if (reorthogonalizeBlocks_) {
1139 MVT::MvNorm (X_in, normsBeforeFirstPass);
1143 rawProject (X_in, Q, C);
1156 B_out =
rcp (
new mat_type (ncols_X, ncols_X));
1160 std::invalid_argument,
1161 "normalizeOne: Input matrix B must be at "
1162 "least " << ncols_X <<
" x " << ncols_X
1163 <<
", but is instead " << B->numRows()
1164 <<
" x " << B->numCols() <<
".");
1174 const int firstPassRank = outOfPlace ?
1175 normalizeOutOfPlace (X_in, X_out, B_out) :
1176 normalize (X_in, B_out);
1185 int rank = firstPassRank;
1201 if (firstPassRank < ncols_X && randomizeNullSpace_) {
1202 const int numNullSpaceCols = ncols_X - firstPassRank;
1203 const Range1D nullSpaceIndices (firstPassRank, ncols_X - 1);
1207 for (
int k = 0; k < num_Q_blocks; ++k) {
1208 const int numColsQk = MVT::GetNumberVecs(*Q[k]);
1209 C_null[k] =
rcp (
new mat_type (numColsQk, numNullSpaceCols));
1212 RCP<mat_type> B_null (
new mat_type (numNullSpaceCols, numNullSpaceCols));
1214 int randomVectorsRank;
1218 RCP<MV> X_out_null = MVT::CloneViewNonConst (X_out, nullSpaceIndices);
1223 RCP<MV> X_in_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
1224 MVT::Assign (*X_out_null, *X_in_null);
1227 rawProject (*X_in_null, Q, C_null);
1228 randomVectorsRank = normalizeOutOfPlace (*X_in_null, *X_out_null, B_null);
1232 RCP<MV> X_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
1235 rawProject (*X_null, Q, C_null);
1236 randomVectorsRank = normalize (*X_null, B_null);
1244 "Belos::TsqrOrthoManagerImpl::projectAndNormalize"
1245 "OutOfPlace(): After projecting and normalizing the "
1246 "random vectors (used to replace the null space "
1247 "basis vectors from normalizing X), they have rank "
1248 << randomVectorsRank <<
", but should have full "
1249 "rank " << numNullSpaceCols <<
".");
1254 if (reorthogonalizeBlocks_) {
1257 std::vector<magnitude_type>
1258 normsAfterFirstPass (firstPassRank, SCTM::zero());
1259 std::vector<magnitude_type>
1260 normsAfterSecondPass (firstPassRank, SCTM::zero());
1275 for (
int j = 0; j < firstPassRank; ++j) {
1276 const Scalar*
const B_j = &(*B_out)(0,j);
1279 normsAfterFirstPass[j] = blas.
NRM2 (firstPassRank, B_j, 1);
1283 bool reorthogonalize =
false;
1284 for (
int j = 0; j < firstPassRank; ++j) {
1291 const magnitude_type curThreshold =
1292 blockReorthogThreshold() * normsBeforeFirstPass[j];
1293 if (normsAfterFirstPass[j] < curThreshold) {
1294 reorthogonalize =
true;
1301 if (! reorthogCallback_.is_null()) {
1302 using Teuchos::arrayViewFromVector;
1303 (*reorthogCallback_) (arrayViewFromVector (normsBeforeFirstPass),
1304 arrayViewFromVector (normsAfterFirstPass));
1318 bool reorthogFault =
false;
1320 std::vector<int> faultIndices;
1321 if (reorthogonalize) {
1329 MVT::Assign (X_out, X_in);
1336 allocateProjectionCoefficients (C2, Q, X_in,
false);
1341 rawProject (X_in, Q, C2);
1344 RCP<mat_type> B2 (
new mat_type (ncols_X, ncols_X));
1347 const int secondPassRank = outOfPlace ?
1348 normalizeOutOfPlace (X_in, X_out, B2) :
1349 normalize (X_in, B2);
1350 rank = secondPassRank;
1355 mat_type B_copy (Copy, *B_out, B_out->numRows(), B_out->numCols());
1357 const int err = B_out->multiply (NO_TRANS, NO_TRANS, SCT::one(),
1358 *B2, B_copy, SCT::zero());
1360 "Teuchos::SerialDenseMatrix::multiply "
1361 "returned err = " << err <<
" != 0");
1365 for (
int k = 0; k < num_Q_blocks; ++k) {
1366 mat_type C_k_copy (Copy, *C[k], C[k]->numRows(), C[k]->numCols());
1369 const int err1 = C[k]->multiply (NO_TRANS, NO_TRANS, SCT::one(),
1370 *C2[k], B_copy, SCT::one());
1372 "Teuchos::SerialDenseMatrix::multiply "
1373 "returned err = " << err1 <<
" != 0");
1378 for (
int j = 0; j < rank; ++j) {
1379 const Scalar*
const B2_j = &(*B2)(0,j);
1380 normsAfterSecondPass[j] = blas.
NRM2 (rank, B2_j, 1);
1385 reorthogFault =
false;
1386 for (
int j = 0; j < rank; ++j) {
1387 const magnitude_type relativeLowerBound =
1388 blockReorthogThreshold() * normsAfterFirstPass[j];
1389 if (normsAfterSecondPass[j] < relativeLowerBound) {
1390 reorthogFault =
true;
1391 faultIndices.push_back (j);
1396 if (reorthogFault) {
1397 if (throwOnReorthogFault_) {
1398 raiseReorthogFault (normsAfterFirstPass,
1399 normsAfterSecondPass,
1408 "TsqrOrthoManagerImpl has not yet implemented"
1409 " recovery from an orthogonalization fault.");
1417 template<
class Scalar,
class MV>
1419 TsqrOrthoManagerImpl<Scalar, MV>::
1420 raiseReorthogFault (
const std::vector<magnitude_type>& normsAfterFirstPass,
1421 const std::vector<magnitude_type>& normsAfterSecondPass,
1422 const std::vector<int>& faultIndices)
1425 typedef std::vector<int>::size_type size_type;
1426 std::ostringstream os;
1428 os <<
"Orthogonalization fault at the following column(s) of X:" << endl;
1429 os <<
"Column\tNorm decrease factor" << endl;
1430 for (size_type k = 0; k < faultIndices.size(); ++k) {
1431 const int index = faultIndices[k];
1432 const magnitude_type decreaseFactor =
1433 normsAfterSecondPass[index] / normsAfterFirstPass[index];
1434 os << index <<
"\t" << decreaseFactor << endl;
1436 throw TsqrOrthoFault (os.str());
1439 template<
class Scalar,
class MV>
1444 using Teuchos::parameterList;
1447 if (defaultParams_.is_null()) {
1448 RCP<ParameterList> params = parameterList (
"TsqrOrthoManagerImpl");
1452 params->set (
"TSQR implementation", *(tsqrAdaptor_.getValidParameters()),
1453 "TSQR implementation parameters.");
1457 const bool defaultRandomizeNullSpace =
true;
1458 params->set (
"randomizeNullSpace", defaultRandomizeNullSpace,
1459 "Whether to fill in null space vectors with random data.");
1461 const bool defaultReorthogonalizeBlocks =
true;
1462 params->set (
"reorthogonalizeBlocks", defaultReorthogonalizeBlocks,
1463 "Whether to do block reorthogonalization as necessary.");
1469 params->set (
"blockReorthogThreshold", defaultBlockReorthogThreshold,
1470 "If reorthogonalizeBlocks==true, and if the norm of "
1471 "any column within a block decreases by this much or "
1472 "more after orthogonalization, we reorthogonalize.");
1477 Teuchos::as<magnitude_type>(10) * SCTM::eps();
1482 params->set (
"relativeRankTolerance", defaultRelativeRankTolerance,
1483 "Relative tolerance to determine the numerical rank of a "
1484 "block when normalizing.");
1488 const bool defaultThrowOnReorthogFault =
true;
1489 params->set (
"throwOnReorthogFault", defaultThrowOnReorthogFault,
1490 "Whether to throw an exception if an orthogonalization "
1491 "fault occurs. This only matters if reorthogonalization "
1492 "is enabled (reorthogonalizeBlocks==true).");
1494 const bool defaultForceNonnegativeDiagonal =
false;
1495 params->set (
"forceNonnegativeDiagonal", defaultForceNonnegativeDiagonal,
1496 "Whether to force the R factor produced by the normalization "
1497 "step to have a nonnegative diagonal.");
1499 defaultParams_ = params;
1501 return defaultParams_;
1504 template<
class Scalar,
class MV>
1512 RCP<const ParameterList> defaultParams = getValidParameters();
1514 RCP<ParameterList> params =
rcp (
new ParameterList (*defaultParams));
1523 const bool randomizeNullSpace =
false;
1524 params->set (
"randomizeNullSpace", randomizeNullSpace);
1525 const bool reorthogonalizeBlocks =
false;
1526 params->set (
"reorthogonalizeBlocks", reorthogonalizeBlocks);
1531 template<
class Scalar,
class MV>
1543 tsqrAdaptor_.factorExplicit (X, Q, B, forceNonnegativeDiagonal_);
1546 rank = tsqrAdaptor_.revealRank (Q, B, relativeRankTolerance_);
1547 }
catch (std::exception& e) {
1548 throw TsqrOrthoError (e.what());
1553 template<
class Scalar,
class MV>
1555 TsqrOrthoManagerImpl<Scalar, MV>::
1556 normalizeOne (MV& X,
1566 const int theNumRows = B->
numRows ();
1567 const int theNumCols = B->
numCols ();
1569 theNumRows < 1 || theNumCols < 1, std::invalid_argument,
1570 "normalizeOne: Input matrix B must be at least 1 x 1, but "
1571 "is instead " << theNumRows <<
" x " << theNumCols <<
".");
1577 std::vector<magnitude_type> theNorm (1, SCTM::zero());
1578 MVT::MvNorm (X, theNorm);
1579 (*B_out)(0,0) = theNorm[0];
1593 if (theNorm[0] == SCTM::zero()) {
1596 if (randomizeNullSpace_) {
1598 MVT::MvNorm (X, theNorm);
1599 if (theNorm[0] == SCTM::zero()) {
1604 throw TsqrOrthoError(
"normalizeOne: a supposedly random "
1605 "vector has norm zero!");
1610 const Scalar alpha = SCT::one() / theNorm[0];
1611 MVT::MvScale (X, alpha);
1618 const Scalar alpha = SCT::one() / theNorm[0];
1619 MVT::MvScale (X, alpha);
1625 template<
class Scalar,
class MV>
1627 TsqrOrthoManagerImpl<Scalar, MV>::
1632 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1634 #endif // BELOS_TEUCHOS_TIME_MONITOR
1637 const int num_Q_blocks = Q.size();
1638 for (
int i = 0; i < num_Q_blocks; ++i)
1646 mat_type& Ci = *C[i];
1647 const MV& Qi = *Q[i];
1648 innerProd (Qi, X, Ci);
1649 MVT::MvTimesMatAddMv (-SCT::one(), Qi, Ci, SCT::one(), X);
1654 template<
class Scalar,
class MV>
1656 TsqrOrthoManagerImpl<Scalar, MV>::
1661 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1663 #endif // BELOS_TEUCHOS_TIME_MONITOR
1666 innerProd (*Q, X, *C);
1667 MVT::MvTimesMatAddMv (-SCT::one(), *Q, *C, SCT::one(), X);
1670 template<
class Scalar,
class MV>
1672 TsqrOrthoManagerImpl<Scalar, MV>::
1673 normalizeImpl (MV& X,
1676 const bool outOfPlace)
1682 using Teuchos::tuple;
1684 const int numCols = MVT::GetNumberVecs (X);
1692 MVT::GetNumberVecs (Q) < numCols, std::invalid_argument,
1693 "TsqrOrthoManagerImpl::normalizeImpl: Q has "
1694 << MVT::GetNumberVecs(Q) <<
" columns. This is too "
1695 "few, since X has " << numCols <<
" columns.");
1699 RCP<MV> Q_view = MVT::CloneViewNonConst (Q, Range1D (0, numCols-1));
1707 B_out =
rcp (
new mat_type (numCols, numCols));
1711 B->
numRows() < numCols || B->
numCols() < numCols, std::invalid_argument,
1712 "TsqrOrthoManagerImpl::normalizeImpl: Input matrix B must be at least "
1713 << numCols <<
" x " << numCols <<
", but is instead " << B->
numRows ()
1714 <<
" x " << B->
numCols() <<
".");
1729 const int rank = rawNormalize (X, *Q_view, *B_out);
1740 rank < 0 || rank > numCols, std::logic_error,
1741 "Belos::TsqrOrthoManagerImpl::normalizeImpl: rawNormalize returned rank "
1742 " = " << rank <<
" for a matrix X with " << numCols <<
" columns. "
1743 "Please report this bug to the Belos developers.");
1747 if (rank == numCols || ! randomizeNullSpace_) {
1751 MVT::Assign (*Q_view, X);
1756 if (randomizeNullSpace_ && rank < numCols) {
1763 const int nullSpaceNumCols = numCols - rank;
1766 Range1D nullSpaceIndices (rank, numCols-1);
1773 RCP<MV> Q_null = MVT::CloneViewNonConst (*Q_view, nullSpaceIndices);
1775 MVT::MvRandom (*Q_null);
1781 std::vector<magnitude_type> norms (MVT::GetNumberVecs (*Q_null));
1782 MVT::MvNorm (*Q_null, norms);
1784 bool anyZero =
false;
1785 typedef typename std::vector<magnitude_type>::const_iterator iter_type;
1786 for (iter_type it = norms.begin(); it != norms.end(); ++it) {
1787 if (*it == SCTM::zero()) {
1792 std::ostringstream os;
1793 os <<
"TsqrOrthoManagerImpl::normalizeImpl: "
1794 "We are being asked to randomize the null space, for a matrix "
1795 "with " << numCols <<
" columns and reported column rank "
1796 << rank <<
". The inclusive range of columns to fill with "
1797 "random data is [" << nullSpaceIndices.lbound() <<
","
1798 << nullSpaceIndices.ubound() <<
"]. After filling the null "
1799 "space vectors with random numbers, at least one of the vectors"
1800 " has norm zero. Here are the norms of all the null space "
1802 for (iter_type it = norms.begin(); it != norms.end(); ++it) {
1804 if (it+1 != norms.end())
1807 os <<
"].) There is a tiny probability that this could happen "
1808 "randomly, but it is likely a bug. Please report it to the "
1809 "Belos developers, especially if you are able to reproduce the "
1822 RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D (0, rank-1));
1827 mat_ptr C_null (
new mat_type (rank, nullSpaceNumCols));
1828 rawProject (*Q_null, Q_col, C_null);
1837 RCP<MV> X_null = MVT::CloneViewNonConst (X, nullSpaceIndices);
1840 mat_type B_null (nullSpaceNumCols, nullSpaceNumCols);
1842 const int nullSpaceBasisRank = rawNormalize (*Q_null, *X_null, B_null);
1856 if (nullSpaceBasisRank < nullSpaceNumCols) {
1857 std::vector<magnitude_type> norms (MVT::GetNumberVecs(*X_null));
1858 MVT::MvNorm (*X_null, norms);
1859 std::ostringstream os;
1860 os <<
"TsqrOrthoManagerImpl::normalizeImpl: "
1861 <<
"We are being asked to randomize the null space, "
1862 <<
"for a matrix with " << numCols <<
" columns and "
1863 <<
"column rank " << rank <<
". After projecting and "
1864 <<
"normalizing the generated random vectors, they "
1865 <<
"only have rank " << nullSpaceBasisRank <<
". They"
1866 <<
" should have full rank " << nullSpaceNumCols
1867 <<
". (The inclusive range of columns to fill with "
1868 <<
"random data is [" << nullSpaceIndices.lbound()
1869 <<
"," << nullSpaceIndices.ubound() <<
"]. The "
1870 <<
"column norms of the resulting Q factor are: [";
1871 for (
typename std::vector<magnitude_type>::size_type k = 0;
1872 k < norms.size(); ++k) {
1874 if (k != norms.size()-1) {
1878 os <<
"].) There is a tiny probability that this could "
1879 <<
"happen randomly, but it is likely a bug. Please "
1880 <<
"report it to the Belos developers, especially if "
1881 <<
"you are able to reproduce the behavior.";
1884 TsqrOrthoError, os.str ());
1894 MVT::Assign (*X_null, *Q_null);
1895 }
else if (rank > 0) {
1897 RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D (0, rank-1));
1898 RCP<MV> X_col = MVT::CloneViewNonConst (X, Range1D (0, rank-1));
1899 MVT::Assign (*Q_col, *X_col);
1906 template<
class Scalar,
class MV>
1908 TsqrOrthoManagerImpl<Scalar, MV>::
1909 checkProjectionDims (
int& ncols_X,
1921 int the_ncols_X, the_num_Q_blocks, the_ncols_Q_total;
1922 the_num_Q_blocks = Q.size();
1923 the_ncols_X = MVT::GetNumberVecs (X);
1926 the_ncols_Q_total = 0;
1931 typedef typename ArrayView<RCP<const MV> >::const_iterator iter_type;
1932 for (iter_type it = Q.begin(); it != Q.end(); ++it) {
1933 const MV& Qi = **it;
1934 the_ncols_Q_total += MVT::GetNumberVecs (Qi);
1938 ncols_X = the_ncols_X;
1939 num_Q_blocks = the_num_Q_blocks;
1940 ncols_Q_total = the_ncols_Q_total;
1945 #endif // __BelosTsqrOrthoManagerImpl_hpp
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.
Interface of callback invoked by TsqrOrthoManager on reorthogonalization.
bool is_null(const boost::shared_ptr< T > &p)
void project(MV &X, Teuchos::Array< mat_ptr > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Compute and .
void norm(const MV &X, std::vector< magnitude_type > &normVec) const
Compute the 2-norm of each column j of X.
TsqrOrthoManagerImpl(const Teuchos::RCP< Teuchos::ParameterList > ¶ms, const std::string &label)
Constructor (that sets user-specified parameters).
magnitude_type orthogError(const MV &X1, const MV &X2) const
Return the Frobenius norm of the inner product of X1 with itself.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
Type of the projection and normalization coefficients.
T & get(const std::string &name, T def_value)
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Set parameters from the given parameter list.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
Error in TsqrOrthoManager or TsqrOrthoManagerImpl.
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
int normalize(MV &X, mat_ptr B)
Orthogonalize the columns of X in place.
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 .
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
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.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Default valid parameter list.
Traits class which defines basic operations on multivectors.
Scalar scalar_type
The template parameter of this class; the type of an inner product result.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual void operator()(Teuchos::ArrayView< magnitude_type > normsBeforeFirstPass, Teuchos::ArrayView< magnitude_type > normsAfterFirstPass)=0
Callback invoked by TsqrOrthoManager on reorthogonalization.
const std::string & getLabel() const
Get the label for timers (if timers are enabled).
magnitude_type orthonormError(const MV &X) const
Return .
Exception thrown to signal error in an orthogonalization manager method.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
magnitude_type relativeRankTolerance() const
Relative tolerance for determining (via the SVD) whether a block is of full numerical rank...
void resize(size_type new_size, const value_type &x=value_type())
TsqrOrthoError(const std::string &what_arg)
OrdinalType numCols() const
virtual ~ReorthogonalizationCallback()
Destructor (virtual for memory safety of derived classes)
Templated virtual class for providing orthogonalization/orthonormalization methods.
TsqrOrthoFault(const std::string &what_arg)
magnitude_type blockReorthogThreshold() const
Relative tolerance for triggering a block 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.
int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B)
Normalize X into Q*B, overwriting X.
void innerProd(const MV &X, const MV &Y, mat_type &Z) const
Euclidean inner product.
Belos header file which uses auto-configuration information to include necessary C++ headers...
OrdinalType numRows() const