32 #ifndef __AnasaziTsqrOrthoManagerImpl_hpp 33 #define __AnasaziTsqrOrthoManagerImpl_hpp 42 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp" 114 template<
class Scalar,
class MV>
118 typedef Scalar scalar_type;
120 typedef MV multivector_type;
130 typedef typename MVT::tsqr_adaptor_type tsqr_adaptor_type;
174 const std::string& label);
190 if (label != label_) {
196 const std::string&
getLabel ()
const {
return label_; }
230 norm (
const MV& X, std::vector<magnitude_type>& normvec)
const;
303 return projectAndNormalizeImpl (X, X,
false, C, B, Q);
334 return projectAndNormalizeImpl (X_in, X_out,
true, C, B, Q);
348 for (
int k = 0; k < ncols; ++k) {
361 mat_type X1_T_X2 (ncols_X1, ncols_X2);
386 tsqr_adaptor_type tsqrAdaptor_;
406 bool randomizeNullSpace_;
413 bool reorthogonalizeBlocks_;
418 bool throwOnReorthogFault_;
421 magnitude_type blockReorthogThreshold_;
424 magnitude_type relativeRankTolerance_;
432 bool forceNonnegativeDiagonal_;
436 raiseReorthogFault (
const std::vector<magnitude_type>& normsAfterFirstPass,
437 const std::vector<magnitude_type>& normsAfterSecondPass,
438 const std::vector<int>& faultIndices);
450 checkProjectionDims (
int& ncols_X,
470 const bool attemptToRecycle =
true)
const;
481 projectAndNormalizeImpl (MV& X_in,
483 const bool outOfPlace,
502 const mat_ptr& C)
const;
530 int rawNormalize (MV& X, MV& Q,
mat_type& B);
549 int normalizeOne (MV& X, mat_ptr B)
const;
578 int normalizeImpl (MV& X, MV& Q, mat_ptr B,
const bool outOfPlace);
581 template<
class Scalar,
class MV>
587 using Teuchos::parameterList;
589 using Teuchos::sublist;
590 typedef magnitude_type M;
592 RCP<const ParameterList> defaultParams = getValidParameters ();
594 RCP<ParameterList> tsqrParams;
596 RCP<ParameterList> theParams;
598 theParams = parameterList (*defaultParams);
606 randomizeNullSpace_ =
607 theParams->
get<
bool> (
"randomizeNullSpace",
608 defaultParams->get<
bool> (
"randomizeNullSpace"));
609 reorthogonalizeBlocks_ =
610 theParams->get<
bool> (
"reorthogonalizeBlocks",
611 defaultParams->get<
bool> (
"reorthogonalizeBlocks"));
612 throwOnReorthogFault_ =
613 theParams->get<
bool> (
"throwOnReorthogFault",
614 defaultParams->get<
bool> (
"throwOnReorthogFault"));
615 blockReorthogThreshold_ =
616 theParams->get<M> (
"blockReorthogThreshold",
617 defaultParams->get<M> (
"blockReorthogThreshold"));
618 relativeRankTolerance_ =
619 theParams->get<M> (
"relativeRankTolerance",
620 defaultParams->get<M> (
"relativeRankTolerance"));
621 forceNonnegativeDiagonal_ =
622 theParams->get<
bool> (
"forceNonnegativeDiagonal",
623 defaultParams->get<
bool> (
"forceNonnegativeDiagonal"));
627 if (! theParams->isSublist (
"TSQR implementation")) {
628 theParams->set (
"TSQR implementation",
629 defaultParams->sublist (
"TSQR implementation"));
631 tsqrParams = sublist (theParams,
"TSQR implementation",
true);
635 tsqrAdaptor_.setParameterList (tsqrParams);
638 setMyParamList (theParams);
641 template<
class Scalar,
class MV>
644 const std::string& label) :
648 randomizeNullSpace_ (true),
649 reorthogonalizeBlocks_ (true),
650 throwOnReorthogFault_ (false),
651 blockReorthogThreshold_ (0),
652 relativeRankTolerance_ (0),
653 forceNonnegativeDiagonal_ (false)
658 template<
class Scalar,
class MV>
664 randomizeNullSpace_ (true),
665 reorthogonalizeBlocks_ (true),
666 throwOnReorthogFault_ (false),
667 blockReorthogThreshold_ (0),
668 relativeRankTolerance_ (0),
669 forceNonnegativeDiagonal_ (false)
674 template<
class Scalar,
class MV>
677 norm (
const MV& X, std::vector<magnitude_type>& normVec)
const 679 const int numCols = MVT::GetNumberVecs (X);
682 if (normVec.size() <
static_cast<size_t>(numCols))
683 normVec.resize (numCols);
684 MVT::MvNorm (X, normVec);
687 template<
class Scalar,
class MV>
693 int ncols_X, num_Q_blocks, ncols_Q_total;
694 checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X, Q);
698 if (ncols_X == 0 || num_Q_blocks == 0 || ncols_Q_total == 0)
702 allocateProjectionCoefficients (C, Q, X,
true);
706 std::vector<magnitude_type> columnNormsBefore (ncols_X, magnitude_type(0));
707 if (reorthogonalizeBlocks_)
708 MVT::MvNorm (X, columnNormsBefore);
712 rawProject (X, Q, C);
716 if (reorthogonalizeBlocks_)
718 std::vector<magnitude_type> columnNormsAfter (ncols_X, magnitude_type(0));
719 MVT::MvNorm (X, columnNormsAfter);
722 const magnitude_type relThres = blockReorthogThreshold();
728 bool reorthogonalize =
false;
729 for (
int j = 0; j < ncols_X; ++j)
730 if (columnNormsAfter[j] < relThres * columnNormsBefore[j])
732 reorthogonalize =
true;
739 allocateProjectionCoefficients (C2, Q, X,
false);
743 rawProject (X, Q, C2);
745 for (
int k = 0; k < num_Q_blocks; ++k)
753 template<
class Scalar,
class MV>
763 const int numCols = MVT::GetNumberVecs (X);
772 return normalizeOne (X, B);
802 MVT::GetGlobalLength(*Q_) != MVT::GetGlobalLength(X) ||
803 numCols > MVT::GetNumberVecs (*Q_)) {
804 Q_ = MVT::Clone (X, numCols);
813 if (MVT::GetNumberVecs(*Q_) == numCols) {
814 return normalizeImpl (X, *Q_, B,
false);
816 RCP<MV> Q_view = MVT::CloneViewNonConst (*Q_, Range1D(0, numCols-1));
817 return normalizeImpl (X, *Q_view, B,
false);
821 template<
class Scalar,
class MV>
827 const bool attemptToRecycle)
const 829 const int num_Q_blocks = Q.size();
830 const int ncols_X = MVT::GetNumberVecs (X);
832 if (attemptToRecycle)
834 for (
int i = 0; i < num_Q_blocks; ++i)
836 const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
840 C[i] =
rcp (
new mat_type (ncols_Qi, ncols_X));
843 mat_type& Ci = *C[i];
844 if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X)
845 Ci.shape (ncols_Qi, ncols_X);
847 Ci.putScalar (SCT::zero());
853 for (
int i = 0; i < num_Q_blocks; ++i)
855 const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
856 C[i] =
rcp (
new mat_type (ncols_Qi, ncols_X));
861 template<
class Scalar,
class MV>
866 const int numVecs = MVT::GetNumberVecs(X);
869 }
else if (numVecs == 1) {
876 const int rank = normalizeOne (X, B);
878 RCP<MV> Q_0 = MVT::CloneViewNonConst (Q, Range1D(0,0));
879 MVT::Assign (X, *Q_0);
885 return normalizeImpl (X, Q, B,
true);
889 template<
class Scalar,
class MV>
894 const bool outOfPlace,
906 std::invalid_argument,
907 "Belos::TsqrOrthoManagerImpl::" 908 "projectAndNormalizeOutOfPlace(...):" 909 "X_out has " << MVT::GetNumberVecs(X_out)
910 <<
" columns, but X_in has " 911 << MVT::GetNumberVecs(X_in) <<
" columns.");
915 int ncols_X, num_Q_blocks, ncols_Q_total;
916 checkProjectionDims (ncols_X, num_Q_blocks, ncols_Q_total, X_in, Q);
923 if (num_Q_blocks == 0 || ncols_Q_total == 0) {
925 return normalizeOutOfPlace (X_in, X_out, B);
927 return normalize (X_in, B);
934 allocateProjectionCoefficients (C, Q, X_in,
true);
939 std::vector<magnitude_type> normsBeforeFirstPass (ncols_X, SCTM::zero());
940 if (reorthogonalizeBlocks_) {
941 MVT::MvNorm (X_in, normsBeforeFirstPass);
945 rawProject (X_in, Q, C);
958 B_out =
rcp (
new mat_type (ncols_X, ncols_X));
962 std::invalid_argument,
963 "normalizeOne: Input matrix B must be at " 964 "least " << ncols_X <<
" x " << ncols_X
965 <<
", but is instead " << B->numRows()
966 <<
" x " << B->numCols() <<
".");
976 const int firstPassRank = outOfPlace ?
977 normalizeOutOfPlace (X_in, X_out, B_out) :
978 normalize (X_in, B_out);
987 int rank = firstPassRank;
1003 if (firstPassRank < ncols_X && randomizeNullSpace_) {
1004 const int numNullSpaceCols = ncols_X - firstPassRank;
1005 const Range1D nullSpaceIndices (firstPassRank, ncols_X - 1);
1009 for (
int k = 0; k < num_Q_blocks; ++k) {
1010 const int numColsQk = MVT::GetNumberVecs(*Q[k]);
1011 C_null[k] =
rcp (
new mat_type (numColsQk, numNullSpaceCols));
1014 RCP<mat_type> B_null (
new mat_type (numNullSpaceCols, numNullSpaceCols));
1016 int randomVectorsRank;
1020 RCP<MV> X_out_null = MVT::CloneViewNonConst (X_out, nullSpaceIndices);
1025 RCP<MV> X_in_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
1026 MVT::Assign (*X_out_null, *X_in_null);
1029 rawProject (*X_in_null, Q, C_null);
1030 randomVectorsRank = normalizeOutOfPlace (*X_in_null, *X_out_null, B_null);
1034 RCP<MV> X_null = MVT::CloneViewNonConst (X_in, nullSpaceIndices);
1037 rawProject (*X_null, Q, C_null);
1038 randomVectorsRank = normalize (*X_null, B_null);
1046 "Belos::TsqrOrthoManagerImpl::projectAndNormalize" 1047 "OutOfPlace(): After projecting and normalizing the " 1048 "random vectors (used to replace the null space " 1049 "basis vectors from normalizing X), they have rank " 1050 << randomVectorsRank <<
", but should have full " 1051 "rank " << numNullSpaceCols <<
".");
1056 if (reorthogonalizeBlocks_) {
1059 std::vector<magnitude_type>
1060 normsAfterFirstPass (firstPassRank, SCTM::zero());
1061 std::vector<magnitude_type>
1062 normsAfterSecondPass (firstPassRank, SCTM::zero());
1077 for (
int j = 0; j < firstPassRank; ++j) {
1078 const Scalar*
const B_j = &(*B_out)(0,j);
1081 normsAfterFirstPass[j] = blas.
NRM2 (firstPassRank, B_j, 1);
1085 bool reorthogonalize =
false;
1086 for (
int j = 0; j < firstPassRank; ++j) {
1093 const magnitude_type curThreshold =
1094 blockReorthogThreshold() * normsBeforeFirstPass[j];
1095 if (normsAfterFirstPass[j] < curThreshold) {
1096 reorthogonalize =
true;
1111 bool reorthogFault =
false;
1113 std::vector<int> faultIndices;
1114 if (reorthogonalize) {
1122 MVT::Assign (X_out, X_in);
1128 allocateProjectionCoefficients (C2, Q, X_in,
false);
1133 rawProject (X_in, Q, C2);
1136 RCP<mat_type> B2 (
new mat_type (ncols_X, ncols_X));
1139 const int secondPassRank = outOfPlace ?
1140 normalizeOutOfPlace (X_in, X_out, B2) :
1141 normalize (X_in, B2);
1142 rank = secondPassRank;
1147 mat_type B_copy (Copy, *B_out, B_out->numRows(), B_out->numCols());
1149 const int err = B_out->multiply (NO_TRANS, NO_TRANS, SCT::one(),
1150 *B2, B_copy, SCT::zero());
1152 "Teuchos::SerialDenseMatrix::multiply " 1153 "returned err = " << err <<
" != 0");
1157 for (
int k = 0; k < num_Q_blocks; ++k) {
1158 mat_type C_k_copy (Copy, *C[k], C[k]->numRows(), C[k]->numCols());
1161 const int err2 = C[k]->multiply (NO_TRANS, NO_TRANS, SCT::one(),
1162 *C2[k], B_copy, SCT::one());
1164 "Teuchos::SerialDenseMatrix::multiply " 1165 "returned err = " << err <<
" != 0");
1170 for (
int j = 0; j < rank; ++j) {
1171 const Scalar*
const B2_j = &(*B2)(0,j);
1172 normsAfterSecondPass[j] = blas.
NRM2 (rank, B2_j, 1);
1177 reorthogFault =
false;
1178 for (
int j = 0; j < rank; ++j) {
1179 const magnitude_type relativeLowerBound =
1180 blockReorthogThreshold() * normsAfterFirstPass[j];
1181 if (normsAfterSecondPass[j] < relativeLowerBound) {
1182 reorthogFault =
true;
1183 faultIndices.push_back (j);
1188 if (reorthogFault) {
1189 if (throwOnReorthogFault_) {
1190 raiseReorthogFault (normsAfterFirstPass,
1191 normsAfterSecondPass,
1200 "TsqrOrthoManagerImpl has not yet implemented" 1201 " recovery from an orthogonalization fault.");
1209 template<
class Scalar,
class MV>
1211 TsqrOrthoManagerImpl<Scalar, MV>::
1212 raiseReorthogFault (
const std::vector<magnitude_type>& normsAfterFirstPass,
1213 const std::vector<magnitude_type>& normsAfterSecondPass,
1214 const std::vector<int>& faultIndices)
1217 typedef std::vector<int>::size_type size_type;
1218 std::ostringstream os;
1220 os <<
"Orthogonalization fault at the following column(s) of X:" << endl;
1221 os <<
"Column\tNorm decrease factor" << endl;
1222 for (size_type k = 0; k < faultIndices.size(); ++k)
1224 const int index = faultIndices[k];
1225 const magnitude_type decreaseFactor =
1226 normsAfterSecondPass[index] / normsAfterFirstPass[index];
1227 os << index <<
"\t" << decreaseFactor << endl;
1229 throw TsqrOrthoFault (os.str());
1232 template<
class Scalar,
class MV>
1237 using Teuchos::parameterList;
1240 if (defaultParams_.is_null()) {
1241 RCP<ParameterList> params = parameterList (
"TsqrOrthoManagerImpl");
1245 params->set (
"TSQR implementation", *(tsqrAdaptor_.getValidParameters()),
1246 "TSQR implementation parameters.");
1250 const bool defaultRandomizeNullSpace =
true;
1251 params->set (
"randomizeNullSpace", defaultRandomizeNullSpace,
1252 "Whether to fill in null space vectors with random data.");
1254 const bool defaultReorthogonalizeBlocks =
true;
1255 params->set (
"reorthogonalizeBlocks", defaultReorthogonalizeBlocks,
1256 "Whether to do block reorthogonalization as necessary.");
1260 const magnitude_type defaultBlockReorthogThreshold =
1261 magnitude_type(10) * SCTM::squareroot (SCTM::eps());
1262 params->set (
"blockReorthogThreshold", defaultBlockReorthogThreshold,
1263 "If reorthogonalizeBlocks==true, and if the norm of " 1264 "any column within a block decreases by this much or " 1265 "more after orthogonalization, we reorthogonalize.");
1269 const magnitude_type defaultRelativeRankTolerance =
1270 Teuchos::as<magnitude_type>(10) * SCTM::eps();
1275 params->set (
"relativeRankTolerance", defaultRelativeRankTolerance,
1276 "Relative tolerance to determine the numerical rank of a " 1277 "block when normalizing.");
1281 const bool defaultThrowOnReorthogFault =
true;
1282 params->set (
"throwOnReorthogFault", defaultThrowOnReorthogFault,
1283 "Whether to throw an exception if an orthogonalization " 1284 "fault occurs. This only matters if reorthogonalization " 1285 "is enabled (reorthogonalizeBlocks==true).");
1287 const bool defaultForceNonnegativeDiagonal =
false;
1288 params->set (
"forceNonnegativeDiagonal", defaultForceNonnegativeDiagonal,
1289 "Whether to force the R factor produced by the normalization " 1290 "step to have a nonnegative diagonal.");
1292 defaultParams_ = params;
1294 return defaultParams_;
1297 template<
class Scalar,
class MV>
1305 RCP<const ParameterList> defaultParams = getValidParameters();
1307 RCP<ParameterList> params =
rcp (
new ParameterList (*defaultParams));
1316 const bool randomizeNullSpace =
false;
1317 params->set (
"randomizeNullSpace", randomizeNullSpace);
1318 const bool reorthogonalizeBlocks =
false;
1319 params->set (
"reorthogonalizeBlocks", reorthogonalizeBlocks);
1324 template<
class Scalar,
class MV>
1336 tsqrAdaptor_.factorExplicit (X, Q, B, forceNonnegativeDiagonal_);
1339 rank = tsqrAdaptor_.revealRank (Q, B, relativeRankTolerance_);
1340 }
catch (std::exception& e) {
1341 throw TsqrOrthoError (e.what());
1346 template<
class Scalar,
class MV>
1348 TsqrOrthoManagerImpl<Scalar, MV>::
1349 normalizeOne (MV& X,
1357 B_out =
rcp (
new mat_type (1, 1));
1359 const int numRows = B->
numRows();
1360 const int numCols = B->
numCols();
1362 std::invalid_argument,
1363 "normalizeOne: Input matrix B must be at " 1364 "least 1 x 1, but is instead " << numRows
1365 <<
" x " << numCols <<
".");
1371 std::vector<magnitude_type> theNorm (1, SCTM::zero());
1372 MVT::MvNorm (X, theNorm);
1373 (*B_out)(0,0) = theNorm[0];
1387 if (theNorm[0] == SCTM::zero()) {
1390 if (randomizeNullSpace_) {
1392 MVT::MvNorm (X, theNorm);
1393 if (theNorm[0] == SCTM::zero()) {
1398 throw TsqrOrthoError(
"normalizeOne: a supposedly random " 1399 "vector has norm zero!");
1404 const Scalar alpha = SCT::one() / theNorm[0];
1405 MVT::MvScale (X, alpha);
1412 const Scalar alpha = SCT::one() / theNorm[0];
1413 MVT::MvScale (X, alpha);
1419 template<
class Scalar,
class MV>
1421 TsqrOrthoManagerImpl<Scalar, MV>::
1427 const int num_Q_blocks = Q.size();
1428 for (
int i = 0; i < num_Q_blocks; ++i)
1436 mat_type& Ci = *C[i];
1437 const MV& Qi = *Q[i];
1438 innerProd (Qi, X, Ci);
1439 MVT::MvTimesMatAddMv (-SCT::one(), Qi, Ci, SCT::one(), X);
1444 template<
class Scalar,
class MV>
1446 TsqrOrthoManagerImpl<Scalar, MV>::
1452 innerProd (*Q, X, *C);
1453 MVT::MvTimesMatAddMv (-SCT::one(), *Q, *C, SCT::one(), X);
1456 template<
class Scalar,
class MV>
1458 TsqrOrthoManagerImpl<Scalar, MV>::
1459 normalizeImpl (MV& X,
1462 const bool outOfPlace)
1468 using Teuchos::tuple;
1473 const bool extraDebug =
false;
1475 const int numCols = MVT::GetNumberVecs (X);
1483 std::invalid_argument,
1484 "TsqrOrthoManagerImpl::normalizeImpl(X,Q,B): Q has " 1485 << MVT::GetNumberVecs(Q) <<
" columns. This is too " 1486 "few, since X has " << numCols <<
" columns.");
1490 RCP<MV> Q_view = MVT::CloneViewNonConst (Q, Range1D(0, numCols-1));
1498 B_out =
rcp (
new mat_type (numCols, numCols));
1502 std::invalid_argument,
1503 "normalizeOne: Input matrix B must be at " 1504 "least " << numCols <<
" x " << numCols
1505 <<
", but is instead " << B->
numRows()
1506 <<
" x " << B->
numCols() <<
".");
1513 std::vector<magnitude_type> norms (numCols);
1514 MVT::MvNorm (X, norms);
1515 cerr <<
"Column norms of X before orthogonalization: ";
1516 typedef typename std::vector<magnitude_type>::const_iterator iter_type;
1517 for (iter_type iter = norms.begin(); iter != norms.end(); ++iter) {
1519 if (iter+1 != norms.end())
1533 const int rank = rawNormalize (X, *Q_view, *B_out);
1544 std::vector<magnitude_type> norms (numCols);
1545 MVT::MvNorm (*Q_view, norms);
1546 cerr <<
"Column norms of Q_view after orthogonalization: ";
1547 for (
typename std::vector<magnitude_type>::const_iterator iter = norms.begin();
1548 iter != norms.end(); ++iter) {
1550 if (iter+1 != norms.end())
1555 "Belos::TsqrOrthoManagerImpl::normalizeImpl: " 1556 "rawNormalize() returned rank = " << rank <<
" for a " 1557 "matrix X with " << numCols <<
" columns. Please report" 1558 " this bug to the Belos developers.");
1559 if (extraDebug && rank == 0) {
1562 const mat_type& B_ref = *B;
1563 std::vector<magnitude_type> norms (B_ref.numCols());
1564 for (
typename mat_type::ordinalType j = 0; j < B_ref.numCols(); ++j) {
1565 typedef typename mat_type::scalarType mat_scalar_type;
1566 mat_scalar_type sumOfSquares = ScalarTraits<mat_scalar_type>::zero();
1567 for (
typename mat_type::ordinalType i = 0; i <= j; ++i) {
1568 const mat_scalar_type B_ij =
1569 ScalarTraits<mat_scalar_type>::magnitude (B_ref(i,j));
1570 sumOfSquares += B_ij*B_ij;
1572 norms[j] = ScalarTraits<mat_scalar_type>::squareroot (sumOfSquares);
1576 cerr <<
"Norms of columns of B after orthogonalization: ";
1577 for (
typename mat_type::ordinalType j = 0; j < B_ref.numCols(); ++j) {
1579 if (j != B_ref.numCols() - 1)
1587 if (rank == numCols || ! randomizeNullSpace_) {
1591 MVT::Assign (*Q_view, X);
1596 if (randomizeNullSpace_ && rank < numCols) {
1603 const int nullSpaceNumCols = numCols - rank;
1606 Range1D nullSpaceIndices (rank, numCols-1);
1613 RCP<MV> Q_null = MVT::CloneViewNonConst (*Q_view, nullSpaceIndices);
1615 MVT::MvRandom (*Q_null);
1621 std::vector<magnitude_type> norms (MVT::GetNumberVecs(*Q_null));
1622 MVT::MvNorm (*Q_null, norms);
1624 bool anyZero =
false;
1625 typedef typename std::vector<magnitude_type>::const_iterator iter_type;
1626 for (iter_type it = norms.begin(); it != norms.end(); ++it) {
1627 if (*it == SCTM::zero()) {
1632 std::ostringstream os;
1633 os <<
"TsqrOrthoManagerImpl::normalizeImpl: " 1634 "We are being asked to randomize the null space, for a matrix " 1635 "with " << numCols <<
" columns and reported column rank " 1636 << rank <<
". The inclusive range of columns to fill with " 1637 "random data is [" << nullSpaceIndices.lbound() <<
"," 1638 << nullSpaceIndices.ubound() <<
"]. After filling the null " 1639 "space vectors with random numbers, at least one of the vectors" 1640 " has norm zero. Here are the norms of all the null space " 1642 for (iter_type it = norms.begin(); it != norms.end(); ++it) {
1644 if (it+1 != norms.end())
1647 os <<
"].) There is a tiny probability that this could happen " 1648 "randomly, but it is likely a bug. Please report it to the " 1649 "Belos developers, especially if you are able to reproduce the " 1662 RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D(0, rank-1));
1667 mat_ptr C_null (
new mat_type (rank, nullSpaceNumCols));
1668 rawProject (*Q_null, Q_col, C_null);
1677 RCP<MV> X_null = MVT::CloneViewNonConst (X, nullSpaceIndices);
1680 mat_type B_null (nullSpaceNumCols, nullSpaceNumCols);
1682 const int nullSpaceBasisRank = rawNormalize (*Q_null, *X_null, B_null);
1696 if (nullSpaceBasisRank < nullSpaceNumCols) {
1697 std::vector<magnitude_type> norms (MVT::GetNumberVecs(*X_null));
1698 MVT::MvNorm (*X_null, norms);
1699 std::ostringstream os;
1700 os <<
"TsqrOrthoManagerImpl::normalizeImpl: " 1701 <<
"We are being asked to randomize the null space, " 1702 <<
"for a matrix with " << numCols <<
" columns and " 1703 <<
"column rank " << rank <<
". After projecting and " 1704 <<
"normalizing the generated random vectors, they " 1705 <<
"only have rank " << nullSpaceBasisRank <<
". They" 1706 <<
" should have full rank " << nullSpaceNumCols
1707 <<
". (The inclusive range of columns to fill with " 1708 <<
"random data is [" << nullSpaceIndices.lbound()
1709 <<
"," << nullSpaceIndices.ubound() <<
"]. The " 1710 <<
"column norms of the resulting Q factor are: [";
1711 for (
typename std::vector<magnitude_type>::size_type k = 0;
1712 k < norms.size(); ++k) {
1714 if (k != norms.size()-1) {
1718 os <<
"].) There is a tiny probability that this could " 1719 <<
"happen randomly, but it is likely a bug. Please " 1720 <<
"report it to the Belos developers, especially if " 1721 <<
"you are able to reproduce the behavior.";
1724 TsqrOrthoError, os.str());
1734 MVT::Assign (*X_null, *Q_null);
1735 }
else if (rank > 0) {
1737 RCP<const MV> Q_col = MVT::CloneView (*Q_view, Range1D(0, rank-1));
1738 RCP<MV> X_col = MVT::CloneViewNonConst (X, Range1D(0, rank-1));
1739 MVT::Assign (*Q_col, *X_col);
1746 template<
class Scalar,
class MV>
1748 TsqrOrthoManagerImpl<Scalar, MV>::
1749 checkProjectionDims (
int& ncols_X,
1761 int the_ncols_X, the_num_Q_blocks, the_ncols_Q_total;
1762 the_num_Q_blocks = Q.size();
1763 the_ncols_X = MVT::GetNumberVecs (X);
1766 the_ncols_Q_total = 0;
1771 typedef typename ArrayView<RCP<const MV> >::const_iterator iter_type;
1772 for (iter_type it = Q.begin(); it != Q.end(); ++it) {
1773 const MV& Qi = **it;
1774 the_ncols_Q_total += MVT::GetNumberVecs (Qi);
1778 ncols_X = the_ncols_X;
1779 num_Q_blocks = the_num_Q_blocks;
1780 ncols_Q_total = the_ncols_Q_total;
1785 #endif // __AnasaziTsqrOrthoManagerImpl_hpp ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
void norm(const MV &X, std::vector< magnitude_type > &normvec) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Default valid parameter list.
magnitude_type relativeRankTolerance() const
bool is_null(const boost::shared_ptr< T > &p)
TSQR-based OrthoManager subclass implementation.
magnitude_type orthonormError(const MV &X) const
Return .
magnitude_type blockReorthogThreshold() const
magnitude_type orthogError(const MV &X1, const MV &X2) const
Return the Frobenius norm of the inner product of X1 with itself.
Declaration of basic traits for the multivector type.
T & get(const std::string &name, T def_value)
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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType n, const ScalarType *x, const OrdinalType incx) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Set parameters from the given parameter list.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
TsqrOrthoManagerImpl(const Teuchos::RCP< Teuchos::ParameterList > ¶ms, const std::string &label)
Constructor (that sets user-specified parameters).
TsqrOrthoManager(Impl) error.
const std::string & getLabel() const
Get the label for timers (if timers are enabled).
void setLabel(const std::string &label)
Set the label for timers.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
OrdinalType numRows() const
Templated virtual class for providing orthogonalization/orthonormalization methods.
Traits class which defines basic operations on multivectors.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters()
Get "fast" parameters for TsqrOrthoManagerImpl.
void project(MV &X, Teuchos::Array< mat_ptr > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Compute and .
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
void resize(size_type new_size, const value_type &x=value_type())
Exception thrown to signal error in an orthogonalization manager method.
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.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
int normalize(MV &X, mat_ptr B)
Orthogonalize the columns of X in place.
OrdinalType numCols() const