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;
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;
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 <<
".");
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());
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);
1291 mat_type X_star (X.numCols(), X.numRows());
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;
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);
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) {
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 =
1455 blockCosines, blockSines, startCol, endCol);
1461 for (
int startCol = 0; startCol < numCols; ++startCol) {
1462 residualNormBlockGivens =
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);
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);
1500 const magnitude_type blockGivensSolutionError = testBlockGivens ?
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 = " 1544 if (testBlockGivens) {
1545 out <<
"||H y_blockGivens - z||_2 / ||H||_F = " 1548 out <<
"||H y_lapack - z||_2 / ||H||_F = " 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 ();
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)
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());
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(),
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(),
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());
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.");
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>
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));
1975 return std::make_pair (sigmas[0], sigmas[std::min(numRows, numCols) - 1]);
1994 const std::pair<magnitude_type, magnitude_type> sigmaMaxMin =
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;
2041 mat_type r (b.numRows(), b.numCols());
2046 ops.
matMatMult (STS::one(), r, -STS::one(), A, x);
2047 return r.normFrobenius ();
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);
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;
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) );
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) {
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;
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.
void computeGivensRotation(const Scalar &x, const Scalar &y, Scalar &theCosine, Scalar &theSine, Scalar &result)
Compute the Givens rotation corresponding to [x; y].
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.
Teuchos::ScalarTraits< scalar_type > STS
void axpy(mat_type &Y, const scalar_type &alpha, const mat_type &X) const
Y := Y + alpha * X.
magnitude_type updateColumnGivens(const mat_type &H, mat_type &R, mat_type &y, mat_type &z, Teuchos::ArrayView< scalar_type > theCosines, Teuchos::ArrayView< scalar_type > theSines, const int curCol)
Update current column using Givens rotations.
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.
std::pair< magnitude_type, magnitude_type > extremeSingularValues(const mat_type &A)
The (largest, smallest) singular values of the given matrix.
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.
Teuchos::ScalarTraits< scalar_type > STS
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.
std::ostream & warn_
Stream to which to output warnings.
ERobustness defaultRobustness_
Default robustness level, for things like triangular solves.
magnitude_type updateColumnsGivens(const mat_type &H, mat_type &R, mat_type &y, mat_type &z, Teuchos::ArrayView< scalar_type > theCosines, Teuchos::ArrayView< scalar_type > theSines, const int startCol, const int endCol)
Update columns [startCol,endCol] of the projected least-squares problem.
Teuchos::BLAS< int, scalar_type > blas_type
Teuchos::LAPACK< int, scalar_type > lapack_type
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 solveGivens(mat_type &y, mat_type &R, const mat_type &z, const int curCol)
Solve the projected least-squares problem, assuming Givens rotations updates.
void makeRandomProblem(mat_type &H, mat_type &z)
Make a random 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.
Teuchos::ScalarTraits< magnitude_type > STM
magnitude_type updateColumnsGivensBlock(const mat_type &H, mat_type &R, mat_type &y, mat_type &z, Teuchos::ArrayView< scalar_type > theCosines, Teuchos::ArrayView< scalar_type > theSines, const int startCol, const int endCol)
Update columns [startCol,endCol] of the projected least-squares problem.
std::pair< bool, std::pair< magnitude_type, magnitude_type > > isUpperHessenberg(const mat_type &A) const
Is the matrix A upper Hessenberg?
magnitude_type solutionError(const mat_type &x_approx, const mat_type &x_exact)
||x_approx - x_exact||_2 // ||x_exact||_2.
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.
magnitude_type leastSquaresConditionNumber(const mat_type &A, const mat_type &b, const magnitude_type &residualNorm)
Normwise 2-norm condition number of the least-squares problem.
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.
magnitude_type solveLapack(const mat_type &H, mat_type &R, mat_type &y, mat_type &z, const int curCol)
Solve the least-squares problem using LAPACK.
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.
int solveLeastSquaresUsingSVD(mat_type &A, mat_type &X)
Solve using the SVD.
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.
Teuchos::ScalarTraits< magnitude_type > STM
Teuchos::BLAS< int, scalar_type > blas_type
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 leastSquaresResidualNorm(const mat_type &A, const mat_type &x, const mat_type &b)
(Frobenius norm if b has more than one column).
void singularValues(const mat_type &A, Teuchos::ArrayView< magnitude_type > sigmas)
Compute the singular values of A. Store them in the sigmas array.
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.
Teuchos::LAPACK< int, scalar_type > lapack_type