42 #ifndef __Belos_ProjectedLeastSquaresSolver_hpp 43 #define __Belos_ProjectedLeastSquaresSolver_hpp 47 #include "Teuchos_Array.hpp" 48 #include "Teuchos_BLAS.hpp" 49 #include "Teuchos_LAPACK.hpp" 50 #include "Teuchos_oblackholestream.hpp" 51 #include "Teuchos_ScalarTraits.hpp" 52 #include "Teuchos_SerialDenseMatrix.hpp" 53 #include "Teuchos_StandardParameterEntryValidators.hpp" 75 template<
class Scalar>
77 printMatrix (std::ostream& out,
78 const std::string& name,
79 const Teuchos::SerialDenseMatrix<int, Scalar>& A)
83 const int numRows = A.numRows();
84 const int numCols = A.numCols();
86 out << name <<
" = " << endl <<
'[';
89 for (
int i = 0; i < numRows; ++i) {
96 for (
int i = 0; i < numRows; ++i) {
97 for (
int j = 0; j < numCols; ++j) {
101 }
else if (i < numRows-1) {
112 template<
class Scalar>
114 print (std::ostream& out,
115 const Teuchos::SerialDenseMatrix<int, Scalar>& A,
116 const std::string& linePrefix)
120 const int numRows = A.numRows();
121 const int numCols = A.numCols();
123 out << linePrefix <<
'[';
126 for (
int i = 0; i < numRows; ++i) {
133 for (
int i = 0; i < numRows; ++i) {
134 for (
int j = 0; j < numCols; ++j) {
136 out << linePrefix <<
" ";
141 }
else if (i < numRows-1) {
147 out << linePrefix <<
']' << endl;
158 template<
class Scalar>
175 Teuchos::SerialDenseMatrix<int,Scalar>
H;
188 Teuchos::SerialDenseMatrix<int,Scalar>
R;
200 Teuchos::SerialDenseMatrix<int,Scalar>
y;
210 Teuchos::SerialDenseMatrix<int,Scalar>
z;
238 H (maxNumIterations+1, maxNumIterations),
239 R (maxNumIterations+1, maxNumIterations),
240 y (maxNumIterations+1, 1),
241 z (maxNumIterations+1, 1),
265 reset (
const typename Teuchos::ScalarTraits<Scalar>::magnitudeType beta)
267 typedef Teuchos::ScalarTraits<Scalar> STS;
270 z.putScalar (STS::zero());
274 const Scalar initialResidualNorm (beta);
275 z(0,0) = initialResidualNorm;
293 const int maxNumIterations)
295 typedef Teuchos::ScalarTraits<Scalar> STS;
296 typedef Teuchos::ScalarTraits<magnitude_type> STM;
298 TEUCHOS_TEST_FOR_EXCEPTION(beta < STM::zero(), std::invalid_argument,
299 "ProjectedLeastSquaresProblem::reset: initial " 300 "residual beta = " << beta <<
" < 0.");
301 TEUCHOS_TEST_FOR_EXCEPTION(maxNumIterations <= 0, std::invalid_argument,
302 "ProjectedLeastSquaresProblem::reset: maximum number " 303 "of iterations " << maxNumIterations <<
" <= 0.");
305 if (
H.numRows() < maxNumIterations+1 ||
H.numCols() < maxNumIterations) {
306 const int errcode =
H.reshape (maxNumIterations+1, maxNumIterations);
307 TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error,
308 "Failed to reshape H into a " << (maxNumIterations+1)
309 <<
" x " << maxNumIterations <<
" matrix.");
311 (void)
H.putScalar (STS::zero());
313 if (
R.numRows() < maxNumIterations+1 ||
R.numCols() < maxNumIterations) {
314 const int errcode =
R.reshape (maxNumIterations+1, maxNumIterations);
315 TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error,
316 "Failed to reshape R into a " << (maxNumIterations+1)
317 <<
" x " << maxNumIterations <<
" matrix.");
319 (void)
R.putScalar (STS::zero());
321 if (
y.numRows() < maxNumIterations+1 ||
y.numCols() < 1) {
322 const int errcode =
y.reshape (maxNumIterations+1, 1);
323 TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error,
324 "Failed to reshape y into a " << (maxNumIterations+1)
325 <<
" x " << 1 <<
" matrix.");
327 (void)
y.putScalar (STS::zero());
329 if (
z.numRows() < maxNumIterations+1 ||
z.numCols() < 1) {
330 const int errcode =
z.reshape (maxNumIterations+1, 1);
331 TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error,
332 "Failed to reshape z into a " << (maxNumIterations+1)
333 <<
" x " << 1 <<
" matrix.");
348 template<
class Scalar>
359 typedef Teuchos::SerialDenseMatrix<int,Scalar>
mat_type;
362 typedef Teuchos::ScalarTraits<scalar_type> STS;
363 typedef Teuchos::ScalarTraits<magnitude_type> STM;
364 typedef Teuchos::BLAS<int, scalar_type> blas_type;
365 typedef Teuchos::LAPACK<int, scalar_type> lapack_type;
372 for (
int i = 0; i < A.numRows(); ++i) {
373 for (
int j = 0; j < A.numCols(); ++j) {
374 A_star(j,i) = STS::conjugate (A(i,j));
383 const int N = R.numCols();
385 for (
int j = 0; j < N; ++j) {
386 for (
int i = 0; i <= j; ++i) {
387 L(j,i) = STS::conjugate (R(i,j));
396 const int N = std::min (A.numRows(), A.numCols());
398 for (
int j = 0; j < N; ++j) {
399 for (
int i = j+1; i < A.numRows(); ++i) {
400 A(i,j) = STS::zero();
415 Teuchos::RCP<mat_type>& A_21,
416 Teuchos::RCP<mat_type>& A_12,
417 Teuchos::RCP<mat_type>& A_22,
427 A_11 = rcp (
new mat_type (View, A, numRows1, numCols1, 0, 0));
428 A_21 = rcp (
new mat_type (View, A, numRows2, numCols1, numRows1, 0));
429 A_12 = rcp (
new mat_type (View, A, numRows1, numCols2, 0, numCols1));
430 A_22 = rcp (
new mat_type (View, A, numRows2, numCols2, numRows1, numCols1));
438 const int numRows = A.numRows();
439 const int numCols = A.numCols();
441 if (numRows == 0 || numCols == 0) {
444 for (
int j = 0; j < numCols; ++j) {
447 for (
int i = 0; i < numRows; ++i) {
463 const int numRows = Y.numRows();
464 const int numCols = Y.numCols();
466 TEUCHOS_TEST_FOR_EXCEPTION(numRows != X.numRows() || numCols != X.numCols(),
467 std::invalid_argument,
"Dimensions of X and Y don't " 468 "match. X is " << X.numRows() <<
" x " << X.numCols()
469 <<
", and Y is " << numRows <<
" x " << numCols <<
".");
470 for (
int j = 0; j < numCols; ++j) {
471 for (
int i = 0; i < numRows; ++i) {
472 Y(i,j) += alpha * X(i,j);
481 const int numRows = A.numRows();
482 const int numCols = A.numCols();
484 TEUCHOS_TEST_FOR_EXCEPTION(
485 B.numRows() != numRows || B.numCols() != numCols,
486 std::invalid_argument,
487 "matAdd: The input matrices A and B have incompatible dimensions. " 488 "A is " << numRows <<
" x " << numCols <<
", but B is " <<
489 B.numRows () <<
" x " << B.numCols () <<
".");
490 if (numRows == 0 || numCols == 0) {
493 for (
int j = 0; j < numCols; ++j) {
497 for (
int i = 0; i < numRows; ++i) {
508 const int numRows = A.numRows();
509 const int numCols = A.numCols();
511 TEUCHOS_TEST_FOR_EXCEPTION(
512 B.numRows() != numRows || B.numCols() != numCols,
513 std::invalid_argument,
514 "matSub: The input matrices A and B have incompatible dimensions. " 515 "A is " << numRows <<
" x " << numCols <<
", but B is " <<
516 B.numRows () <<
" x " << B.numCols () <<
".");
517 if (numRows == 0 || numCols == 0) {
520 for (
int j = 0; j < numCols; ++j) {
524 for (
int i = 0; i < numRows; ++i) {
539 TEUCHOS_TEST_FOR_EXCEPTION(B.numCols() != R.numRows(),
540 std::invalid_argument,
541 "rightUpperTriSolve: R and B have incompatible " 542 "dimensions. B has " << B.numCols() <<
" columns, " 543 "but R has " << R.numRows() <<
" rows.");
545 blas.TRSM (Teuchos::RIGHT_SIDE, Teuchos::UPPER_TRI,
546 Teuchos::NO_TRANS, Teuchos::NON_UNIT_DIAG,
547 R.numCols(), B.numCols(),
548 STS::one(), R.values(), R.stride(),
549 B.values(), B.stride());
564 using Teuchos::NO_TRANS;
566 TEUCHOS_TEST_FOR_EXCEPTION(A.numCols() != B.numRows(),
567 std::invalid_argument,
568 "matMatMult: The input matrices A and B have " 569 "incompatible dimensions. A is " << A.numRows()
570 <<
" x " << A.numCols() <<
", but B is " 571 << B.numRows() <<
" x " << B.numCols() <<
".");
572 TEUCHOS_TEST_FOR_EXCEPTION(A.numRows() != C.numRows(),
573 std::invalid_argument,
574 "matMatMult: The input matrix A and the output " 575 "matrix C have incompatible dimensions. A has " 576 << A.numRows() <<
" rows, but C has " << C.numRows()
578 TEUCHOS_TEST_FOR_EXCEPTION(B.numCols() != C.numCols(),
579 std::invalid_argument,
580 "matMatMult: The input matrix B and the output " 581 "matrix C have incompatible dimensions. B has " 582 << B.numCols() <<
" columns, but C has " 583 << C.numCols() <<
" columns.");
585 blas.GEMM (NO_TRANS, NO_TRANS, C.numRows(), C.numCols(), A.numCols(),
586 alpha, A.values(), A.stride(), B.values(), B.stride(),
587 beta, C.values(), C.stride());
601 for (
int j = 0; j < A.numCols(); ++j) {
602 if (upperTriangular) {
603 for (
int i = 0; i <= j && i < A.numRows(); ++i) {
604 if (STS::isnaninf (A(i,j))) {
609 for (
int i = 0; i < A.numRows(); ++i) {
610 if (STS::isnaninf (A(i,j))) {
624 std::pair<bool, std::pair<magnitude_type, magnitude_type> >
631 for (
int j = 0; j < A.numCols(); ++j) {
635 for (
int i = 0; i <= j && i < A.numRows(); ++i) {
637 upperTri += A_ij_mag * A_ij_mag;
640 for (
int i = j+1; i < A.numRows(); ++i) {
642 lowerTri += A_ij_mag * A_ij_mag;
643 if (A_ij_mag != STM::zero()) {
648 return std::make_pair (count == 0, std::make_pair (lowerTri, upperTri));
657 std::pair<bool, std::pair<magnitude_type, magnitude_type> >
664 for (
int j = 0; j < A.numCols(); ++j) {
668 for (
int i = 0; i <= j+1 && i < A.numRows(); ++i) {
670 upper += A_ij_mag * A_ij_mag;
673 for (
int i = j+2; i < A.numRows(); ++i) {
675 lower += A_ij_mag * A_ij_mag;
676 if (A_ij_mag != STM::zero()) {
681 return std::make_pair (count == 0, std::make_pair (lower, upper));
692 const char*
const matrixName)
const 694 std::pair<bool, std::pair<magnitude_type, magnitude_type> > result =
697 TEUCHOS_TEST_FOR_EXCEPTION(! result.first, std::invalid_argument,
698 "The " << A.numRows() <<
" x " << A.numCols()
699 <<
" matrix " << matrixName <<
" is not upper " 700 "triangular. ||tril(A)||_F = " 701 << result.second.first <<
" and ||A||_F = " 702 << result.second.second <<
".");
713 const char*
const matrixName)
const 715 std::pair<bool, std::pair<magnitude_type, magnitude_type> > result =
718 TEUCHOS_TEST_FOR_EXCEPTION(! result.first, std::invalid_argument,
719 "The " << A.numRows() <<
" x " << A.numCols()
720 <<
" matrix " << matrixName <<
" is not upper " 721 "triangular. ||tril(A(2:end, :))||_F = " 722 << result.second.first <<
" and ||A||_F = " 723 << result.second.second <<
".");
743 const char*
const matrixName,
746 std::pair<bool, std::pair<magnitude_type, magnitude_type> > result =
752 result.second.first :
753 result.second.first / result.second.second);
754 TEUCHOS_TEST_FOR_EXCEPTION(err > relativeTolerance, std::invalid_argument,
755 "The " << A.numRows() <<
" x " << A.numCols()
756 <<
" matrix " << matrixName <<
" is not upper " 757 "triangular. ||tril(A(2:end, :))||_F " 758 << (result.second.second == STM::zero() ?
"" :
" / ||A||_F")
759 <<
" = " << err <<
" > " << relativeTolerance <<
".");
776 const char*
const matrixName,
777 const int minNumRows,
778 const int minNumCols)
const 780 TEUCHOS_TEST_FOR_EXCEPTION(A.numRows() < minNumRows || A.numCols() < minNumCols,
781 std::invalid_argument,
782 "The matrix " << matrixName <<
" is " << A.numRows()
783 <<
" x " << A.numCols() <<
", and therefore does not " 784 "satisfy the minimum dimensions " << minNumRows
785 <<
" x " << minNumCols <<
".");
801 const char*
const matrixName,
803 const int numCols)
const 805 TEUCHOS_TEST_FOR_EXCEPTION(A.numRows() != numRows || A.numCols() != numCols,
806 std::invalid_argument,
807 "The matrix " << matrixName <<
" is supposed to be " 808 << numRows <<
" x " << numCols <<
", but is " 809 << A.numRows() <<
" x " << A.numCols() <<
" instead.");
843 const char* strings[] = {
"None",
"Some",
"Lots"};
845 std::invalid_argument,
846 "Invalid enum value " << x <<
".");
847 return std::string (strings[x]);
854 const char* strings[] = {
"None",
"Some",
"Lots"};
856 if (x == strings[r]) {
860 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
861 "Invalid robustness string " << x <<
".");
874 inline Teuchos::RCP<Teuchos::ParameterEntryValidator>
877 using Teuchos::stringToIntegralParameterEntryValidator;
879 Teuchos::Array<std::string> strs (3);
883 Teuchos::Array<std::string> docs (3);
884 docs[0] =
"Use the BLAS' triangular solve. This may result in Inf or " 885 "NaN output if the triangular matrix is rank deficient.";
886 docs[1] =
"Robustness somewhere between \"None\" and \"Lots\".";
887 docs[2] =
"Solve the triangular system in a least-squares sense, using " 888 "an SVD-based algorithm. This will always succeed, though the " 889 "solution may not make sense for GMRES.";
890 Teuchos::Array<ERobustness> ints (3);
894 const std::string pname (
"Robustness of Projected Least-Squares Solve");
896 return stringToIntegralParameterEntryValidator<ERobustness> (strs, docs,
962 template<
class Scalar>
979 typedef Teuchos::SerialDenseMatrix<int,Scalar>
mat_type;
982 typedef Teuchos::ScalarTraits<scalar_type> STS;
983 typedef Teuchos::ScalarTraits<magnitude_type> STM;
984 typedef Teuchos::BLAS<int, scalar_type> blas_type;
985 typedef Teuchos::LAPACK<int, scalar_type> lapack_type;
1004 defaultRobustness_ (defaultRobustness)
1024 return updateColumnGivens (problem.
H, problem.
R, problem.
y, problem.
z,
1047 return updateColumnsGivens (problem.
H, problem.
R, problem.
y, problem.
z,
1069 solveGivens (problem.
y, problem.
R, problem.
z, curCol);
1076 std::pair<int, bool>
1123 std::pair<int, bool>
1130 TEUCHOS_TEST_FOR_EXCEPTION(X.numRows() != B.numRows(), std::invalid_argument,
1131 "The output X and right-hand side B have different " 1132 "numbers of rows. X has " << X.numRows() <<
" rows" 1133 ", and B has " << B.numRows() <<
" rows.");
1138 TEUCHOS_TEST_FOR_EXCEPTION(X.numCols() > B.numCols(), std::invalid_argument,
1139 "The output X has more columns than the " 1140 "right-hand side B. X has " << X.numCols()
1141 <<
" columns and B has " << B.numCols()
1144 mat_type B_view (Teuchos::View, B, B.numRows(), X.numCols());
1159 std::pair<int, bool>
1174 std::pair<int, bool>
1180 using Teuchos::Array;
1181 using Teuchos::Copy;
1182 using Teuchos::LEFT_SIDE;
1183 using Teuchos::RIGHT_SIDE;
1186 const int M = R.numRows();
1187 const int N = R.numCols();
1188 TEUCHOS_TEST_FOR_EXCEPTION(M < N, std::invalid_argument,
1189 "The input matrix R has fewer columns than rows. " 1190 "R is " << M <<
" x " << N <<
".");
1192 mat_type R_view (Teuchos::View, R, N, N);
1194 if (side == LEFT_SIDE) {
1195 TEUCHOS_TEST_FOR_EXCEPTION(X.numRows() < N, std::invalid_argument,
1196 "The input/output matrix X has only " 1197 << X.numRows() <<
" rows, but needs at least " 1198 << N <<
" rows to match the matrix for a " 1199 "left-side solve R \\ X.");
1200 }
else if (side == RIGHT_SIDE) {
1201 TEUCHOS_TEST_FOR_EXCEPTION(X.numCols() < N, std::invalid_argument,
1202 "The input/output matrix X has only " 1203 << X.numCols() <<
" columns, but needs at least " 1204 << N <<
" columns to match the matrix for a " 1205 "right-side solve X / R.");
1209 std::invalid_argument,
1210 "Invalid robustness value " << robustness <<
".");
1217 TEUCHOS_TEST_FOR_EXCEPTION(count > 0, std::runtime_error,
1218 "There " << (count != 1 ?
"are" :
"is")
1219 <<
" " << count <<
" Inf or NaN entr" 1220 << (count != 1 ?
"ies" :
"y")
1221 <<
" in the upper triangle of R.");
1223 TEUCHOS_TEST_FOR_EXCEPTION(count > 0, std::runtime_error,
1224 "There " << (count != 1 ?
"are" :
"is")
1225 <<
" " << count <<
" Inf or NaN entr" 1226 << (count != 1 ?
"ies" :
"y") <<
" in the " 1227 "right-hand side(s) X.");
1232 bool foundRankDeficiency =
false;
1240 blas.TRSM(side, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
1241 Teuchos::NON_UNIT_DIAG, X.numRows(), X.numCols(),
1242 STS::one(), R.values(), R.stride(),
1243 X.values(), X.stride());
1247 mat_type B (Copy, X, X.numRows(), X.numCols());
1251 blas.TRSM(side, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
1252 Teuchos::NON_UNIT_DIAG, X.numRows(), X.numCols(),
1253 STS::one(), R.values(), R.stride(),
1254 X.values(), X.stride());
1260 warn_ <<
"Upper triangular solve: Found Infs and/or NaNs in the " 1261 "solution after using the fast algorithm. Retrying using a more " 1262 "robust algorithm." << std::endl;
1271 if (side == LEFT_SIDE) {
1273 mat_type R_copy (Teuchos::Copy, R_view, N, N);
1282 rank = solveLeastSquaresUsingSVD (R_copy, X);
1291 mat_type X_star (X.numCols(), X.numRows());
1297 rank = solveLeastSquaresUsingSVD (R_star, X_star);
1303 foundRankDeficiency =
true;
1307 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
1308 "Should never get here! Invalid robustness value " 1309 << robustness <<
". Please report this bug to the " 1310 "Belos developers.");
1312 return std::make_pair (rank, foundRankDeficiency);
1330 out <<
"Testing Givens rotations:" << endl;
1331 Scalar x = STS::random();
1332 Scalar y = STS::random();
1333 out <<
" x = " << x <<
", y = " << y << endl;
1335 Scalar theCosine, theSine, result;
1337 computeGivensRotation (x, y, theCosine, theSine, result);
1338 out <<
"-- After computing rotation:" << endl;
1339 out <<
"---- cos,sin = " << theCosine <<
"," << theSine << endl;
1340 out <<
"---- x = " << x <<
", y = " << y
1341 <<
", result = " << result << endl;
1343 blas.ROT (1, &x, 1, &y, 1, &theCosine, &theSine);
1344 out <<
"-- After applying rotation:" << endl;
1345 out <<
"---- cos,sin = " << theCosine <<
"," << theSine << endl;
1346 out <<
"---- x = " << x <<
", y = " << y << endl;
1349 if (STS::magnitude(y) > 2*STS::eps())
1387 const bool testBlockGivens=
false,
1388 const bool extraVerbose=
false)
1390 using Teuchos::Array;
1393 TEUCHOS_TEST_FOR_EXCEPTION(numCols <= 0, std::invalid_argument,
1394 "numCols = " << numCols <<
" <= 0.");
1395 const int numRows = numCols + 1;
1400 mat_type R_givens (numRows, numCols);
1403 Array<Scalar> theCosines (numCols);
1404 Array<Scalar> theSines (numCols);
1406 mat_type R_blockGivens (numRows, numCols);
1407 mat_type y_blockGivens (numRows, 1);
1408 mat_type z_blockGivens (numRows, 1);
1409 Array<Scalar> blockCosines (numCols);
1410 Array<Scalar> blockSines (numCols);
1411 const int panelWidth = std::min (3, numCols);
1413 mat_type R_lapack (numRows, numCols);
1418 makeRandomProblem (H, z);
1420 printMatrix<Scalar> (out,
"H", H);
1421 printMatrix<Scalar> (out,
"z", z);
1426 z_givens.assign (z);
1427 if (testBlockGivens) {
1428 z_blockGivens.assign (z);
1430 z_lapack.assign (z);
1438 for (
int curCol = 0; curCol < numCols; ++curCol) {
1439 residualNormGivens = updateColumnGivens (H, R_givens, y_givens, z_givens,
1440 theCosines, theSines, curCol);
1442 solveGivens (y_givens, R_givens, z_givens, numCols-1);
1447 if (testBlockGivens) {
1448 const bool testBlocksAtATime =
true;
1449 if (testBlocksAtATime) {
1451 for (
int startCol = 0; startCol < numCols; startCol += panelWidth) {
1452 int endCol = std::min (startCol + panelWidth - 1, numCols - 1);
1453 residualNormBlockGivens =
1454 updateColumnsGivens (H, R_blockGivens, y_blockGivens, z_blockGivens,
1455 blockCosines, blockSines, startCol, endCol);
1461 for (
int startCol = 0; startCol < numCols; ++startCol) {
1462 residualNormBlockGivens =
1463 updateColumnsGivens (H, R_blockGivens, y_blockGivens, z_blockGivens,
1464 blockCosines, blockSines, startCol, startCol);
1471 solveGivens (y_blockGivens, R_blockGivens, z_blockGivens, numCols-1);
1476 solveLapack (H, R_lapack, y_lapack, z_lapack, numCols-1);
1483 leastSquaresConditionNumber (H, z, residualNormLapack);
1494 mat_type y_givens_view (Teuchos::View, y_givens, numCols, 1);
1495 mat_type y_blockGivens_view (Teuchos::View, y_blockGivens, numCols, 1);
1496 mat_type y_lapack_view (Teuchos::View, y_lapack, numCols, 1);
1499 solutionError (y_givens_view, y_lapack_view);
1500 const magnitude_type blockGivensSolutionError = testBlockGivens ?
1501 solutionError (y_blockGivens_view, y_lapack_view) :
1508 mat_type R_factorFromGivens (numCols, numCols);
1509 mat_type R_factorFromBlockGivens (numCols, numCols);
1510 mat_type R_factorFromLapack (numCols, numCols);
1512 for (
int j = 0; j < numCols; ++j) {
1513 for (
int i = 0; i <= j; ++i) {
1514 R_factorFromGivens(i,j) = R_givens(i,j);
1515 if (testBlockGivens) {
1516 R_factorFromBlockGivens(i,j) = R_blockGivens(i,j);
1518 R_factorFromLapack(i,j) = R_lapack(i,j);
1522 printMatrix<Scalar> (out,
"R_givens", R_factorFromGivens);
1523 printMatrix<Scalar> (out,
"y_givens", y_givens_view);
1524 printMatrix<Scalar> (out,
"z_givens", z_givens);
1526 if (testBlockGivens) {
1527 printMatrix<Scalar> (out,
"R_blockGivens", R_factorFromBlockGivens);
1528 printMatrix<Scalar> (out,
"y_blockGivens", y_blockGivens_view);
1529 printMatrix<Scalar> (out,
"z_blockGivens", z_blockGivens);
1532 printMatrix<Scalar> (out,
"R_lapack", R_factorFromLapack);
1533 printMatrix<Scalar> (out,
"y_lapack", y_lapack_view);
1534 printMatrix<Scalar> (out,
"z_lapack", z_lapack);
1540 out <<
"||H||_F = " << H_norm << endl;
1542 out <<
"||H y_givens - z||_2 / ||H||_F = " 1543 << leastSquaresResidualNorm (H, y_givens_view, z) / H_norm << endl;
1544 if (testBlockGivens) {
1545 out <<
"||H y_blockGivens - z||_2 / ||H||_F = " 1546 << leastSquaresResidualNorm (H, y_blockGivens_view, z) / H_norm << endl;
1548 out <<
"||H y_lapack - z||_2 / ||H||_F = " 1549 << leastSquaresResidualNorm (H, y_lapack_view, z) / H_norm << endl;
1551 out <<
"||y_givens - y_lapack||_2 / ||y_lapack||_2 = " 1552 << givensSolutionError << endl;
1553 if (testBlockGivens) {
1554 out <<
"||y_blockGivens - y_lapack||_2 / ||y_lapack||_2 = " 1555 << blockGivensSolutionError << endl;
1558 out <<
"Least-squares condition number = " 1559 << leastSquaresCondNum << endl;
1575 10 * STM::squareroot( numRows*numCols );
1577 wiggleFactor * leastSquaresCondNum;
1579 solutionErrorBoundFactor * STS::eps();
1580 out <<
"Solution error bound: " << solutionErrorBoundFactor
1581 <<
" * eps = " << solutionErrorBound << endl;
1586 if (STM::isnaninf (solutionErrorBound)) {
1592 if (STM::isnaninf (givensSolutionError)) {
1594 }
else if (givensSolutionError > solutionErrorBound) {
1596 }
else if (testBlockGivens) {
1597 if (STM::isnaninf (blockGivensSolutionError)) {
1599 }
else if (blockGivensSolutionError > solutionErrorBound) {
1627 const int testProblemSize,
1629 const bool verbose=
false)
1631 using Teuchos::LEFT_SIDE;
1632 using Teuchos::RIGHT_SIDE;
1634 typedef Teuchos::SerialDenseMatrix<int, scalar_type>
mat_type;
1636 Teuchos::oblackholestream blackHole;
1637 std::ostream& verboseOut = verbose ? out : blackHole;
1639 verboseOut <<
"Testing upper triangular solves" << endl;
1643 verboseOut <<
"-- Generating test matrix" << endl;
1644 const int N = testProblemSize;
1647 for (
int j = 0; j < N; ++j) {
1648 for (
int i = 0; i <= j; ++i) {
1649 R(i,j) = STS::random ();
1661 mat_type R_copy (Teuchos::Copy, R, N, N);
1662 mat_type B_copy (Teuchos::Copy, B, N, 1);
1668 verboseOut <<
"-- Solving RX=B" << endl;
1673 Resid.assign (B_copy);
1675 ops.
matMatMult (STS::one(), Resid, -STS::one(), R_copy, X);
1676 verboseOut <<
"---- ||R*X - B||_F = " << Resid.normFrobenius() << endl;
1677 verboseOut <<
"---- ||R||_F ||X||_F + ||B||_F = " 1678 << (R_norm * X.normFrobenius() + B_norm)
1692 B_star_copy.assign (B_star);
1694 verboseOut <<
"-- Solving YR=B^*" << endl;
1699 Resid2.assign (B_star_copy);
1700 ops.
matMatMult (STS::one(), Resid2, -STS::one(), Y, R_copy);
1701 verboseOut <<
"---- ||Y*R - B^*||_F = " << Resid2.normFrobenius() << endl;
1702 verboseOut <<
"---- ||Y||_F ||R||_F + ||B^*||_F = " 1703 << (Y.normFrobenius() * R_norm + B_norm)
1716 std::ostream& warn_;
1739 using Teuchos::Array;
1744 TEUCHOS_TEST_FOR_EXCEPTION(count != 0, std::invalid_argument,
1745 "solveLeastSquaresUsingSVD: The input matrix A " 1746 "contains " << count <<
"Inf and/or NaN entr" 1747 << (count != 1 ?
"ies" :
"y") <<
".");
1749 TEUCHOS_TEST_FOR_EXCEPTION(count != 0, std::invalid_argument,
1750 "solveLeastSquaresUsingSVD: The input matrix X " 1751 "contains " << count <<
"Inf and/or NaN entr" 1752 << (count != 1 ?
"ies" :
"y") <<
".");
1754 const int N = std::min (A.numRows(), A.numCols());
1765 Array<magnitude_type> singularValues (N);
1773 Array<magnitude_type> rwork (1);
1774 if (STS::isComplex) {
1775 rwork.resize (std::max (1, 5 * N));
1780 Scalar lworkScalar = STS::one();
1782 lapack.GELSS (A.numRows(), A.numCols(), X.numCols(),
1783 A.values(), A.stride(), X.values(), X.stride(),
1784 &singularValues[0], rankTolerance, &rank,
1785 &lworkScalar, -1, &rwork[0], &info);
1786 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1787 "_GELSS workspace query returned INFO = " 1788 << info <<
" != 0.");
1789 const int lwork =
static_cast<int> (STS::real (lworkScalar));
1790 TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
1791 "_GELSS workspace query returned LWORK = " 1792 << lwork <<
" < 0.");
1794 Array<Scalar> work (std::max (1, lwork));
1796 lapack.GELSS (A.numRows(), A.numCols(), X.numCols(),
1797 A.values(), A.stride(), X.values(), X.stride(),
1798 &singularValues[0], rankTolerance, &rank,
1799 &work[0], lwork, &rwork[0], &info);
1800 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
1801 "_GELSS returned INFO = " << info <<
" != 0.");
1817 const int numRows = curCol + 2;
1822 const mat_type R_view (Teuchos::View, R, numRows-1, numRows-1);
1823 const mat_type z_view (Teuchos::View, z, numRows-1, z.numCols());
1824 mat_type y_view (Teuchos::View, y, numRows-1, y.numCols());
1827 R_view, z_view, defaultRobustness_);
1841 for (
int j = 0; j < H.numCols(); ++j) {
1842 for (
int i = j+2; i < H.numRows(); ++i) {
1843 H(i,j) = STS::zero();
1854 const int numTrials = 1000;
1856 for (
int trial = 0; trial < numTrials && z_init == STM::zero(); ++trial) {
1857 z_init = STM::random();
1859 TEUCHOS_TEST_FOR_EXCEPTION(z_init == STM::zero(), std::runtime_error,
1860 "After " << numTrials <<
" trial" 1861 << (numTrials != 1 ?
"s" :
"")
1862 <<
", we were unable to generate a nonzero pseudo" 1863 "random real number. This most likely indicates a " 1864 "broken pseudorandom number generator.");
1877 computeGivensRotation (
const Scalar& x,
1885 const bool useLartg =
false;
1890 lapack.LARTG (x, y, &theCosine, &theSine, &result);
1899 blas.ROTG (&x_temp, &y_temp, &theCosine, &theSine);
1907 Teuchos::ArrayView<magnitude_type> sigmas)
1909 using Teuchos::Array;
1910 using Teuchos::ArrayView;
1912 const int numRows = A.numRows();
1913 const int numCols = A.numCols();
1914 TEUCHOS_TEST_FOR_EXCEPTION(sigmas.size() < std::min (numRows, numCols),
1915 std::invalid_argument,
1916 "The sigmas array is only of length " << sigmas.size()
1917 <<
", but must be of length at least " 1918 << std::min (numRows, numCols)
1919 <<
" in order to hold all the singular values of the " 1925 mat_type A_copy (numRows, numCols);
1931 Scalar lworkScalar = STS::zero();
1932 Array<magnitude_type> rwork (std::max (std::min (numRows, numCols) - 1, 1));
1933 lapack.GESVD (
'N',
'N', numRows, numCols,
1934 A_copy.values(), A_copy.stride(), &sigmas[0],
1935 (Scalar*) NULL, 1, (Scalar*) NULL, 1,
1936 &lworkScalar, -1, &rwork[0], &info);
1938 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1939 "LAPACK _GESVD workspace query failed with INFO = " 1941 const int lwork =
static_cast<int> (STS::real (lworkScalar));
1942 TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
1943 "LAPACK _GESVD workspace query returned LWORK = " 1944 << lwork <<
" < 0.");
1947 Teuchos::Array<Scalar> work (std::max (1, lwork));
1950 lapack.GESVD (
'N',
'N', numRows, numCols,
1951 A_copy.values(), A_copy.stride(), &sigmas[0],
1952 (Scalar*) NULL, 1, (Scalar*) NULL, 1,
1953 &work[0], lwork, &rwork[0], &info);
1954 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1955 "LAPACK _GESVD failed with INFO = " << info <<
".");
1965 std::pair<magnitude_type, magnitude_type>
1966 extremeSingularValues (
const mat_type& A)
1968 using Teuchos::Array;
1970 const int numRows = A.numRows();
1971 const int numCols = A.numCols();
1973 Array<magnitude_type> sigmas (std::min (numRows, numCols));
1974 singularValues (A, sigmas);
1975 return std::make_pair (sigmas[0], sigmas[std::min(numRows, numCols) - 1]);
1989 leastSquaresConditionNumber (
const mat_type& A,
1994 const std::pair<magnitude_type, magnitude_type> sigmaMaxMin =
1995 extremeSingularValues (A);
1999 TEUCHOS_TEST_FOR_EXCEPTION(sigmaMaxMin.second == STM::zero(), std::runtime_error,
2000 "The test matrix is rank deficient; LAPACK's _GESVD " 2001 "routine reports that its smallest singular value is " 2005 const magnitude_type A_cond = sigmaMaxMin.first / sigmaMaxMin.second;
2011 const magnitude_type sinTheta = residualNorm / b.normFrobenius();
2024 STM::zero() : STM::squareroot (1 - sinTheta * sinTheta);
2032 return 2 * A_cond / cosTheta + tanTheta * A_cond * A_cond;
2037 leastSquaresResidualNorm (
const mat_type& A,
2041 mat_type r (b.numRows(), b.numCols());
2045 LocalDenseMatrixOps<Scalar> ops;
2046 ops.
matMatMult (STS::one(), r, -STS::one(), A, x);
2047 return r.normFrobenius ();
2055 solutionError (
const mat_type& x_approx,
2058 const int numRows = x_exact.numRows();
2059 const int numCols = x_exact.numCols();
2061 mat_type x_diff (numRows, numCols);
2062 for (
int j = 0; j < numCols; ++j) {
2063 for (
int i = 0; i < numRows; ++i) {
2064 x_diff(i,j) = x_exact(i,j) - x_approx(i,j);
2070 return x_diff.normFrobenius() /
2071 (scalingFactor == STM::zero() ? STM::one() : scalingFactor);
2121 updateColumnGivens (
const mat_type& H,
2125 Teuchos::ArrayView<scalar_type> theCosines,
2126 Teuchos::ArrayView<scalar_type> theSines,
2132 const int numRows = curCol + 2;
2133 const int LDR = R.stride();
2134 const bool extraDebug =
false;
2137 cerr <<
"updateColumnGivens: curCol = " << curCol << endl;
2143 const mat_type H_col (Teuchos::View, H, numRows, 1, 0, curCol);
2147 mat_type R_col (Teuchos::View, R, numRows, 1, 0, curCol);
2151 R_col.assign (H_col);
2154 printMatrix<Scalar> (cerr,
"R_col before ", R_col);
2160 for (
int j = 0; j < curCol; ++j) {
2163 Scalar theCosine = theCosines[j];
2164 Scalar theSine = theSines[j];
2167 cerr <<
" j = " << j <<
": (cos,sin) = " 2168 << theCosines[j] <<
"," << theSines[j] << endl;
2170 blas.ROT (1, &R_col(j,0), LDR, &R_col(j+1,0), LDR,
2171 &theCosine, &theSine);
2173 if (extraDebug && curCol > 0) {
2174 printMatrix<Scalar> (cerr,
"R_col after applying previous " 2175 "Givens rotations", R_col);
2180 Scalar theCosine, theSine, result;
2181 computeGivensRotation (R_col(curCol,0), R_col(curCol+1,0),
2182 theCosine, theSine, result);
2183 theCosines[curCol] = theCosine;
2184 theSines[curCol] = theSine;
2187 cerr <<
" New cos,sin = " << theCosine <<
"," << theSine << endl;
2193 R_col(curCol, 0) = result;
2194 R_col(curCol+1, 0) = STS::zero();
2197 printMatrix<Scalar> (cerr,
"R_col after applying current " 2198 "Givens rotation", R_col);
2206 const int LDZ = z.stride();
2207 blas.ROT (z.numCols(),
2208 &z(curCol,0), LDZ, &z(curCol+1,0), LDZ,
2209 &theCosine, &theSine);
2215 printMatrix<Scalar> (cerr,
"z_after", z);
2221 return STS::magnitude( z(numRows-1, 0) );
2264 const int numRows = curCol + 2;
2265 const int numCols = curCol + 1;
2266 const int LDR = R.stride();
2269 const mat_type H_view (Teuchos::View, H, numRows, numCols);
2270 mat_type R_view (Teuchos::View, R, numRows, numCols);
2271 R_view.assign (H_view);
2275 mat_type y_view (Teuchos::View, y, numRows, y.numCols());
2276 mat_type z_view (Teuchos::View, z, numRows, y.numCols());
2277 y_view.assign (z_view);
2281 Scalar lworkScalar = STS::zero();
2283 lapack.GELS (
'N', numRows, numCols, y_view.numCols(),
2284 NULL, LDR, NULL, y_view.stride(),
2285 &lworkScalar, -1, &info);
2286 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
2287 "LAPACK _GELS workspace query failed with INFO = " 2288 << info <<
", for a " << numRows <<
" x " << numCols
2289 <<
" matrix with " << y_view.numCols()
2290 <<
" right hand side" 2291 << ((y_view.numCols() != 1) ?
"s" :
"") <<
".");
2292 TEUCHOS_TEST_FOR_EXCEPTION(STS::real(lworkScalar) < STM::zero(),
2294 "LAPACK _GELS workspace query returned an LWORK with " 2295 "negative real part: LWORK = " << lworkScalar
2296 <<
". That should never happen. Please report this " 2297 "to the Belos developers.");
2298 TEUCHOS_TEST_FOR_EXCEPTION(STS::isComplex && STS::imag(lworkScalar) != STM::zero(),
2300 "LAPACK _GELS workspace query returned an LWORK with " 2301 "nonzero imaginary part: LWORK = " << lworkScalar
2302 <<
". That should never happen. Please report this " 2303 "to the Belos developers.");
2309 const int lwork = std::max (1, static_cast<int> (STS::real (lworkScalar)));
2312 Teuchos::Array<Scalar> work (lwork);
2317 lapack.GELS (
'N', numRows, numCols, y_view.numCols(),
2318 R_view.values(), R_view.stride(),
2319 y_view.values(), y_view.stride(),
2320 (lwork > 0 ? &work[0] : (Scalar*) NULL),
2323 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
2324 "Solving projected least-squares problem with LAPACK " 2325 "_GELS failed with INFO = " << info <<
", for a " 2326 << numRows <<
" x " << numCols <<
" matrix with " 2327 << y_view.numCols() <<
" right hand side" 2328 << (y_view.numCols() != 1 ?
"s" :
"") <<
".");
2332 return STS::magnitude( y_view(numRows-1, 0) );
2346 updateColumnsGivens (
const mat_type& H,
2350 Teuchos::ArrayView<scalar_type> theCosines,
2351 Teuchos::ArrayView<scalar_type> theSines,
2355 TEUCHOS_TEST_FOR_EXCEPTION(startCol > endCol, std::invalid_argument,
2356 "updateColumnGivens: startCol = " << startCol
2357 <<
" > endCol = " << endCol <<
".");
2360 for (
int curCol = startCol; curCol <= endCol; ++curCol) {
2361 lastResult = updateColumnGivens (H, R, y, z, theCosines, theSines, curCol);
2379 updateColumnsGivensBlock (
const mat_type& H,
2383 Teuchos::ArrayView<scalar_type> theCosines,
2384 Teuchos::ArrayView<scalar_type> theSines,
2388 const int numRows = endCol + 2;
2389 const int numColsToUpdate = endCol - startCol + 1;
2390 const int LDR = R.stride();
2395 const mat_type H_view (Teuchos::View, H, numRows, numColsToUpdate, 0, startCol);
2396 mat_type R_view (Teuchos::View, R, numRows, numColsToUpdate, 0, startCol);
2397 R_view.assign (H_view);
2405 for (
int j = 0; j < startCol; ++j) {
2406 blas.ROT (numColsToUpdate,
2407 &R(j, startCol), LDR, &R(j+1, startCol), LDR,
2408 &theCosines[j], &theSines[j]);
2412 for (
int curCol = startCol; curCol < endCol; ++curCol) {
2415 for (
int j = startCol; j < curCol; ++j) {
2416 blas.ROT (1, &R(j, curCol), LDR, &R(j+1, curCol), LDR,
2417 &theCosines[j], &theSines[j]);
2421 Scalar theCosine, theSine, result;
2422 computeGivensRotation (R(curCol, curCol), R(curCol+1, curCol),
2423 theCosine, theSine, result);
2424 theCosines[curCol] = theCosine;
2425 theSines[curCol] = theSine;
2430 R(curCol+1, curCol) = result;
2431 R(curCol+1, curCol) = STS::zero();
2438 const int LDZ = z.stride();
2439 blas.ROT (z.numCols(),
2440 &z(curCol,0), LDZ, &z(curCol+1,0), LDZ,
2441 &theCosine, &theSine);
2447 return STS::magnitude( z(numRows-1, 0) );
2453 #endif // __Belos_ProjectedLeastSquaresSolver_hpp Collection of types and exceptions used within the Belos solvers.
Teuchos::SerialDenseMatrix< int, Scalar > z
Current right-hand side of the projected least-squares problem.
Teuchos::SerialDenseMatrix< int, Scalar > y
Current solution of the projected least-squares problem.
Teuchos::Array< Scalar > theCosines
Array of cosines from the computed Givens rotations.
std::pair< int, bool > solveUpperTriangularSystemInPlace(Teuchos::ESide side, mat_type &X, const mat_type &R, const ERobustness robustness)
Solve square upper triangular linear system(s) in place.
ProjectedLeastSquaresSolver(std::ostream &warnStream, const ERobustness defaultRobustness=ROBUSTNESS_NONE)
Constructor.
void rightUpperTriSolve(mat_type &B, const mat_type &R) const
In Matlab notation: B = B / R, where R is upper triangular.
void axpy(mat_type &Y, const scalar_type &alpha, const mat_type &X) const
Y := Y + alpha * X.
std::string robustnessEnumToString(const ERobustness x)
Convert the given ERobustness enum value to a string.
std::pair< bool, std::pair< magnitude_type, magnitude_type > > isUpperTriangular(const mat_type &A) const
Is the matrix A upper triangular / trapezoidal?
ERobustness robustnessStringToEnum(const std::string &x)
Convert the given robustness string value to an ERobustness enum.
bool testGivensRotations(std::ostream &out)
Test Givens rotations.
std::pair< int, bool > solveUpperTriangularSystem(Teuchos::ESide side, mat_type &X, const mat_type &R, const mat_type &B, const ERobustness robustness)
Solve the given square upper triangular linear system(s).
void ensureUpperTriangular(const mat_type &A, const char *const matrixName) const
Throw an exception if A is not upper triangular / trapezoidal.
void conjugateTransposeOfUpperTriangular(mat_type &L, const mat_type &R) const
L := (conjugate) transpose of R (upper triangular).
Methods for solving GMRES' projected least-squares problem.
ProjectedLeastSquaresProblem(const int maxNumIterations)
Constructor.
void ensureEqualDimensions(const mat_type &A, const char *const matrixName, const int numRows, const int numCols) const
Ensure that the matrix A is exactly numRows by numCols.
magnitude_type updateColumns(ProjectedLeastSquaresProblem< Scalar > &problem, const int startCol, const int endCol)
Update columns [startCol,endCol] of the projected least-squares problem.
Namespace containing implementation details of Belos solvers.
Teuchos::SerialDenseMatrix< int, Scalar > H
The upper Hessenberg matrix from GMRES.
"Container" for the GMRES projected least-squares problem.
void ensureUpperHessenberg(const mat_type &A, const char *const matrixName) const
Throw an exception if A is not (strictly) upper Hessenberg.
void reset(const typename Teuchos::ScalarTraits< Scalar >::magnitudeType beta)
Reset the projected least-squares problem.
void ensureMinimumDimensions(const mat_type &A, const char *const matrixName, const int minNumRows, const int minNumCols) const
Ensure that the matrix A is at least minNumRows by minNumCols.
std::pair< int, bool > solveUpperTriangularSystemInPlace(Teuchos::ESide side, mat_type &X, const mat_type &R)
Solve square upper triangular linear system(s) in place.
Scalar scalar_type
The template parameter of this class.
std::pair< bool, std::pair< magnitude_type, magnitude_type > > isUpperHessenberg(const mat_type &A) const
Is the matrix A upper Hessenberg?
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
The type of the magnitude of scalar_type values.
void ensureUpperHessenberg(const mat_type &A, const char *const matrixName, const magnitude_type relativeTolerance) const
Throw an exception if A is not "approximately" upper Hessenberg.
Teuchos::RCP< Teuchos::ParameterEntryValidator > robustnessValidator()
Make a ParameterList validator for ERobustness.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
The type of the magnitude of a scalar_type value.
void matSub(mat_type &A, const mat_type &B) const
A := A - B.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
The type of a dense matrix (or vector) of scalar_type.
void solve(ProjectedLeastSquaresProblem< Scalar > &problem, const int curCol)
Solve the projected least-squares problem.
bool testUpdateColumn(std::ostream &out, const int numCols, const bool testBlockGivens=false, const bool extraVerbose=false)
Test update and solve using Givens rotations.
void zeroOutStrictLowerTriangle(mat_type &A) const
Zero out everything below the diagonal of A.
void conjugateTranspose(mat_type &A_star, const mat_type &A) const
A_star := (conjugate) transpose of A.
void matAdd(mat_type &A, const mat_type &B) const
A := A + B.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
The type of the magnitude of a scalar_type value.
Teuchos::Array< Scalar > theSines
Array of sines from the computed Givens rotations.
void matMatMult(const scalar_type &beta, mat_type &C, const scalar_type &alpha, const mat_type &A, const mat_type &B) const
C := beta*C + alpha*A*B.
Scalar scalar_type
The type of the entries in the projected least-squares problem.
Low-level operations on non-distributed dense matrices.
void reallocateAndReset(const typename Teuchos::ScalarTraits< Scalar >::magnitudeType beta, const int maxNumIterations)
(Re)allocate and reset the projected least-squares problem.
ERobustness
Robustness level of projected least-squares solver operations.
void matScale(mat_type &A, const scalar_type &alpha) const
A := alpha * A.
bool testTriangularSolves(std::ostream &out, const int testProblemSize, const ERobustness robustness, const bool verbose=false)
Test upper triangular solves.
std::pair< int, bool > solveUpperTriangularSystem(Teuchos::ESide side, mat_type &X, const mat_type &R, const mat_type &B)
Solve the given square upper triangular linear system(s).
int infNaNCount(const mat_type &A, const bool upperTriangular=false) const
Return the number of Inf or NaN entries in the matrix A.
void partition(Teuchos::RCP< mat_type > &A_11, Teuchos::RCP< mat_type > &A_21, Teuchos::RCP< mat_type > &A_12, Teuchos::RCP< mat_type > &A_22, mat_type &A, const int numRows1, const int numRows2, const int numCols1, const int numCols2)
A -> [A_11, A_21, A_12, A_22].
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::SerialDenseMatrix< int, Scalar > R
Upper triangular factor from the QR factorization of H.
magnitude_type updateColumn(ProjectedLeastSquaresProblem< Scalar > &problem, const int curCol)
Update column curCol of the projected least-squares problem.
Scalar scalar_type
The template parameter of this class.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
The type of a dense matrix (or vector) of scalar_type.