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) );