Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosProjectedLeastSquaresSolver.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef __Belos_ProjectedLeastSquaresSolver_hpp
11 #define __Belos_ProjectedLeastSquaresSolver_hpp
12 
13 #include "BelosConfigDefs.hpp"
14 #include "BelosTypes.hpp"
15 #include "Teuchos_Array.hpp"
16 #include "Teuchos_BLAS.hpp"
17 #include "Teuchos_LAPACK.hpp"
18 #include "Teuchos_oblackholestream.hpp"
19 #include "Teuchos_ScalarTraits.hpp"
21 #include "Teuchos_StandardParameterEntryValidators.hpp"
22 
26 
27 namespace Belos {
28 
37  namespace details {
38 
39  // Anonymous namespace restricts contents to file scope.
40  namespace {
43  template<class Scalar>
44  void
45  printMatrix (std::ostream& out,
46  const std::string& name,
48  {
49  using std::endl;
50 
51  const int numRows = A.numRows();
52  const int numCols = A.numCols();
53 
54  out << name << " = " << endl << '[';
55  if (numCols == 1) {
56  // Compact form for column vectors; valid Matlab.
57  for (int i = 0; i < numRows; ++i) {
58  out << A(i,0);
59  if (i < numRows-1) {
60  out << "; ";
61  }
62  }
63  } else {
64  for (int i = 0; i < numRows; ++i) {
65  for (int j = 0; j < numCols; ++j) {
66  out << A(i,j);
67  if (j < numCols-1) {
68  out << ", ";
69  } else if (i < numRows-1) {
70  out << ';' << endl;
71  }
72  }
73  }
74  }
75  out << ']' << endl;
76  }
77 
80  template<class Scalar>
81  void
82  print (std::ostream& out,
84  const std::string& linePrefix)
85  {
86  using std::endl;
87 
88  const int numRows = A.numRows();
89  const int numCols = A.numCols();
90 
91  out << linePrefix << '[';
92  if (numCols == 1) {
93  // Compact form for column vectors; valid Matlab.
94  for (int i = 0; i < numRows; ++i) {
95  out << A(i,0);
96  if (i < numRows-1) {
97  out << "; ";
98  }
99  }
100  } else {
101  for (int i = 0; i < numRows; ++i) {
102  for (int j = 0; j < numCols; ++j) {
103  if (numRows > 1) {
104  out << linePrefix << " ";
105  }
106  out << A(i,j);
107  if (j < numCols-1) {
108  out << ", ";
109  } else if (i < numRows-1) {
110  out << ';' << endl;
111  }
112  }
113  }
114  }
115  out << linePrefix << ']' << endl;
116  }
117  } // namespace (anonymous)
118 
126  template<class Scalar>
128  public:
131  typedef Scalar scalar_type;
135 
144 
157 
169 
179 
190 
196 
205  ProjectedLeastSquaresProblem (const int maxNumIterations) :
206  H (maxNumIterations+1, maxNumIterations),
207  R (maxNumIterations+1, maxNumIterations),
208  y (maxNumIterations+1, 1),
209  z (maxNumIterations+1, 1),
210  theCosines (maxNumIterations+1),
211  theSines (maxNumIterations+1)
212  {}
213 
232  void
234  {
235  typedef Teuchos::ScalarTraits<Scalar> STS;
236 
237  // Zero out the right-hand side of the least-squares problem.
238  z.putScalar (STS::zero());
239 
240  // Promote the initial residual norm from a magnitude type to
241  // a scalar type, so we can assign it to the first entry of z.
242  const Scalar initialResidualNorm (beta);
243  z(0,0) = initialResidualNorm;
244  }
245 
259  void
261  const int maxNumIterations)
262  {
263  typedef Teuchos::ScalarTraits<Scalar> STS;
265 
266  TEUCHOS_TEST_FOR_EXCEPTION(beta < STM::zero(), std::invalid_argument,
267  "ProjectedLeastSquaresProblem::reset: initial "
268  "residual beta = " << beta << " < 0.");
269  TEUCHOS_TEST_FOR_EXCEPTION(maxNumIterations <= 0, std::invalid_argument,
270  "ProjectedLeastSquaresProblem::reset: maximum number "
271  "of iterations " << maxNumIterations << " <= 0.");
272 
273  if (H.numRows() < maxNumIterations+1 || H.numCols() < maxNumIterations) {
274  const int errcode = H.reshape (maxNumIterations+1, maxNumIterations);
275  TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error,
276  "Failed to reshape H into a " << (maxNumIterations+1)
277  << " x " << maxNumIterations << " matrix.");
278  }
279  (void) H.putScalar (STS::zero());
280 
281  if (R.numRows() < maxNumIterations+1 || R.numCols() < maxNumIterations) {
282  const int errcode = R.reshape (maxNumIterations+1, maxNumIterations);
283  TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error,
284  "Failed to reshape R into a " << (maxNumIterations+1)
285  << " x " << maxNumIterations << " matrix.");
286  }
287  (void) R.putScalar (STS::zero());
288 
289  if (y.numRows() < maxNumIterations+1 || y.numCols() < 1) {
290  const int errcode = y.reshape (maxNumIterations+1, 1);
291  TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error,
292  "Failed to reshape y into a " << (maxNumIterations+1)
293  << " x " << 1 << " matrix.");
294  }
295  (void) y.putScalar (STS::zero());
296 
297  if (z.numRows() < maxNumIterations+1 || z.numCols() < 1) {
298  const int errcode = z.reshape (maxNumIterations+1, 1);
299  TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error,
300  "Failed to reshape z into a " << (maxNumIterations+1)
301  << " x " << 1 << " matrix.");
302  }
303  reset (beta);
304  }
305 
306  };
307 
308 
316  template<class Scalar>
318  public:
321  typedef Scalar scalar_type;
328 
329  private:
334 
335  public:
337  void
338  conjugateTranspose (mat_type& A_star, const mat_type& A) const
339  {
340  for (int i = 0; i < A.numRows(); ++i) {
341  for (int j = 0; j < A.numCols(); ++j) {
342  A_star(j,i) = STS::conjugate (A(i,j));
343  }
344  }
345  }
346 
348  void
350  {
351  const int N = R.numCols();
352 
353  for (int j = 0; j < N; ++j) {
354  for (int i = 0; i <= j; ++i) {
355  L(j,i) = STS::conjugate (R(i,j));
356  }
357  }
358  }
359 
361  void
363  {
364  const int N = std::min (A.numRows(), A.numCols());
365 
366  for (int j = 0; j < N; ++j) {
367  for (int i = j+1; i < A.numRows(); ++i) {
368  A(i,j) = STS::zero();
369  }
370  }
371  }
372 
381  void
386  mat_type& A,
387  const int numRows1,
388  const int numRows2,
389  const int numCols1,
390  const int numCols2)
391  {
392  using Teuchos::rcp;
393  using Teuchos::View;
394 
395  A_11 = rcp (new mat_type (View, A, numRows1, numCols1, 0, 0));
396  A_21 = rcp (new mat_type (View, A, numRows2, numCols1, numRows1, 0));
397  A_12 = rcp (new mat_type (View, A, numRows1, numCols2, 0, numCols1));
398  A_22 = rcp (new mat_type (View, A, numRows2, numCols2, numRows1, numCols1));
399  }
400 
402  void
403  matScale (mat_type& A, const scalar_type& alpha) const
404  {
405  // const int LDA = A.stride(); // unused
406  const int numRows = A.numRows();
407  const int numCols = A.numCols();
408 
409  if (numRows == 0 || numCols == 0) {
410  return;
411  } else {
412  for (int j = 0; j < numCols; ++j) {
413  scalar_type* const A_j = &A(0,j);
414 
415  for (int i = 0; i < numRows; ++i) {
416  A_j[i] *= alpha;
417  }
418  }
419  }
420  }
421 
426  void
428  const scalar_type& alpha,
429  const mat_type& X) const
430  {
431  const int numRows = Y.numRows();
432  const int numCols = Y.numCols();
433 
434  TEUCHOS_TEST_FOR_EXCEPTION(numRows != X.numRows() || numCols != X.numCols(),
435  std::invalid_argument, "Dimensions of X and Y don't "
436  "match. X is " << X.numRows() << " x " << X.numCols()
437  << ", and Y is " << numRows << " x " << numCols << ".");
438  for (int j = 0; j < numCols; ++j) {
439  for (int i = 0; i < numRows; ++i) {
440  Y(i,j) += alpha * X(i,j);
441  }
442  }
443  }
444 
446  void
447  matAdd (mat_type& A, const mat_type& B) const
448  {
449  const int numRows = A.numRows();
450  const int numCols = A.numCols();
451 
453  B.numRows() != numRows || B.numCols() != numCols,
454  std::invalid_argument,
455  "matAdd: The input matrices A and B have incompatible dimensions. "
456  "A is " << numRows << " x " << numCols << ", but B is " <<
457  B.numRows () << " x " << B.numCols () << ".");
458  if (numRows == 0 || numCols == 0) {
459  return;
460  } else {
461  for (int j = 0; j < numCols; ++j) {
462  scalar_type* const A_j = &A(0,j);
463  const scalar_type* const B_j = &B(0,j);
464 
465  for (int i = 0; i < numRows; ++i) {
466  A_j[i] += B_j[i];
467  }
468  }
469  }
470  }
471 
473  void
474  matSub (mat_type& A, const mat_type& B) const
475  {
476  const int numRows = A.numRows();
477  const int numCols = A.numCols();
478 
480  B.numRows() != numRows || B.numCols() != numCols,
481  std::invalid_argument,
482  "matSub: The input matrices A and B have incompatible dimensions. "
483  "A is " << numRows << " x " << numCols << ", but B is " <<
484  B.numRows () << " x " << B.numCols () << ".");
485  if (numRows == 0 || numCols == 0) {
486  return;
487  } else {
488  for (int j = 0; j < numCols; ++j) {
489  scalar_type* const A_j = &A(0,j);
490  const scalar_type* const B_j = &B(0,j);
491 
492  for (int i = 0; i < numRows; ++i) {
493  A_j[i] -= B_j[i];
494  }
495  }
496  }
497  }
498 
503  void
505  const mat_type& R) const
506  {
508  std::invalid_argument,
509  "rightUpperTriSolve: R and B have incompatible "
510  "dimensions. B has " << B.numCols() << " columns, "
511  "but R has " << R.numRows() << " rows.");
512  blas_type blas;
515  R.numCols(), B.numCols(),
516  STS::one(), R.values(), R.stride(),
517  B.values(), B.stride());
518  }
519 
525  void
526  matMatMult (const scalar_type& beta,
527  mat_type& C,
528  const scalar_type& alpha,
529  const mat_type& A,
530  const mat_type& B) const
531  {
532  using Teuchos::NO_TRANS;
533 
535  std::invalid_argument,
536  "matMatMult: The input matrices A and B have "
537  "incompatible dimensions. A is " << A.numRows()
538  << " x " << A.numCols() << ", but B is "
539  << B.numRows() << " x " << B.numCols() << ".");
541  std::invalid_argument,
542  "matMatMult: The input matrix A and the output "
543  "matrix C have incompatible dimensions. A has "
544  << A.numRows() << " rows, but C has " << C.numRows()
545  << " rows.");
547  std::invalid_argument,
548  "matMatMult: The input matrix B and the output "
549  "matrix C have incompatible dimensions. B has "
550  << B.numCols() << " columns, but C has "
551  << C.numCols() << " columns.");
552  blas_type blas;
553  blas.GEMM (NO_TRANS, NO_TRANS, C.numRows(), C.numCols(), A.numCols(),
554  alpha, A.values(), A.stride(), B.values(), B.stride(),
555  beta, C.values(), C.stride());
556  }
557 
565  int
566  infNaNCount (const mat_type& A, const bool upperTriangular=false) const
567  {
568  int count = 0;
569  for (int j = 0; j < A.numCols(); ++j) {
570  if (upperTriangular) {
571  for (int i = 0; i <= j && i < A.numRows(); ++i) {
572  if (STS::isnaninf (A(i,j))) {
573  ++count;
574  }
575  }
576  } else {
577  for (int i = 0; i < A.numRows(); ++i) {
578  if (STS::isnaninf (A(i,j))) {
579  ++count;
580  }
581  }
582  }
583  }
584  return count;
585  }
586 
592  std::pair<bool, std::pair<magnitude_type, magnitude_type> >
593  isUpperTriangular (const mat_type& A) const
594  {
595  magnitude_type lowerTri = STM::zero();
596  magnitude_type upperTri = STM::zero();
597  int count = 0;
598 
599  for (int j = 0; j < A.numCols(); ++j) {
600  // Compute the Frobenius norm of the upper triangle /
601  // trapezoid of A. The second clause of the loop upper
602  // bound is for matrices with fewer rows than columns.
603  for (int i = 0; i <= j && i < A.numRows(); ++i) {
604  const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
605  upperTri += A_ij_mag * A_ij_mag;
606  }
607  // Scan the strict lower triangle / trapezoid of A.
608  for (int i = j+1; i < A.numRows(); ++i) {
609  const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
610  lowerTri += A_ij_mag * A_ij_mag;
611  if (A_ij_mag != STM::zero()) {
612  ++count;
613  }
614  }
615  }
616  return std::make_pair (count == 0, std::make_pair (lowerTri, upperTri));
617  }
618 
619 
625  std::pair<bool, std::pair<magnitude_type, magnitude_type> >
626  isUpperHessenberg (const mat_type& A) const
627  {
628  magnitude_type lower = STM::zero();
629  magnitude_type upper = STM::zero();
630  int count = 0;
631 
632  for (int j = 0; j < A.numCols(); ++j) {
633  // Compute the Frobenius norm of the upper Hessenberg part
634  // of A. The second clause of the loop upper bound is for
635  // matrices with fewer rows than columns.
636  for (int i = 0; i <= j+1 && i < A.numRows(); ++i) {
637  const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
638  upper += A_ij_mag * A_ij_mag;
639  }
640  // Scan the strict lower part of A.
641  for (int i = j+2; i < A.numRows(); ++i) {
642  const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
643  lower += A_ij_mag * A_ij_mag;
644  if (A_ij_mag != STM::zero()) {
645  ++count;
646  }
647  }
648  }
649  return std::make_pair (count == 0, std::make_pair (lower, upper));
650  }
651 
658  void
660  const char* const matrixName) const
661  {
662  std::pair<bool, std::pair<magnitude_type, magnitude_type> > result =
663  isUpperTriangular (A);
664 
665  TEUCHOS_TEST_FOR_EXCEPTION(! result.first, std::invalid_argument,
666  "The " << A.numRows() << " x " << A.numCols()
667  << " matrix " << matrixName << " is not upper "
668  "triangular. ||tril(A)||_F = "
669  << result.second.first << " and ||A||_F = "
670  << result.second.second << ".");
671  }
672 
679  void
681  const char* const matrixName) const
682  {
683  std::pair<bool, std::pair<magnitude_type, magnitude_type> > result =
684  isUpperHessenberg (A);
685 
686  TEUCHOS_TEST_FOR_EXCEPTION(! result.first, std::invalid_argument,
687  "The " << A.numRows() << " x " << A.numCols()
688  << " matrix " << matrixName << " is not upper "
689  "triangular. ||tril(A(2:end, :))||_F = "
690  << result.second.first << " and ||A||_F = "
691  << result.second.second << ".");
692  }
693 
709  void
711  const char* const matrixName,
712  const magnitude_type relativeTolerance) const
713  {
714  std::pair<bool, std::pair<magnitude_type, magnitude_type> > result =
715  isUpperHessenberg (A);
716 
717  if (result.first) {
718  // Mollified relative departure from upper Hessenberg.
719  const magnitude_type err = (result.second.second == STM::zero() ?
720  result.second.first :
721  result.second.first / result.second.second);
722  TEUCHOS_TEST_FOR_EXCEPTION(err > relativeTolerance, std::invalid_argument,
723  "The " << A.numRows() << " x " << A.numCols()
724  << " matrix " << matrixName << " is not upper "
725  "triangular. ||tril(A(2:end, :))||_F "
726  << (result.second.second == STM::zero() ? "" : " / ||A||_F")
727  << " = " << err << " > " << relativeTolerance << ".");
728  }
729  }
730 
742  void
744  const char* const matrixName,
745  const int minNumRows,
746  const int minNumCols) const
747  {
748  TEUCHOS_TEST_FOR_EXCEPTION(A.numRows() < minNumRows || A.numCols() < minNumCols,
749  std::invalid_argument,
750  "The matrix " << matrixName << " is " << A.numRows()
751  << " x " << A.numCols() << ", and therefore does not "
752  "satisfy the minimum dimensions " << minNumRows
753  << " x " << minNumCols << ".");
754  }
755 
767  void
769  const char* const matrixName,
770  const int numRows,
771  const int numCols) const
772  {
773  TEUCHOS_TEST_FOR_EXCEPTION(A.numRows() != numRows || A.numCols() != numCols,
774  std::invalid_argument,
775  "The matrix " << matrixName << " is supposed to be "
776  << numRows << " x " << numCols << ", but is "
777  << A.numRows() << " x " << A.numCols() << " instead.");
778  }
779 
780  };
781 
800  enum ERobustness {
805  };
806 
808  inline std::string
810  {
811  const char* strings[] = {"None", "Some", "Lots"};
812  TEUCHOS_TEST_FOR_EXCEPTION(x < ROBUSTNESS_NONE || x >= ROBUSTNESS_INVALID,
813  std::invalid_argument,
814  "Invalid enum value " << x << ".");
815  return std::string (strings[x]);
816  }
817 
820  inline robustnessStringToEnum (const std::string& x)
821  {
822  const char* strings[] = {"None", "Some", "Lots"};
823  for (int r = 0; r < static_cast<int> (ROBUSTNESS_INVALID); ++r) {
824  if (x == strings[r]) {
825  return static_cast<ERobustness> (r);
826  }
827  }
828  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
829  "Invalid robustness string " << x << ".");
830  }
831 
844  {
845  using Teuchos::stringToIntegralParameterEntryValidator;
846 
852  docs[0] = "Use the BLAS' triangular solve. This may result in Inf or "
853  "NaN output if the triangular matrix is rank deficient.";
854  docs[1] = "Robustness somewhere between \"None\" and \"Lots\".";
855  docs[2] = "Solve the triangular system in a least-squares sense, using "
856  "an SVD-based algorithm. This will always succeed, though the "
857  "solution may not make sense for GMRES.";
859  ints[0] = ROBUSTNESS_NONE;
860  ints[1] = ROBUSTNESS_SOME;
861  ints[2] = ROBUSTNESS_LOTS;
862  const std::string pname ("Robustness of Projected Least-Squares Solve");
863 
864  return stringToIntegralParameterEntryValidator<ERobustness> (strs, docs,
865  ints, pname);
866  }
867 
930  template<class Scalar>
932  public:
939  typedef Scalar scalar_type;
948 
949  private:
954 
955  public:
969  ProjectedLeastSquaresSolver (std::ostream& warnStream,
970  const ERobustness defaultRobustness=ROBUSTNESS_NONE) :
971  warn_ (warnStream),
972  defaultRobustness_ (defaultRobustness)
973  {}
974 
990  const int curCol)
991  {
992  return updateColumnGivens (problem.H, problem.R, problem.y, problem.z,
993  problem.theCosines, problem.theSines, curCol);
994  }
995 
1012  const int startCol,
1013  const int endCol)
1014  {
1015  return updateColumnsGivens (problem.H, problem.R, problem.y, problem.z,
1016  problem.theCosines, problem.theSines,
1017  startCol, endCol);
1018  }
1019 
1033  void
1035  const int curCol)
1036  {
1037  solveGivens (problem.y, problem.R, problem.z, curCol);
1038  }
1039 
1044  std::pair<int, bool>
1046  mat_type& X,
1047  const mat_type& R,
1048  const mat_type& B)
1049  {
1050  return solveUpperTriangularSystem (side, X, R, B, defaultRobustness_);
1051  }
1052 
1091  std::pair<int, bool>
1093  mat_type& X,
1094  const mat_type& R,
1095  const mat_type& B,
1096  const ERobustness robustness)
1097  {
1098  TEUCHOS_TEST_FOR_EXCEPTION(X.numRows() != B.numRows(), std::invalid_argument,
1099  "The output X and right-hand side B have different "
1100  "numbers of rows. X has " << X.numRows() << " rows"
1101  ", and B has " << B.numRows() << " rows.");
1102  // If B has more columns than X, we ignore the remaining
1103  // columns of B when solving the upper triangular system. If
1104  // B has _fewer_ columns than X, we can't solve for all the
1105  // columns of X, so we throw an exception.
1106  TEUCHOS_TEST_FOR_EXCEPTION(X.numCols() > B.numCols(), std::invalid_argument,
1107  "The output X has more columns than the "
1108  "right-hand side B. X has " << X.numCols()
1109  << " columns and B has " << B.numCols()
1110  << " columns.");
1111  // See above explaining the number of columns in B_view.
1112  mat_type B_view (Teuchos::View, B, B.numRows(), X.numCols());
1113 
1114  // Both the BLAS' _TRSM and LAPACK's _LATRS overwrite the
1115  // right-hand side with the solution, so first copy B_view
1116  // into X.
1117  X.assign (B_view);
1118 
1119  // Solve the upper triangular system.
1120  return solveUpperTriangularSystemInPlace (side, X, R, robustness);
1121  }
1122 
1127  std::pair<int, bool>
1129  mat_type& X,
1130  const mat_type& R)
1131  {
1132  return solveUpperTriangularSystemInPlace (side, X, R, defaultRobustness_);
1133  }
1134 
1142  std::pair<int, bool>
1144  mat_type& X,
1145  const mat_type& R,
1146  const ERobustness robustness)
1147  {
1148  using Teuchos::Array;
1149  using Teuchos::Copy;
1150  using Teuchos::LEFT_SIDE;
1151  using Teuchos::RIGHT_SIDE;
1153 
1154  const int M = R.numRows();
1155  const int N = R.numCols();
1156  TEUCHOS_TEST_FOR_EXCEPTION(M < N, std::invalid_argument,
1157  "The input matrix R has fewer columns than rows. "
1158  "R is " << M << " x " << N << ".");
1159  // Ignore any additional rows of R by working with a square view.
1160  mat_type R_view (Teuchos::View, R, N, N);
1161 
1162  if (side == LEFT_SIDE) {
1163  TEUCHOS_TEST_FOR_EXCEPTION(X.numRows() < N, std::invalid_argument,
1164  "The input/output matrix X has only "
1165  << X.numRows() << " rows, but needs at least "
1166  << N << " rows to match the matrix for a "
1167  "left-side solve R \\ X.");
1168  } else if (side == RIGHT_SIDE) {
1169  TEUCHOS_TEST_FOR_EXCEPTION(X.numCols() < N, std::invalid_argument,
1170  "The input/output matrix X has only "
1171  << X.numCols() << " columns, but needs at least "
1172  << N << " columns to match the matrix for a "
1173  "right-side solve X / R.");
1174  }
1176  robustness >= ROBUSTNESS_INVALID,
1177  std::invalid_argument,
1178  "Invalid robustness value " << robustness << ".");
1179 
1180  // In robust mode, scan the matrix and right-hand side(s) for
1181  // Infs and NaNs. Only look at the upper triangle of the
1182  // matrix.
1183  if (robustness > ROBUSTNESS_NONE) {
1184  int count = ops.infNaNCount (R_view, true);
1185  TEUCHOS_TEST_FOR_EXCEPTION(count > 0, std::runtime_error,
1186  "There " << (count != 1 ? "are" : "is")
1187  << " " << count << " Inf or NaN entr"
1188  << (count != 1 ? "ies" : "y")
1189  << " in the upper triangle of R.");
1190  count = ops.infNaNCount (X, false);
1191  TEUCHOS_TEST_FOR_EXCEPTION(count > 0, std::runtime_error,
1192  "There " << (count != 1 ? "are" : "is")
1193  << " " << count << " Inf or NaN entr"
1194  << (count != 1 ? "ies" : "y") << " in the "
1195  "right-hand side(s) X.");
1196  }
1197 
1198  // Pair of values to return from this method.
1199  int rank = N;
1200  bool foundRankDeficiency = false;
1201 
1202  // Solve for X.
1203  blas_type blas;
1204 
1205  if (robustness == ROBUSTNESS_NONE) {
1206  // Fast triangular solve using the BLAS' _TRSM. This does
1207  // no checking for rank deficiency.
1210  STS::one(), R.values(), R.stride(),
1211  X.values(), X.stride());
1212  } else if (robustness < ROBUSTNESS_INVALID) {
1213  // Save a copy of X, since X contains the right-hand side on
1214  // input.
1215  mat_type B (Copy, X, X.numRows(), X.numCols());
1216 
1217  // Fast triangular solve using the BLAS' _TRSM. This does
1218  // no checking for rank deficiency.
1221  STS::one(), R.values(), R.stride(),
1222  X.values(), X.stride());
1223 
1224  // Check for Infs or NaNs in X. If there are any, then
1225  // assume that TRSM failed, and use a more robust algorithm.
1226  if (ops.infNaNCount (X, false) != 0) {
1227 
1228  warn_ << "Upper triangular solve: Found Infs and/or NaNs in the "
1229  "solution after using the fast algorithm. Retrying using a more "
1230  "robust algorithm." << std::endl;
1231 
1232  // Restore X from the copy.
1233  X.assign (B);
1234 
1235  // Find the minimum-norm solution to the least-squares
1236  // problem $\min_x \|RX - B\|_2$, using the singular value
1237  // decomposition (SVD).
1239  if (side == LEFT_SIDE) {
1240  // _GELSS overwrites its matrix input, so make a copy.
1241  mat_type R_copy (Teuchos::Copy, R_view, N, N);
1242 
1243  // Zero out the lower triangle of R_copy, since the
1244  // mat_type constructor copies all the entries, not just
1245  // the upper triangle. _GELSS will read all the entries
1246  // of the input matrix.
1247  ops.zeroOutStrictLowerTriangle (R_copy);
1248 
1249  // Solve the least-squares problem.
1250  rank = solveLeastSquaresUsingSVD (R_copy, X);
1251  } else {
1252  // If solving with R on the right-hand side, the interface
1253  // requires that instead of solving $\min \|XR - B\|_2$,
1254  // we have to solve $\min \|R^* X^* - B^*\|_2$. We
1255  // compute (conjugate) transposes in newly allocated
1256  // temporary matrices X_star resp. R_star. (B is already
1257  // in X and _GELSS overwrites its input vector X with the
1258  // solution.)
1259  mat_type X_star (X.numCols(), X.numRows());
1260  ops.conjugateTranspose (X_star, X);
1261  mat_type R_star (N, N); // Filled with zeros automatically.
1262  ops.conjugateTransposeOfUpperTriangular (R_star, R);
1263 
1264  // Solve the least-squares problem.
1265  rank = solveLeastSquaresUsingSVD (R_star, X_star);
1266 
1267  // Copy the transpose of X_star back into X.
1268  ops.conjugateTranspose (X, X_star);
1269  }
1270  if (rank < N) {
1271  foundRankDeficiency = true;
1272  }
1273  }
1274  } else {
1275  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
1276  "Should never get here! Invalid robustness value "
1277  << robustness << ". Please report this bug to the "
1278  "Belos developers.");
1279  }
1280  return std::make_pair (rank, foundRankDeficiency);
1281  }
1282 
1283 
1284  public:
1293  bool
1294  testGivensRotations (std::ostream& out)
1295  {
1296  using std::endl;
1297 
1298  out << "Testing Givens rotations:" << endl;
1299  Scalar x = STS::random();
1300  Scalar y = STS::random();
1301  out << " x = " << x << ", y = " << y << endl;
1302 
1303  Scalar theCosine, theSine, result;
1304  blas_type blas;
1305  computeGivensRotation (x, y, theCosine, theSine, result);
1306  out << "-- After computing rotation:" << endl;
1307  out << "---- cos,sin = " << theCosine << "," << theSine << endl;
1308  out << "---- x = " << x << ", y = " << y
1309  << ", result = " << result << endl;
1310 
1311  blas.ROT (1, &x, 1, &y, 1, &theCosine, &theSine);
1312  out << "-- After applying rotation:" << endl;
1313  out << "---- cos,sin = " << theCosine << "," << theSine << endl;
1314  out << "---- x = " << x << ", y = " << y << endl;
1315 
1316  // Allow only a tiny bit of wiggle room for zeroing-out of y.
1317  if (STS::magnitude(y) > 2*STS::eps())
1318  return false;
1319  else
1320  return true;
1321  }
1322 
1352  bool
1353  testUpdateColumn (std::ostream& out,
1354  const int numCols,
1355  const bool testBlockGivens=false,
1356  const bool extraVerbose=false)
1357  {
1358  using Teuchos::Array;
1359  using std::endl;
1360 
1361  TEUCHOS_TEST_FOR_EXCEPTION(numCols <= 0, std::invalid_argument,
1362  "numCols = " << numCols << " <= 0.");
1363  const int numRows = numCols + 1;
1364 
1365  mat_type H (numRows, numCols);
1366  mat_type z (numRows, 1);
1367 
1368  mat_type R_givens (numRows, numCols);
1369  mat_type y_givens (numRows, 1);
1370  mat_type z_givens (numRows, 1);
1371  Array<Scalar> theCosines (numCols);
1372  Array<Scalar> theSines (numCols);
1373 
1374  mat_type R_blockGivens (numRows, numCols);
1375  mat_type y_blockGivens (numRows, 1);
1376  mat_type z_blockGivens (numRows, 1);
1377  Array<Scalar> blockCosines (numCols);
1378  Array<Scalar> blockSines (numCols);
1379  const int panelWidth = std::min (3, numCols);
1380 
1381  mat_type R_lapack (numRows, numCols);
1382  mat_type y_lapack (numRows, 1);
1383  mat_type z_lapack (numRows, 1);
1384 
1385  // Make a random least-squares problem.
1386  makeRandomProblem (H, z);
1387  if (extraVerbose) {
1388  printMatrix<Scalar> (out, "H", H);
1389  printMatrix<Scalar> (out, "z", z);
1390  }
1391 
1392  // Set up the right-hand side copies for each of the methods.
1393  // Each method is free to overwrite its given right-hand side.
1394  z_givens.assign (z);
1395  if (testBlockGivens) {
1396  z_blockGivens.assign (z);
1397  }
1398  z_lapack.assign (z);
1399 
1400  //
1401  // Imitate how one would update the least-squares problem in a
1402  // typical GMRES implementation, for each updating method.
1403  //
1404  // Update using Givens rotations, one at a time.
1405  magnitude_type residualNormGivens = STM::zero();
1406  for (int curCol = 0; curCol < numCols; ++curCol) {
1407  residualNormGivens = updateColumnGivens (H, R_givens, y_givens, z_givens,
1408  theCosines, theSines, curCol);
1409  }
1410  solveGivens (y_givens, R_givens, z_givens, numCols-1);
1411 
1412  // Update using the "panel left-looking" Givens approach, with
1413  // the given panel width.
1414  magnitude_type residualNormBlockGivens = STM::zero();
1415  if (testBlockGivens) {
1416  const bool testBlocksAtATime = true;
1417  if (testBlocksAtATime) {
1418  // Blocks of columns at a time.
1419  for (int startCol = 0; startCol < numCols; startCol += panelWidth) {
1420  int endCol = std::min (startCol + panelWidth - 1, numCols - 1);
1421  residualNormBlockGivens =
1422  updateColumnsGivens (H, R_blockGivens, y_blockGivens, z_blockGivens,
1423  blockCosines, blockSines, startCol, endCol);
1424  }
1425  } else {
1426  // One column at a time. This is good as a sanity check
1427  // to make sure updateColumnsGivens() with a single column
1428  // does the same thing as updateColumnGivens().
1429  for (int startCol = 0; startCol < numCols; ++startCol) {
1430  residualNormBlockGivens =
1431  updateColumnsGivens (H, R_blockGivens, y_blockGivens, z_blockGivens,
1432  blockCosines, blockSines, startCol, startCol);
1433  }
1434  }
1435  // The panel version of Givens should compute the same
1436  // cosines and sines as the non-panel version, and should
1437  // update the right-hand side z in the same way. Thus, we
1438  // should be able to use the same triangular solver.
1439  solveGivens (y_blockGivens, R_blockGivens, z_blockGivens, numCols-1);
1440  }
1441 
1442  // Solve using LAPACK's least-squares solver.
1443  const magnitude_type residualNormLapack =
1444  solveLapack (H, R_lapack, y_lapack, z_lapack, numCols-1);
1445 
1446  // Compute the condition number of the least-squares problem.
1447  // This requires a residual, so use the residual from the
1448  // LAPACK method. All that the method needs for an accurate
1449  // residual norm is forward stability.
1450  const magnitude_type leastSquaresCondNum =
1451  leastSquaresConditionNumber (H, z, residualNormLapack);
1452 
1453  // Compute the relative least-squares solution error for both
1454  // Givens methods. We assume that the LAPACK solution is
1455  // "exact" and compare against the Givens rotations solution.
1456  // This is taking liberties with the definition of condition
1457  // number, but it's the best we can do, since we don't know
1458  // the exact solution and don't have an extended-precision
1459  // solver.
1460 
1461  // The solution lives only in y[0 .. numCols-1].
1462  mat_type y_givens_view (Teuchos::View, y_givens, numCols, 1);
1463  mat_type y_blockGivens_view (Teuchos::View, y_blockGivens, numCols, 1);
1464  mat_type y_lapack_view (Teuchos::View, y_lapack, numCols, 1);
1465 
1466  const magnitude_type givensSolutionError =
1467  solutionError (y_givens_view, y_lapack_view);
1468  const magnitude_type blockGivensSolutionError = testBlockGivens ?
1469  solutionError (y_blockGivens_view, y_lapack_view) :
1470  STM::zero();
1471 
1472  // If printing out the matrices, copy out the upper triangular
1473  // factors for printing. (Both methods are free to leave data
1474  // below the lower triangle.)
1475  if (extraVerbose) {
1476  mat_type R_factorFromGivens (numCols, numCols);
1477  mat_type R_factorFromBlockGivens (numCols, numCols);
1478  mat_type R_factorFromLapack (numCols, numCols);
1479 
1480  for (int j = 0; j < numCols; ++j) {
1481  for (int i = 0; i <= j; ++i) {
1482  R_factorFromGivens(i,j) = R_givens(i,j);
1483  if (testBlockGivens) {
1484  R_factorFromBlockGivens(i,j) = R_blockGivens(i,j);
1485  }
1486  R_factorFromLapack(i,j) = R_lapack(i,j);
1487  }
1488  }
1489 
1490  printMatrix<Scalar> (out, "R_givens", R_factorFromGivens);
1491  printMatrix<Scalar> (out, "y_givens", y_givens_view);
1492  printMatrix<Scalar> (out, "z_givens", z_givens);
1493 
1494  if (testBlockGivens) {
1495  printMatrix<Scalar> (out, "R_blockGivens", R_factorFromBlockGivens);
1496  printMatrix<Scalar> (out, "y_blockGivens", y_blockGivens_view);
1497  printMatrix<Scalar> (out, "z_blockGivens", z_blockGivens);
1498  }
1499 
1500  printMatrix<Scalar> (out, "R_lapack", R_factorFromLapack);
1501  printMatrix<Scalar> (out, "y_lapack", y_lapack_view);
1502  printMatrix<Scalar> (out, "z_lapack", z_lapack);
1503  }
1504 
1505  // Compute the (Frobenius) norm of the original matrix H.
1506  const magnitude_type H_norm = H.normFrobenius();
1507 
1508  out << "||H||_F = " << H_norm << endl;
1509 
1510  out << "||H y_givens - z||_2 / ||H||_F = "
1511  << leastSquaresResidualNorm (H, y_givens_view, z) / H_norm << endl;
1512  if (testBlockGivens) {
1513  out << "||H y_blockGivens - z||_2 / ||H||_F = "
1514  << leastSquaresResidualNorm (H, y_blockGivens_view, z) / H_norm << endl;
1515  }
1516  out << "||H y_lapack - z||_2 / ||H||_F = "
1517  << leastSquaresResidualNorm (H, y_lapack_view, z) / H_norm << endl;
1518 
1519  out << "||y_givens - y_lapack||_2 / ||y_lapack||_2 = "
1520  << givensSolutionError << endl;
1521  if (testBlockGivens) {
1522  out << "||y_blockGivens - y_lapack||_2 / ||y_lapack||_2 = "
1523  << blockGivensSolutionError << endl;
1524  }
1525 
1526  out << "Least-squares condition number = "
1527  << leastSquaresCondNum << endl;
1528 
1529  // Now for the controversial part of the test: judging whether
1530  // we succeeded. This includes the problem's condition
1531  // number, which is a measure of the maximum perturbation in
1532  // the solution for which we can still say that the solution
1533  // is valid. We include a little wiggle room by including a
1534  // factor proportional to the square root of the number of
1535  // floating-point operations that influence the last entry
1536  // (the conventional Wilkinsonian heuristic), times 10 for
1537  // good measure.
1538  //
1539  // (The square root looks like it has something to do with an
1540  // average-case probabilistic argument, but doesn't really.
1541  // What's an "average problem"?)
1542  const magnitude_type wiggleFactor =
1543  10 * STM::squareroot( numRows*numCols );
1544  const magnitude_type solutionErrorBoundFactor =
1545  wiggleFactor * leastSquaresCondNum;
1546  const magnitude_type solutionErrorBound =
1547  solutionErrorBoundFactor * STS::eps();
1548  out << "Solution error bound: " << solutionErrorBoundFactor
1549  << " * eps = " << solutionErrorBound << endl;
1550 
1551  // Remember that NaN is not greater than, not less than, and
1552  // not equal to any other number, including itself. Some
1553  // compilers will rudely optimize away the "x != x" test.
1554  if (STM::isnaninf (solutionErrorBound)) {
1555  // Hm, the solution error bound is Inf or NaN. This
1556  // probably means that the test problem was generated
1557  // incorrectly. We should return false in this case.
1558  return false;
1559  } else { // solution error bound is finite.
1560  if (STM::isnaninf (givensSolutionError)) {
1561  return false;
1562  } else if (givensSolutionError > solutionErrorBound) {
1563  return false;
1564  } else if (testBlockGivens) {
1565  if (STM::isnaninf (blockGivensSolutionError)) {
1566  return false;
1567  } else if (blockGivensSolutionError > solutionErrorBound) {
1568  return false;
1569  } else { // Givens and Block Givens tests succeeded.
1570  return true;
1571  }
1572  } else { // Not testing block Givens; Givens test succeeded.
1573  return true;
1574  }
1575  }
1576  }
1577 
1593  bool
1594  testTriangularSolves (std::ostream& out,
1595  const int testProblemSize,
1596  const ERobustness robustness,
1597  const bool verbose=false)
1598  {
1599  using Teuchos::LEFT_SIDE;
1600  using Teuchos::RIGHT_SIDE;
1601  using std::endl;
1603 
1604  Teuchos::oblackholestream blackHole;
1605  std::ostream& verboseOut = verbose ? out : blackHole;
1606 
1607  verboseOut << "Testing upper triangular solves" << endl;
1608  //
1609  // Construct an upper triangular linear system to solve.
1610  //
1611  verboseOut << "-- Generating test matrix" << endl;
1612  const int N = testProblemSize;
1613  mat_type R (N, N);
1614  // Fill the upper triangle of R with random numbers.
1615  for (int j = 0; j < N; ++j) {
1616  for (int i = 0; i <= j; ++i) {
1617  R(i,j) = STS::random ();
1618  }
1619  }
1620  // Compute the Frobenius norm of R for later use.
1621  const magnitude_type R_norm = R.normFrobenius ();
1622  // Fill the right-hand side B with random numbers.
1623  mat_type B (N, 1);
1624  B.random ();
1625  // Compute the Frobenius norm of B for later use.
1626  const magnitude_type B_norm = B.normFrobenius ();
1627 
1628  // Save a copy of the original upper triangular system.
1629  mat_type R_copy (Teuchos::Copy, R, N, N);
1630  mat_type B_copy (Teuchos::Copy, B, N, 1);
1631 
1632  // Solution vector.
1633  mat_type X (N, 1);
1634 
1635  // Solve RX = B.
1636  verboseOut << "-- Solving RX=B" << endl;
1637  // We're ignoring the return values for now.
1638  (void) solveUpperTriangularSystem (LEFT_SIDE, X, R, B, robustness);
1639  // Test the residual error.
1640  mat_type Resid (N, 1);
1641  Resid.assign (B_copy);
1643  ops.matMatMult (STS::one(), Resid, -STS::one(), R_copy, X);
1644  verboseOut << "---- ||R*X - B||_F = " << Resid.normFrobenius() << endl;
1645  verboseOut << "---- ||R||_F ||X||_F + ||B||_F = "
1646  << (R_norm * X.normFrobenius() + B_norm)
1647  << endl;
1648 
1649  // Restore R and B.
1650  R.assign (R_copy);
1651  B.assign (B_copy);
1652 
1653  //
1654  // Set up a right-side test problem: YR = B^*.
1655  //
1656  mat_type Y (1, N);
1657  mat_type B_star (1, N);
1658  ops.conjugateTranspose (B_star, B);
1659  mat_type B_star_copy (1, N);
1660  B_star_copy.assign (B_star);
1661  // Solve YR = B^*.
1662  verboseOut << "-- Solving YR=B^*" << endl;
1663  // We're ignoring the return values for now.
1664  (void) solveUpperTriangularSystem (RIGHT_SIDE, Y, R, B_star, robustness);
1665  // Test the residual error.
1666  mat_type Resid2 (1, N);
1667  Resid2.assign (B_star_copy);
1668  ops.matMatMult (STS::one(), Resid2, -STS::one(), Y, R_copy);
1669  verboseOut << "---- ||Y*R - B^*||_F = " << Resid2.normFrobenius() << endl;
1670  verboseOut << "---- ||Y||_F ||R||_F + ||B^*||_F = "
1671  << (Y.normFrobenius() * R_norm + B_norm)
1672  << endl;
1673 
1674  // FIXME (mfh 14 Oct 2011) The test always "passes" for now;
1675  // you have to inspect the output in order to see whether it
1676  // succeeded. We really should fix the above to use the
1677  // infinity-norm bounds in Higham's book for triangular
1678  // solves. That would automate the test.
1679  return true;
1680  }
1681 
1682  private:
1684  std::ostream& warn_;
1685 
1690  ERobustness defaultRobustness_;
1691 
1692  private:
1704  int
1705  solveLeastSquaresUsingSVD (mat_type& A, mat_type& X)
1706  {
1707  using Teuchos::Array;
1709 
1710  if (defaultRobustness_ > ROBUSTNESS_SOME) {
1711  int count = ops.infNaNCount (A);
1712  TEUCHOS_TEST_FOR_EXCEPTION(count != 0, std::invalid_argument,
1713  "solveLeastSquaresUsingSVD: The input matrix A "
1714  "contains " << count << "Inf and/or NaN entr"
1715  << (count != 1 ? "ies" : "y") << ".");
1716  count = ops.infNaNCount (X);
1717  TEUCHOS_TEST_FOR_EXCEPTION(count != 0, std::invalid_argument,
1718  "solveLeastSquaresUsingSVD: The input matrix X "
1719  "contains " << count << "Inf and/or NaN entr"
1720  << (count != 1 ? "ies" : "y") << ".");
1721  }
1722  const int N = std::min (A.numRows(), A.numCols());
1723  lapack_type lapack;
1724 
1725  // Rank of A; to be computed by _GELSS and returned.
1726  int rank = N;
1727 
1728  // Use Scalar's machine precision for the rank tolerance,
1729  // not magnitude_type's machine precision.
1730  const magnitude_type rankTolerance = STS::eps();
1731 
1732  // Array of singular values.
1733  Array<magnitude_type> singularValues (N);
1734 
1735  // Extra workspace. This is only used by _GELSS if Scalar is
1736  // complex. Teuchos::LAPACK presents a unified interface to
1737  // _GELSS that always includes the RWORK argument, even though
1738  // LAPACK's SGELSS and DGELSS don't have the RWORK argument.
1739  // We always allocate at least one entry so that &rwork[0]
1740  // makes sense.
1741  Array<magnitude_type> rwork (1);
1742  if (STS::isComplex) {
1743  rwork.resize (std::max (1, 5 * N));
1744  }
1745  //
1746  // Workspace query
1747  //
1748  Scalar lworkScalar = STS::one(); // To be set by workspace query
1749  int info = 0;
1750  lapack.GELSS (A.numRows(), A.numCols(), X.numCols(),
1751  A.values(), A.stride(), X.values(), X.stride(),
1752  &singularValues[0], rankTolerance, &rank,
1753  &lworkScalar, -1, &rwork[0], &info);
1754  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1755  "_GELSS workspace query returned INFO = "
1756  << info << " != 0.");
1757  const int lwork = static_cast<int> (STS::real (lworkScalar));
1758  TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
1759  "_GELSS workspace query returned LWORK = "
1760  << lwork << " < 0.");
1761  // Allocate workspace. Size > 0 means &work[0] makes sense.
1762  Array<Scalar> work (std::max (1, lwork));
1763  // Solve the least-squares problem.
1764  lapack.GELSS (A.numRows(), A.numCols(), X.numCols(),
1765  A.values(), A.stride(), X.values(), X.stride(),
1766  &singularValues[0], rankTolerance, &rank,
1767  &work[0], lwork, &rwork[0], &info);
1768  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
1769  "_GELSS returned INFO = " << info << " != 0.");
1770  return rank;
1771  }
1772 
1779  void
1780  solveGivens (mat_type& y,
1781  mat_type& R,
1782  const mat_type& z,
1783  const int curCol)
1784  {
1785  const int numRows = curCol + 2;
1786 
1787  // Now that we have the updated R factor of H, and the updated
1788  // right-hand side z, solve the least-squares problem by
1789  // solving the upper triangular linear system Ry=z for y.
1790  const mat_type R_view (Teuchos::View, R, numRows-1, numRows-1);
1791  const mat_type z_view (Teuchos::View, z, numRows-1, z.numCols());
1792  mat_type y_view (Teuchos::View, y, numRows-1, y.numCols());
1793 
1795  R_view, z_view, defaultRobustness_);
1796  }
1797 
1799  void
1800  makeRandomProblem (mat_type& H, mat_type& z)
1801  {
1802  // In GMRES, z always starts out with only the first entry
1803  // being nonzero. That entry always has nonnegative real part
1804  // and zero imaginary part, since it is the initial residual
1805  // norm.
1806  H.random ();
1807  // Zero out the entries below the subdiagonal of H, so that it
1808  // is upper Hessenberg.
1809  for (int j = 0; j < H.numCols(); ++j) {
1810  for (int i = j+2; i < H.numRows(); ++i) {
1811  H(i,j) = STS::zero();
1812  }
1813  }
1814  // Initialize z, the right-hand side of the least-squares
1815  // problem. Make the first entry of z nonzero.
1816  {
1817  // It's still possible that a random number will come up
1818  // zero after 1000 trials, but unlikely. Nevertheless, it's
1819  // still important not to allow an infinite loop, for
1820  // example if the pseudorandom number generator is broken
1821  // and always returns zero.
1822  const int numTrials = 1000;
1823  magnitude_type z_init = STM::zero();
1824  for (int trial = 0; trial < numTrials && z_init == STM::zero(); ++trial) {
1825  z_init = STM::random();
1826  }
1827  TEUCHOS_TEST_FOR_EXCEPTION(z_init == STM::zero(), std::runtime_error,
1828  "After " << numTrials << " trial"
1829  << (numTrials != 1 ? "s" : "")
1830  << ", we were unable to generate a nonzero pseudo"
1831  "random real number. This most likely indicates a "
1832  "broken pseudorandom number generator.");
1833  const magnitude_type z_first = (z_init < 0) ? -z_init : z_init;
1834 
1835  // NOTE I'm assuming here that "scalar_type = magnitude_type"
1836  // assignments make sense.
1837  z(0,0) = z_first;
1838  }
1839  }
1840 
1844  void
1845  computeGivensRotation (const Scalar& x,
1846  const Scalar& y,
1847  Scalar& theCosine,
1848  Scalar& theSine,
1849  Scalar& result)
1850  {
1851  // _LARTG, an LAPACK aux routine, is slower but more accurate
1852  // than the BLAS' _ROTG.
1853  const bool useLartg = false;
1854 
1855  if (useLartg) {
1856  lapack_type lapack;
1857  // _LARTG doesn't clobber its input arguments x and y.
1858  lapack.LARTG (x, y, &theCosine, &theSine, &result);
1859  } else {
1860  // _ROTG clobbers its first two arguments. x is overwritten
1861  // with the result of applying the Givens rotation: [x; y] ->
1862  // [x (on output); 0]. y is overwritten with the "fast"
1863  // Givens transform (see Golub and Van Loan, 3rd ed.).
1864  Scalar x_temp = x;
1865  Scalar y_temp = y;
1866  blas_type blas;
1867  blas.ROTG (&x_temp, &y_temp, &theCosine, &theSine);
1868  result = x_temp;
1869  }
1870  }
1871 
1873  void
1874  singularValues (const mat_type& A,
1876  {
1877  using Teuchos::Array;
1878  using Teuchos::ArrayView;
1879 
1880  const int numRows = A.numRows();
1881  const int numCols = A.numCols();
1882  TEUCHOS_TEST_FOR_EXCEPTION(sigmas.size() < std::min (numRows, numCols),
1883  std::invalid_argument,
1884  "The sigmas array is only of length " << sigmas.size()
1885  << ", but must be of length at least "
1886  << std::min (numRows, numCols)
1887  << " in order to hold all the singular values of the "
1888  "matrix A.");
1889 
1890  // Compute the condition number of the matrix A, using a singular
1891  // value decomposition (SVD). LAPACK's SVD routine overwrites the
1892  // input matrix, so make a copy.
1893  mat_type A_copy (numRows, numCols);
1894  A_copy.assign (A);
1895 
1896  // Workspace query.
1897  lapack_type lapack;
1898  int info = 0;
1899  Scalar lworkScalar = STS::zero();
1900  Array<magnitude_type> rwork (std::max (std::min (numRows, numCols) - 1, 1));
1901  lapack.GESVD ('N', 'N', numRows, numCols,
1902  A_copy.values(), A_copy.stride(), &sigmas[0],
1903  (Scalar*) NULL, 1, (Scalar*) NULL, 1,
1904  &lworkScalar, -1, &rwork[0], &info);
1905 
1906  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1907  "LAPACK _GESVD workspace query failed with INFO = "
1908  << info << ".");
1909  const int lwork = static_cast<int> (STS::real (lworkScalar));
1910  TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
1911  "LAPACK _GESVD workspace query returned LWORK = "
1912  << lwork << " < 0.");
1913  // Make sure that the workspace array always has positive
1914  // length, so that &work[0] makes sense.
1915  Teuchos::Array<Scalar> work (std::max (1, lwork));
1916 
1917  // Compute the singular values of A.
1918  lapack.GESVD ('N', 'N', numRows, numCols,
1919  A_copy.values(), A_copy.stride(), &sigmas[0],
1920  (Scalar*) NULL, 1, (Scalar*) NULL, 1,
1921  &work[0], lwork, &rwork[0], &info);
1922  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1923  "LAPACK _GESVD failed with INFO = " << info << ".");
1924  }
1925 
1933  std::pair<magnitude_type, magnitude_type>
1934  extremeSingularValues (const mat_type& A)
1935  {
1936  using Teuchos::Array;
1937 
1938  const int numRows = A.numRows();
1939  const int numCols = A.numCols();
1940 
1941  Array<magnitude_type> sigmas (std::min (numRows, numCols));
1942  singularValues (A, sigmas);
1943  return std::make_pair (sigmas[0], sigmas[std::min(numRows, numCols) - 1]);
1944  }
1945 
1956  magnitude_type
1957  leastSquaresConditionNumber (const mat_type& A,
1958  const mat_type& b,
1959  const magnitude_type& residualNorm)
1960  {
1961  // Extreme singular values of A.
1962  const std::pair<magnitude_type, magnitude_type> sigmaMaxMin =
1963  extremeSingularValues (A);
1964 
1965  // Our solvers currently assume that H has full rank. If the
1966  // test matrix doesn't have full rank, we stop right away.
1967  TEUCHOS_TEST_FOR_EXCEPTION(sigmaMaxMin.second == STM::zero(), std::runtime_error,
1968  "The test matrix is rank deficient; LAPACK's _GESVD "
1969  "routine reports that its smallest singular value is "
1970  "zero.");
1971  // 2-norm condition number of A. We checked above that the
1972  // denominator is nonzero.
1973  const magnitude_type A_cond = sigmaMaxMin.first / sigmaMaxMin.second;
1974 
1975  // "Theta" in the variable names below refers to the angle between
1976  // the vectors b and A*x, where x is the computed solution. It
1977  // measures whether the residual norm is large (near ||b||) or
1978  // small (near 0).
1979  const magnitude_type sinTheta = residualNorm / b.normFrobenius();
1980 
1981  // \sin^2 \theta + \cos^2 \theta = 1.
1982  //
1983  // The range of sine is [-1,1], so squaring it won't overflow.
1984  // We still have to check whether sinTheta > 1, though. This
1985  // is impossible in exact arithmetic, assuming that the
1986  // least-squares solver worked (b-A*0 = b and x minimizes
1987  // ||b-A*x||_2, so ||b-A*0||_2 >= ||b-A*x||_2). However, it
1988  // might just be possible in floating-point arithmetic. We're
1989  // just looking for an estimate, so if sinTheta > 1, we cap it
1990  // at 1.
1991  const magnitude_type cosTheta = (sinTheta > STM::one()) ?
1992  STM::zero() : STM::squareroot (1 - sinTheta * sinTheta);
1993 
1994  // This may result in Inf, if cosTheta is zero. That's OK; in
1995  // that case, the condition number of the (full-rank)
1996  // least-squares problem is rightfully infinite.
1997  const magnitude_type tanTheta = sinTheta / cosTheta;
1998 
1999  // Condition number for the full-rank least-squares problem.
2000  return 2 * A_cond / cosTheta + tanTheta * A_cond * A_cond;
2001  }
2002 
2004  magnitude_type
2005  leastSquaresResidualNorm (const mat_type& A,
2006  const mat_type& x,
2007  const mat_type& b)
2008  {
2009  mat_type r (b.numRows(), b.numCols());
2010 
2011  // r := b - A*x
2012  r.assign (b);
2013  LocalDenseMatrixOps<Scalar> ops;
2014  ops.matMatMult (STS::one(), r, -STS::one(), A, x);
2015  return r.normFrobenius ();
2016  }
2017 
2022  magnitude_type
2023  solutionError (const mat_type& x_approx,
2024  const mat_type& x_exact)
2025  {
2026  const int numRows = x_exact.numRows();
2027  const int numCols = x_exact.numCols();
2028 
2029  mat_type x_diff (numRows, numCols);
2030  for (int j = 0; j < numCols; ++j) {
2031  for (int i = 0; i < numRows; ++i) {
2032  x_diff(i,j) = x_exact(i,j) - x_approx(i,j);
2033  }
2034  }
2035  const magnitude_type scalingFactor = x_exact.normFrobenius();
2036 
2037  // If x_exact has zero norm, just use the absolute difference.
2038  return x_diff.normFrobenius() /
2039  (scalingFactor == STM::zero() ? STM::one() : scalingFactor);
2040  }
2041 
2088  magnitude_type
2089  updateColumnGivens (const mat_type& H,
2090  mat_type& R,
2091  mat_type& y,
2092  mat_type& z,
2095  const int curCol)
2096  {
2097  using std::cerr;
2098  using std::endl;
2099 
2100  const int numRows = curCol + 2; // curCol is zero-based
2101  const int LDR = R.stride();
2102  const bool extraDebug = false;
2103 
2104  if (extraDebug) {
2105  cerr << "updateColumnGivens: curCol = " << curCol << endl;
2106  }
2107 
2108  // View of H( 1:curCol+1, curCol ) (in Matlab notation, if
2109  // curCol were a one-based index, as it would be in Matlab but
2110  // is not here).
2111  const mat_type H_col (Teuchos::View, H, numRows, 1, 0, curCol);
2112 
2113  // View of R( 1:curCol+1, curCol ) (again, in Matlab notation,
2114  // if curCol were a one-based index).
2115  mat_type R_col (Teuchos::View, R, numRows, 1, 0, curCol);
2116 
2117  // 1. Copy the current column from H into R, where it will be
2118  // modified.
2119  R_col.assign (H_col);
2120 
2121  if (extraDebug) {
2122  printMatrix<Scalar> (cerr, "R_col before ", R_col);
2123  }
2124 
2125  // 2. Apply all the previous Givens rotations, if any, to the
2126  // current column of the matrix.
2127  blas_type blas;
2128  for (int j = 0; j < curCol; ++j) {
2129  // BLAS::ROT really should take "const Scalar*" for these
2130  // arguments, but it wants a "Scalar*" instead, alas.
2131  Scalar theCosine = theCosines[j];
2132  Scalar theSine = theSines[j];
2133 
2134  if (extraDebug) {
2135  cerr << " j = " << j << ": (cos,sin) = "
2136  << theCosines[j] << "," << theSines[j] << endl;
2137  }
2138  blas.ROT (1, &R_col(j,0), LDR, &R_col(j+1,0), LDR,
2139  &theCosine, &theSine);
2140  }
2141  if (extraDebug && curCol > 0) {
2142  printMatrix<Scalar> (cerr, "R_col after applying previous "
2143  "Givens rotations", R_col);
2144  }
2145 
2146  // 3. Calculate new Givens rotation for R(curCol, curCol),
2147  // R(curCol+1, curCol).
2148  Scalar theCosine, theSine, result;
2149  computeGivensRotation (R_col(curCol,0), R_col(curCol+1,0),
2150  theCosine, theSine, result);
2151  theCosines[curCol] = theCosine;
2152  theSines[curCol] = theSine;
2153 
2154  if (extraDebug) {
2155  cerr << " New cos,sin = " << theCosine << "," << theSine << endl;
2156  }
2157 
2158  // 4. _Apply_ the new Givens rotation. We don't need to
2159  // invoke _ROT here, because computeGivensRotation()
2160  // already gives us the result: [x; y] -> [result; 0].
2161  R_col(curCol, 0) = result;
2162  R_col(curCol+1, 0) = STS::zero();
2163 
2164  if (extraDebug) {
2165  printMatrix<Scalar> (cerr, "R_col after applying current "
2166  "Givens rotation", R_col);
2167  }
2168 
2169  // 5. Apply the resulting Givens rotation to z (the right-hand
2170  // side of the projected least-squares problem).
2171  //
2172  // We prefer overgeneralization to undergeneralization by assuming
2173  // here that z may have more than one column.
2174  const int LDZ = z.stride();
2175  blas.ROT (z.numCols(),
2176  &z(curCol,0), LDZ, &z(curCol+1,0), LDZ,
2177  &theCosine, &theSine);
2178 
2179  if (extraDebug) {
2180  //mat_type R_after (Teuchos::View, R, numRows, numRows-1);
2181  //printMatrix<Scalar> (cerr, "R_after", R_after);
2182  //mat_type z_after (Teuchos::View, z, numRows, z.numCols());
2183  printMatrix<Scalar> (cerr, "z_after", z);
2184  }
2185 
2186  // The last entry of z is the nonzero part of the residual of the
2187  // least-squares problem. Its magnitude gives the residual 2-norm
2188  // of the least-squares problem.
2189  return STS::magnitude( z(numRows-1, 0) );
2190  }
2191 
2225  magnitude_type
2226  solveLapack (const mat_type& H,
2227  mat_type& R,
2228  mat_type& y,
2229  mat_type& z,
2230  const int curCol)
2231  {
2232  const int numRows = curCol + 2;
2233  const int numCols = curCol + 1;
2234  const int LDR = R.stride();
2235 
2236  // Copy H( 1:curCol+1, 1:curCol ) into R( 1:curCol+1, 1:curCol ).
2237  const mat_type H_view (Teuchos::View, H, numRows, numCols);
2238  mat_type R_view (Teuchos::View, R, numRows, numCols);
2239  R_view.assign (H_view);
2240 
2241  // The LAPACK least-squares solver overwrites the right-hand side
2242  // vector with the solution, so first copy z into y.
2243  mat_type y_view (Teuchos::View, y, numRows, y.numCols());
2244  mat_type z_view (Teuchos::View, z, numRows, y.numCols());
2245  y_view.assign (z_view);
2246 
2247  // Workspace query for the least-squares routine.
2248  int info = 0;
2249  Scalar lworkScalar = STS::zero();
2250  lapack_type lapack;
2251  lapack.GELS ('N', numRows, numCols, y_view.numCols(),
2252  NULL, LDR, NULL, y_view.stride(),
2253  &lworkScalar, -1, &info);
2254  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
2255  "LAPACK _GELS workspace query failed with INFO = "
2256  << info << ", for a " << numRows << " x " << numCols
2257  << " matrix with " << y_view.numCols()
2258  << " right hand side"
2259  << ((y_view.numCols() != 1) ? "s" : "") << ".");
2261  std::logic_error,
2262  "LAPACK _GELS workspace query returned an LWORK with "
2263  "negative real part: LWORK = " << lworkScalar
2264  << ". That should never happen. Please report this "
2265  "to the Belos developers.");
2267  std::logic_error,
2268  "LAPACK _GELS workspace query returned an LWORK with "
2269  "nonzero imaginary part: LWORK = " << lworkScalar
2270  << ". That should never happen. Please report this "
2271  "to the Belos developers.");
2272  // Cast workspace from Scalar to int. Scalar may be complex,
2273  // hence the request for the real part. Don't ask for the
2274  // magnitude, since computing the magnitude may overflow due
2275  // to squaring and square root to int. Hopefully LAPACK
2276  // doesn't ever overflow int this way.
2277  const int lwork = std::max (1, static_cast<int> (STS::real (lworkScalar)));
2278 
2279  // Allocate workspace for solving the least-squares problem.
2280  Teuchos::Array<Scalar> work (lwork);
2281 
2282  // Solve the least-squares problem. The ?: operator prevents
2283  // accessing the first element of the work array, if it has
2284  // length zero.
2285  lapack.GELS ('N', numRows, numCols, y_view.numCols(),
2286  R_view.values(), R_view.stride(),
2287  y_view.values(), y_view.stride(),
2288  (lwork > 0 ? &work[0] : (Scalar*) NULL),
2289  lwork, &info);
2290 
2291  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
2292  "Solving projected least-squares problem with LAPACK "
2293  "_GELS failed with INFO = " << info << ", for a "
2294  << numRows << " x " << numCols << " matrix with "
2295  << y_view.numCols() << " right hand side"
2296  << (y_view.numCols() != 1 ? "s" : "") << ".");
2297  // Extract the projected least-squares problem's residual error.
2298  // It's the magnitude of the last entry of y_view on output from
2299  // LAPACK's least-squares solver.
2300  return STS::magnitude( y_view(numRows-1, 0) );
2301  }
2302 
2313  magnitude_type
2314  updateColumnsGivens (const mat_type& H,
2315  mat_type& R,
2316  mat_type& y,
2317  mat_type& z,
2320  const int startCol,
2321  const int endCol)
2322  {
2323  TEUCHOS_TEST_FOR_EXCEPTION(startCol > endCol, std::invalid_argument,
2324  "updateColumnGivens: startCol = " << startCol
2325  << " > endCol = " << endCol << ".");
2326  magnitude_type lastResult = STM::zero();
2327  // [startCol, endCol] is an inclusive range.
2328  for (int curCol = startCol; curCol <= endCol; ++curCol) {
2329  lastResult = updateColumnGivens (H, R, y, z, theCosines, theSines, curCol);
2330  }
2331  return lastResult;
2332  }
2333 
2346  magnitude_type
2347  updateColumnsGivensBlock (const mat_type& H,
2348  mat_type& R,
2349  mat_type& y,
2350  mat_type& z,
2353  const int startCol,
2354  const int endCol)
2355  {
2356  const int numRows = endCol + 2;
2357  const int numColsToUpdate = endCol - startCol + 1;
2358  const int LDR = R.stride();
2359 
2360  // 1. Copy columns [startCol, endCol] from H into R, where they
2361  // will be modified.
2362  {
2363  const mat_type H_view (Teuchos::View, H, numRows, numColsToUpdate, 0, startCol);
2364  mat_type R_view (Teuchos::View, R, numRows, numColsToUpdate, 0, startCol);
2365  R_view.assign (H_view);
2366  }
2367 
2368  // 2. Apply all the previous Givens rotations, if any, to
2369  // columns [startCol, endCol] of the matrix. (Remember
2370  // that we're using a left-looking QR factorization
2371  // approach; we haven't yet touched those columns.)
2372  blas_type blas;
2373  for (int j = 0; j < startCol; ++j) {
2374  blas.ROT (numColsToUpdate,
2375  &R(j, startCol), LDR, &R(j+1, startCol), LDR,
2376  &theCosines[j], &theSines[j]);
2377  }
2378 
2379  // 3. Update each column in turn of columns [startCol, endCol].
2380  for (int curCol = startCol; curCol < endCol; ++curCol) {
2381  // a. Apply the Givens rotations computed in previous
2382  // iterations of this loop to the current column of R.
2383  for (int j = startCol; j < curCol; ++j) {
2384  blas.ROT (1, &R(j, curCol), LDR, &R(j+1, curCol), LDR,
2385  &theCosines[j], &theSines[j]);
2386  }
2387  // b. Calculate new Givens rotation for R(curCol, curCol),
2388  // R(curCol+1, curCol).
2389  Scalar theCosine, theSine, result;
2390  computeGivensRotation (R(curCol, curCol), R(curCol+1, curCol),
2391  theCosine, theSine, result);
2392  theCosines[curCol] = theCosine;
2393  theSines[curCol] = theSine;
2394 
2395  // c. _Apply_ the new Givens rotation. We don't need to
2396  // invoke _ROT here, because computeGivensRotation()
2397  // already gives us the result: [x; y] -> [result; 0].
2398  R(curCol+1, curCol) = result;
2399  R(curCol+1, curCol) = STS::zero();
2400 
2401  // d. Apply the resulting Givens rotation to z (the right-hand
2402  // side of the projected least-squares problem).
2403  //
2404  // We prefer overgeneralization to undergeneralization by
2405  // assuming here that z may have more than one column.
2406  const int LDZ = z.stride();
2407  blas.ROT (z.numCols(),
2408  &z(curCol,0), LDZ, &z(curCol+1,0), LDZ,
2409  &theCosine, &theSine);
2410  }
2411 
2412  // The last entry of z is the nonzero part of the residual of the
2413  // least-squares problem. Its magnitude gives the residual 2-norm
2414  // of the least-squares problem.
2415  return STS::magnitude( z(numRows-1, 0) );
2416  }
2417  }; // class ProjectedLeastSquaresSolver
2418  } // namespace details
2419 } // namespace Belos
2420 
2421 #endif // __Belos_ProjectedLeastSquaresSolver_hpp
ScalarType * values() const
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.
void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb) const
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 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, const ERobustness robustness)
Solve square upper triangular linear system(s) in place.
ProjectedLeastSquaresSolver(std::ostream &warnStream, const ERobustness defaultRobustness=ROBUSTNESS_NONE)
Constructor.
void matSub(mat_type &A, const mat_type &B) const
A := A - B.
static magnitudeType eps()
void ROT(const OrdinalType &n, ScalarType *dx, const OrdinalType &incx, ScalarType *dy, const OrdinalType &incy, MagnitudeType *c, ScalarType *s) const
void conjugateTransposeOfUpperTriangular(mat_type &L, const mat_type &R) const
L := (conjugate) transpose of R (upper triangular).
int infNaNCount(const mat_type &A, const bool upperTriangular=false) const
Return the number of Inf or NaN entries in the matrix A.
static T squareroot(T x)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static magnitudeType real(T a)
size_type size() const
std::string robustnessEnumToString(const ERobustness x)
Convert the given ERobustness enum value to a string.
ERobustness robustnessStringToEnum(const std::string &x)
Convert the given robustness string value to an ERobustness enum.
void GEMM(ETransp transa, ETransp transb, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
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).
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Methods for solving GMRES&#39; projected least-squares problem.
static T conjugate(T a)
void zeroOutStrictLowerTriangle(mat_type &A) const
Zero out everything below the diagonal of A.
ProjectedLeastSquaresProblem(const int maxNumIterations)
Constructor.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
magnitude_type updateColumns(ProjectedLeastSquaresProblem< Scalar > &problem, const int startCol, const int endCol)
Update columns [startCol,endCol] of the projected least-squares problem.
Teuchos::SerialDenseMatrix< int, Scalar > H
The upper Hessenberg matrix from GMRES.
&quot;Container&quot; for the GMRES projected least-squares problem.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
void reset(const typename Teuchos::ScalarTraits< Scalar >::magnitudeType beta)
Reset the projected least-squares problem.
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.
void matScale(mat_type &A, const scalar_type &alpha) const
A := alpha * A.
void axpy(mat_type &Y, const scalar_type &alpha, const mat_type &X) const
Y := Y + alpha * X.
static magnitudeType imag(T a)
static bool isnaninf(const T &x)
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
The type of the magnitude of scalar_type values.
Teuchos::RCP< Teuchos::ParameterEntryValidator > robustnessValidator()
Make a ParameterList validator for ERobustness.
void conjugateTranspose(mat_type &A_star, const mat_type &A) const
A_star := (conjugate) transpose of A.
static magnitudeType magnitude(T a)
OrdinalType numCols() const
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
The type of the magnitude of a scalar_type value.
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.
int reshape(OrdinalType numRows, OrdinalType numCols)
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.
Scalar scalar_type
The type of the entries in 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?
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.
void matAdd(mat_type &A, const mat_type &B) const
A := A + B.
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.
bool testTriangularSolves(std::ostream &out, const int testProblemSize, const ERobustness robustness, const bool verbose=false)
Test upper triangular solves.
static const bool isComplex
void ensureUpperHessenberg(const mat_type &A, const char *const matrixName) const
Throw an exception if A is not (strictly) upper Hessenberg.
std::pair< bool, std::pair< magnitude_type, magnitude_type > > isUpperTriangular(const mat_type &A) const
Is the matrix A upper triangular / trapezoidal?
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).
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 -&gt; [A_11, A_21, A_12, A_22].
Belos header file which uses auto-configuration information to include necessary C++ headers...
void rightUpperTriSolve(mat_type &B, const mat_type &R) const
In Matlab notation: B = B / R, where R is upper triangular.
void ensureUpperTriangular(const mat_type &A, const char *const matrixName) const
Throw an exception if A is not upper triangular / trapezoidal.
Teuchos::SerialDenseMatrix< int, Scalar > R
Upper triangular factor from the QR factorization of H.
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
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.
OrdinalType stride() const
magnitude_type updateColumn(ProjectedLeastSquaresProblem< Scalar > &problem, const int curCol)
Update column curCol of the projected least-squares problem.
OrdinalType numRows() const
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.
void ensureUpperHessenberg(const mat_type &A, const char *const matrixName, const magnitude_type relativeTolerance) const
Throw an exception if A is not &quot;approximately&quot; upper Hessenberg.

Generated on Wed Dec 11 2024 09:24:06 for Belos by doxygen 1.8.5