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 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef __Belos_ProjectedLeastSquaresSolver_hpp
43 #define __Belos_ProjectedLeastSquaresSolver_hpp
44 
45 #include "BelosConfigDefs.hpp"
46 #include "BelosTypes.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"
53 #include "Teuchos_StandardParameterEntryValidators.hpp"
54 
58 
59 namespace Belos {
60 
69  namespace details {
70 
71  // Anonymous namespace restricts contents to file scope.
72  namespace {
75  template<class Scalar>
76  void
77  printMatrix (std::ostream& out,
78  const std::string& name,
80  {
81  using std::endl;
82 
83  const int numRows = A.numRows();
84  const int numCols = A.numCols();
85 
86  out << name << " = " << endl << '[';
87  if (numCols == 1) {
88  // Compact form for column vectors; valid Matlab.
89  for (int i = 0; i < numRows; ++i) {
90  out << A(i,0);
91  if (i < numRows-1) {
92  out << "; ";
93  }
94  }
95  } else {
96  for (int i = 0; i < numRows; ++i) {
97  for (int j = 0; j < numCols; ++j) {
98  out << A(i,j);
99  if (j < numCols-1) {
100  out << ", ";
101  } else if (i < numRows-1) {
102  out << ';' << endl;
103  }
104  }
105  }
106  }
107  out << ']' << endl;
108  }
109 
112  template<class Scalar>
113  void
114  print (std::ostream& out,
116  const std::string& linePrefix)
117  {
118  using std::endl;
119 
120  const int numRows = A.numRows();
121  const int numCols = A.numCols();
122 
123  out << linePrefix << '[';
124  if (numCols == 1) {
125  // Compact form for column vectors; valid Matlab.
126  for (int i = 0; i < numRows; ++i) {
127  out << A(i,0);
128  if (i < numRows-1) {
129  out << "; ";
130  }
131  }
132  } else {
133  for (int i = 0; i < numRows; ++i) {
134  for (int j = 0; j < numCols; ++j) {
135  if (numRows > 1) {
136  out << linePrefix << " ";
137  }
138  out << A(i,j);
139  if (j < numCols-1) {
140  out << ", ";
141  } else if (i < numRows-1) {
142  out << ';' << endl;
143  }
144  }
145  }
146  }
147  out << linePrefix << ']' << endl;
148  }
149  } // namespace (anonymous)
150 
158  template<class Scalar>
160  public:
163  typedef Scalar scalar_type;
167 
176 
189 
201 
211 
222 
228 
237  ProjectedLeastSquaresProblem (const int maxNumIterations) :
238  H (maxNumIterations+1, maxNumIterations),
239  R (maxNumIterations+1, maxNumIterations),
240  y (maxNumIterations+1, 1),
241  z (maxNumIterations+1, 1),
242  theCosines (maxNumIterations+1),
243  theSines (maxNumIterations+1)
244  {}
245 
264  void
266  {
267  typedef Teuchos::ScalarTraits<Scalar> STS;
268 
269  // Zero out the right-hand side of the least-squares problem.
270  z.putScalar (STS::zero());
271 
272  // Promote the initial residual norm from a magnitude type to
273  // a scalar type, so we can assign it to the first entry of z.
274  const Scalar initialResidualNorm (beta);
275  z(0,0) = initialResidualNorm;
276  }
277 
291  void
293  const int maxNumIterations)
294  {
295  typedef Teuchos::ScalarTraits<Scalar> STS;
297 
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.");
304 
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.");
310  }
311  (void) H.putScalar (STS::zero());
312 
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.");
318  }
319  (void) R.putScalar (STS::zero());
320 
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.");
326  }
327  (void) y.putScalar (STS::zero());
328 
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.");
334  }
335  reset (beta);
336  }
337 
338  };
339 
340 
348  template<class Scalar>
350  public:
353  typedef Scalar scalar_type;
360 
361  private:
366 
367  public:
369  void
370  conjugateTranspose (mat_type& A_star, const mat_type& A) const
371  {
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));
375  }
376  }
377  }
378 
380  void
382  {
383  const int N = R.numCols();
384 
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));
388  }
389  }
390  }
391 
393  void
395  {
396  const int N = std::min (A.numRows(), A.numCols());
397 
398  for (int j = 0; j < N; ++j) {
399  for (int i = j+1; i < A.numRows(); ++i) {
400  A(i,j) = STS::zero();
401  }
402  }
403  }
404 
413  void
418  mat_type& A,
419  const int numRows1,
420  const int numRows2,
421  const int numCols1,
422  const int numCols2)
423  {
424  using Teuchos::rcp;
425  using Teuchos::View;
426 
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));
431  }
432 
434  void
435  matScale (mat_type& A, const scalar_type& alpha) const
436  {
437  // const int LDA = A.stride(); // unused
438  const int numRows = A.numRows();
439  const int numCols = A.numCols();
440 
441  if (numRows == 0 || numCols == 0) {
442  return;
443  } else {
444  for (int j = 0; j < numCols; ++j) {
445  scalar_type* const A_j = &A(0,j);
446 
447  for (int i = 0; i < numRows; ++i) {
448  A_j[i] *= alpha;
449  }
450  }
451  }
452  }
453 
458  void
460  const scalar_type& alpha,
461  const mat_type& X) const
462  {
463  const int numRows = Y.numRows();
464  const int numCols = Y.numCols();
465 
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);
473  }
474  }
475  }
476 
478  void
479  matAdd (mat_type& A, const mat_type& B) const
480  {
481  const int numRows = A.numRows();
482  const int numCols = A.numCols();
483 
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) {
491  return;
492  } else {
493  for (int j = 0; j < numCols; ++j) {
494  scalar_type* const A_j = &A(0,j);
495  const scalar_type* const B_j = &B(0,j);
496 
497  for (int i = 0; i < numRows; ++i) {
498  A_j[i] += B_j[i];
499  }
500  }
501  }
502  }
503 
505  void
506  matSub (mat_type& A, const mat_type& B) const
507  {
508  const int numRows = A.numRows();
509  const int numCols = A.numCols();
510 
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) {
518  return;
519  } else {
520  for (int j = 0; j < numCols; ++j) {
521  scalar_type* const A_j = &A(0,j);
522  const scalar_type* const B_j = &B(0,j);
523 
524  for (int i = 0; i < numRows; ++i) {
525  A_j[i] -= B_j[i];
526  }
527  }
528  }
529  }
530 
535  void
537  const mat_type& R) const
538  {
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.");
544  blas_type blas;
547  R.numCols(), B.numCols(),
548  STS::one(), R.values(), R.stride(),
549  B.values(), B.stride());
550  }
551 
557  void
558  matMatMult (const scalar_type& beta,
559  mat_type& C,
560  const scalar_type& alpha,
561  const mat_type& A,
562  const mat_type& B) const
563  {
564  using Teuchos::NO_TRANS;
565 
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() << ".");
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()
577  << " rows.");
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.");
584  blas_type blas;
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());
588  }
589 
597  int
598  infNaNCount (const mat_type& A, const bool upperTriangular=false) const
599  {
600  int count = 0;
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))) {
605  ++count;
606  }
607  }
608  } else {
609  for (int i = 0; i < A.numRows(); ++i) {
610  if (STS::isnaninf (A(i,j))) {
611  ++count;
612  }
613  }
614  }
615  }
616  return count;
617  }
618 
624  std::pair<bool, std::pair<magnitude_type, magnitude_type> >
625  isUpperTriangular (const mat_type& A) const
626  {
627  magnitude_type lowerTri = STM::zero();
628  magnitude_type upperTri = STM::zero();
629  int count = 0;
630 
631  for (int j = 0; j < A.numCols(); ++j) {
632  // Compute the Frobenius norm of the upper triangle /
633  // trapezoid of A. The second clause of the loop upper
634  // bound is for matrices with fewer rows than columns.
635  for (int i = 0; i <= j && i < A.numRows(); ++i) {
636  const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
637  upperTri += A_ij_mag * A_ij_mag;
638  }
639  // Scan the strict lower triangle / trapezoid of A.
640  for (int i = j+1; i < A.numRows(); ++i) {
641  const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
642  lowerTri += A_ij_mag * A_ij_mag;
643  if (A_ij_mag != STM::zero()) {
644  ++count;
645  }
646  }
647  }
648  return std::make_pair (count == 0, std::make_pair (lowerTri, upperTri));
649  }
650 
651 
657  std::pair<bool, std::pair<magnitude_type, magnitude_type> >
658  isUpperHessenberg (const mat_type& A) const
659  {
660  magnitude_type lower = STM::zero();
661  magnitude_type upper = STM::zero();
662  int count = 0;
663 
664  for (int j = 0; j < A.numCols(); ++j) {
665  // Compute the Frobenius norm of the upper Hessenberg part
666  // of A. The second clause of the loop upper bound is for
667  // matrices with fewer rows than columns.
668  for (int i = 0; i <= j+1 && i < A.numRows(); ++i) {
669  const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
670  upper += A_ij_mag * A_ij_mag;
671  }
672  // Scan the strict lower part of A.
673  for (int i = j+2; i < A.numRows(); ++i) {
674  const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
675  lower += A_ij_mag * A_ij_mag;
676  if (A_ij_mag != STM::zero()) {
677  ++count;
678  }
679  }
680  }
681  return std::make_pair (count == 0, std::make_pair (lower, upper));
682  }
683 
690  void
692  const char* const matrixName) const
693  {
694  std::pair<bool, std::pair<magnitude_type, magnitude_type> > result =
695  isUpperTriangular (A);
696 
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 << ".");
703  }
704 
711  void
713  const char* const matrixName) const
714  {
715  std::pair<bool, std::pair<magnitude_type, magnitude_type> > result =
716  isUpperHessenberg (A);
717 
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 << ".");
724  }
725 
741  void
743  const char* const matrixName,
744  const magnitude_type relativeTolerance) const
745  {
746  std::pair<bool, std::pair<magnitude_type, magnitude_type> > result =
747  isUpperHessenberg (A);
748 
749  if (result.first) {
750  // Mollified relative departure from upper Hessenberg.
751  const magnitude_type err = (result.second.second == STM::zero() ?
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 << ".");
760  }
761  }
762 
774  void
776  const char* const matrixName,
777  const int minNumRows,
778  const int minNumCols) const
779  {
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 << ".");
786  }
787 
799  void
801  const char* const matrixName,
802  const int numRows,
803  const int numCols) const
804  {
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.");
810  }
811 
812  };
813 
832  enum ERobustness {
837  };
838 
840  inline std::string
842  {
843  const char* strings[] = {"None", "Some", "Lots"};
844  TEUCHOS_TEST_FOR_EXCEPTION(x < ROBUSTNESS_NONE || x >= ROBUSTNESS_INVALID,
845  std::invalid_argument,
846  "Invalid enum value " << x << ".");
847  return std::string (strings[x]);
848  }
849 
852  inline robustnessStringToEnum (const std::string& x)
853  {
854  const char* strings[] = {"None", "Some", "Lots"};
855  for (int r = 0; r < static_cast<int> (ROBUSTNESS_INVALID); ++r) {
856  if (x == strings[r]) {
857  return static_cast<ERobustness> (r);
858  }
859  }
860  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
861  "Invalid robustness string " << x << ".");
862  }
863 
876  {
877  using Teuchos::stringToIntegralParameterEntryValidator;
878 
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.";
891  ints[0] = ROBUSTNESS_NONE;
892  ints[1] = ROBUSTNESS_SOME;
893  ints[2] = ROBUSTNESS_LOTS;
894  const std::string pname ("Robustness of Projected Least-Squares Solve");
895 
896  return stringToIntegralParameterEntryValidator<ERobustness> (strs, docs,
897  ints, pname);
898  }
899 
962  template<class Scalar>
964  public:
971  typedef Scalar scalar_type;
980 
981  private:
986 
987  public:
1001  ProjectedLeastSquaresSolver (std::ostream& warnStream,
1002  const ERobustness defaultRobustness=ROBUSTNESS_NONE) :
1003  warn_ (warnStream),
1004  defaultRobustness_ (defaultRobustness)
1005  {}
1006 
1022  const int curCol)
1023  {
1024  return updateColumnGivens (problem.H, problem.R, problem.y, problem.z,
1025  problem.theCosines, problem.theSines, curCol);
1026  }
1027 
1044  const int startCol,
1045  const int endCol)
1046  {
1047  return updateColumnsGivens (problem.H, problem.R, problem.y, problem.z,
1048  problem.theCosines, problem.theSines,
1049  startCol, endCol);
1050  }
1051 
1065  void
1067  const int curCol)
1068  {
1069  solveGivens (problem.y, problem.R, problem.z, curCol);
1070  }
1071 
1076  std::pair<int, bool>
1078  mat_type& X,
1079  const mat_type& R,
1080  const mat_type& B)
1081  {
1082  return solveUpperTriangularSystem (side, X, R, B, defaultRobustness_);
1083  }
1084 
1123  std::pair<int, bool>
1125  mat_type& X,
1126  const mat_type& R,
1127  const mat_type& B,
1128  const ERobustness robustness)
1129  {
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.");
1134  // If B has more columns than X, we ignore the remaining
1135  // columns of B when solving the upper triangular system. If
1136  // B has _fewer_ columns than X, we can't solve for all the
1137  // columns of X, so we throw an exception.
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()
1142  << " columns.");
1143  // See above explaining the number of columns in B_view.
1144  mat_type B_view (Teuchos::View, B, B.numRows(), X.numCols());
1145 
1146  // Both the BLAS' _TRSM and LAPACK's _LATRS overwrite the
1147  // right-hand side with the solution, so first copy B_view
1148  // into X.
1149  X.assign (B_view);
1150 
1151  // Solve the upper triangular system.
1152  return solveUpperTriangularSystemInPlace (side, X, R, robustness);
1153  }
1154 
1159  std::pair<int, bool>
1161  mat_type& X,
1162  const mat_type& R)
1163  {
1164  return solveUpperTriangularSystemInPlace (side, X, R, defaultRobustness_);
1165  }
1166 
1174  std::pair<int, bool>
1176  mat_type& X,
1177  const mat_type& R,
1178  const ERobustness robustness)
1179  {
1180  using Teuchos::Array;
1181  using Teuchos::Copy;
1182  using Teuchos::LEFT_SIDE;
1183  using Teuchos::RIGHT_SIDE;
1185 
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 << ".");
1191  // Ignore any additional rows of R by working with a square view.
1192  mat_type R_view (Teuchos::View, R, N, N);
1193 
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.");
1206  }
1208  robustness >= ROBUSTNESS_INVALID,
1209  std::invalid_argument,
1210  "Invalid robustness value " << robustness << ".");
1211 
1212  // In robust mode, scan the matrix and right-hand side(s) for
1213  // Infs and NaNs. Only look at the upper triangle of the
1214  // matrix.
1215  if (robustness > ROBUSTNESS_NONE) {
1216  int count = ops.infNaNCount (R_view, true);
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.");
1222  count = ops.infNaNCount (X, false);
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.");
1228  }
1229 
1230  // Pair of values to return from this method.
1231  int rank = N;
1232  bool foundRankDeficiency = false;
1233 
1234  // Solve for X.
1235  blas_type blas;
1236 
1237  if (robustness == ROBUSTNESS_NONE) {
1238  // Fast triangular solve using the BLAS' _TRSM. This does
1239  // no checking for rank deficiency.
1242  STS::one(), R.values(), R.stride(),
1243  X.values(), X.stride());
1244  } else if (robustness < ROBUSTNESS_INVALID) {
1245  // Save a copy of X, since X contains the right-hand side on
1246  // input.
1247  mat_type B (Copy, X, X.numRows(), X.numCols());
1248 
1249  // Fast triangular solve using the BLAS' _TRSM. This does
1250  // no checking for rank deficiency.
1253  STS::one(), R.values(), R.stride(),
1254  X.values(), X.stride());
1255 
1256  // Check for Infs or NaNs in X. If there are any, then
1257  // assume that TRSM failed, and use a more robust algorithm.
1258  if (ops.infNaNCount (X, false) != 0) {
1259 
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;
1263 
1264  // Restore X from the copy.
1265  X.assign (B);
1266 
1267  // Find the minimum-norm solution to the least-squares
1268  // problem $\min_x \|RX - B\|_2$, using the singular value
1269  // decomposition (SVD).
1271  if (side == LEFT_SIDE) {
1272  // _GELSS overwrites its matrix input, so make a copy.
1273  mat_type R_copy (Teuchos::Copy, R_view, N, N);
1274 
1275  // Zero out the lower triangle of R_copy, since the
1276  // mat_type constructor copies all the entries, not just
1277  // the upper triangle. _GELSS will read all the entries
1278  // of the input matrix.
1279  ops.zeroOutStrictLowerTriangle (R_copy);
1280 
1281  // Solve the least-squares problem.
1282  rank = solveLeastSquaresUsingSVD (R_copy, X);
1283  } else {
1284  // If solving with R on the right-hand side, the interface
1285  // requires that instead of solving $\min \|XR - B\|_2$,
1286  // we have to solve $\min \|R^* X^* - B^*\|_2$. We
1287  // compute (conjugate) transposes in newly allocated
1288  // temporary matrices X_star resp. R_star. (B is already
1289  // in X and _GELSS overwrites its input vector X with the
1290  // solution.)
1291  mat_type X_star (X.numCols(), X.numRows());
1292  ops.conjugateTranspose (X_star, X);
1293  mat_type R_star (N, N); // Filled with zeros automatically.
1294  ops.conjugateTransposeOfUpperTriangular (R_star, R);
1295 
1296  // Solve the least-squares problem.
1297  rank = solveLeastSquaresUsingSVD (R_star, X_star);
1298 
1299  // Copy the transpose of X_star back into X.
1300  ops.conjugateTranspose (X, X_star);
1301  }
1302  if (rank < N) {
1303  foundRankDeficiency = true;
1304  }
1305  }
1306  } else {
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.");
1311  }
1312  return std::make_pair (rank, foundRankDeficiency);
1313  }
1314 
1315 
1316  public:
1325  bool
1326  testGivensRotations (std::ostream& out)
1327  {
1328  using std::endl;
1329 
1330  out << "Testing Givens rotations:" << endl;
1331  Scalar x = STS::random();
1332  Scalar y = STS::random();
1333  out << " x = " << x << ", y = " << y << endl;
1334 
1335  Scalar theCosine, theSine, result;
1336  blas_type blas;
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;
1342 
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;
1347 
1348  // Allow only a tiny bit of wiggle room for zeroing-out of y.
1349  if (STS::magnitude(y) > 2*STS::eps())
1350  return false;
1351  else
1352  return true;
1353  }
1354 
1384  bool
1385  testUpdateColumn (std::ostream& out,
1386  const int numCols,
1387  const bool testBlockGivens=false,
1388  const bool extraVerbose=false)
1389  {
1390  using Teuchos::Array;
1391  using std::endl;
1392 
1393  TEUCHOS_TEST_FOR_EXCEPTION(numCols <= 0, std::invalid_argument,
1394  "numCols = " << numCols << " <= 0.");
1395  const int numRows = numCols + 1;
1396 
1397  mat_type H (numRows, numCols);
1398  mat_type z (numRows, 1);
1399 
1400  mat_type R_givens (numRows, numCols);
1401  mat_type y_givens (numRows, 1);
1402  mat_type z_givens (numRows, 1);
1403  Array<Scalar> theCosines (numCols);
1404  Array<Scalar> theSines (numCols);
1405 
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);
1412 
1413  mat_type R_lapack (numRows, numCols);
1414  mat_type y_lapack (numRows, 1);
1415  mat_type z_lapack (numRows, 1);
1416 
1417  // Make a random least-squares problem.
1418  makeRandomProblem (H, z);
1419  if (extraVerbose) {
1420  printMatrix<Scalar> (out, "H", H);
1421  printMatrix<Scalar> (out, "z", z);
1422  }
1423 
1424  // Set up the right-hand side copies for each of the methods.
1425  // Each method is free to overwrite its given right-hand side.
1426  z_givens.assign (z);
1427  if (testBlockGivens) {
1428  z_blockGivens.assign (z);
1429  }
1430  z_lapack.assign (z);
1431 
1432  //
1433  // Imitate how one would update the least-squares problem in a
1434  // typical GMRES implementation, for each updating method.
1435  //
1436  // Update using Givens rotations, one at a time.
1437  magnitude_type residualNormGivens = STM::zero();
1438  for (int curCol = 0; curCol < numCols; ++curCol) {
1439  residualNormGivens = updateColumnGivens (H, R_givens, y_givens, z_givens,
1440  theCosines, theSines, curCol);
1441  }
1442  solveGivens (y_givens, R_givens, z_givens, numCols-1);
1443 
1444  // Update using the "panel left-looking" Givens approach, with
1445  // the given panel width.
1446  magnitude_type residualNormBlockGivens = STM::zero();
1447  if (testBlockGivens) {
1448  const bool testBlocksAtATime = true;
1449  if (testBlocksAtATime) {
1450  // Blocks of columns at a time.
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);
1456  }
1457  } else {
1458  // One column at a time. This is good as a sanity check
1459  // to make sure updateColumnsGivens() with a single column
1460  // does the same thing as updateColumnGivens().
1461  for (int startCol = 0; startCol < numCols; ++startCol) {
1462  residualNormBlockGivens =
1463  updateColumnsGivens (H, R_blockGivens, y_blockGivens, z_blockGivens,
1464  blockCosines, blockSines, startCol, startCol);
1465  }
1466  }
1467  // The panel version of Givens should compute the same
1468  // cosines and sines as the non-panel version, and should
1469  // update the right-hand side z in the same way. Thus, we
1470  // should be able to use the same triangular solver.
1471  solveGivens (y_blockGivens, R_blockGivens, z_blockGivens, numCols-1);
1472  }
1473 
1474  // Solve using LAPACK's least-squares solver.
1475  const magnitude_type residualNormLapack =
1476  solveLapack (H, R_lapack, y_lapack, z_lapack, numCols-1);
1477 
1478  // Compute the condition number of the least-squares problem.
1479  // This requires a residual, so use the residual from the
1480  // LAPACK method. All that the method needs for an accurate
1481  // residual norm is forward stability.
1482  const magnitude_type leastSquaresCondNum =
1483  leastSquaresConditionNumber (H, z, residualNormLapack);
1484 
1485  // Compute the relative least-squares solution error for both
1486  // Givens methods. We assume that the LAPACK solution is
1487  // "exact" and compare against the Givens rotations solution.
1488  // This is taking liberties with the definition of condition
1489  // number, but it's the best we can do, since we don't know
1490  // the exact solution and don't have an extended-precision
1491  // solver.
1492 
1493  // The solution lives only in y[0 .. 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);
1497 
1498  const magnitude_type givensSolutionError =
1499  solutionError (y_givens_view, y_lapack_view);
1500  const magnitude_type blockGivensSolutionError = testBlockGivens ?
1501  solutionError (y_blockGivens_view, y_lapack_view) :
1502  STM::zero();
1503 
1504  // If printing out the matrices, copy out the upper triangular
1505  // factors for printing. (Both methods are free to leave data
1506  // below the lower triangle.)
1507  if (extraVerbose) {
1508  mat_type R_factorFromGivens (numCols, numCols);
1509  mat_type R_factorFromBlockGivens (numCols, numCols);
1510  mat_type R_factorFromLapack (numCols, numCols);
1511 
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);
1517  }
1518  R_factorFromLapack(i,j) = R_lapack(i,j);
1519  }
1520  }
1521 
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);
1525 
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);
1530  }
1531 
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);
1535  }
1536 
1537  // Compute the (Frobenius) norm of the original matrix H.
1538  const magnitude_type H_norm = H.normFrobenius();
1539 
1540  out << "||H||_F = " << H_norm << endl;
1541 
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;
1547  }
1548  out << "||H y_lapack - z||_2 / ||H||_F = "
1549  << leastSquaresResidualNorm (H, y_lapack_view, z) / H_norm << endl;
1550 
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;
1556  }
1557 
1558  out << "Least-squares condition number = "
1559  << leastSquaresCondNum << endl;
1560 
1561  // Now for the controversial part of the test: judging whether
1562  // we succeeded. This includes the problem's condition
1563  // number, which is a measure of the maximum perturbation in
1564  // the solution for which we can still say that the solution
1565  // is valid. We include a little wiggle room by including a
1566  // factor proportional to the square root of the number of
1567  // floating-point operations that influence the last entry
1568  // (the conventional Wilkinsonian heuristic), times 10 for
1569  // good measure.
1570  //
1571  // (The square root looks like it has something to do with an
1572  // average-case probabilistic argument, but doesn't really.
1573  // What's an "average problem"?)
1574  const magnitude_type wiggleFactor =
1575  10 * STM::squareroot( numRows*numCols );
1576  const magnitude_type solutionErrorBoundFactor =
1577  wiggleFactor * leastSquaresCondNum;
1578  const magnitude_type solutionErrorBound =
1579  solutionErrorBoundFactor * STS::eps();
1580  out << "Solution error bound: " << solutionErrorBoundFactor
1581  << " * eps = " << solutionErrorBound << endl;
1582 
1583  // Remember that NaN is not greater than, not less than, and
1584  // not equal to any other number, including itself. Some
1585  // compilers will rudely optimize away the "x != x" test.
1586  if (STM::isnaninf (solutionErrorBound)) {
1587  // Hm, the solution error bound is Inf or NaN. This
1588  // probably means that the test problem was generated
1589  // incorrectly. We should return false in this case.
1590  return false;
1591  } else { // solution error bound is finite.
1592  if (STM::isnaninf (givensSolutionError)) {
1593  return false;
1594  } else if (givensSolutionError > solutionErrorBound) {
1595  return false;
1596  } else if (testBlockGivens) {
1597  if (STM::isnaninf (blockGivensSolutionError)) {
1598  return false;
1599  } else if (blockGivensSolutionError > solutionErrorBound) {
1600  return false;
1601  } else { // Givens and Block Givens tests succeeded.
1602  return true;
1603  }
1604  } else { // Not testing block Givens; Givens test succeeded.
1605  return true;
1606  }
1607  }
1608  }
1609 
1625  bool
1626  testTriangularSolves (std::ostream& out,
1627  const int testProblemSize,
1628  const ERobustness robustness,
1629  const bool verbose=false)
1630  {
1631  using Teuchos::LEFT_SIDE;
1632  using Teuchos::RIGHT_SIDE;
1633  using std::endl;
1635 
1636  Teuchos::oblackholestream blackHole;
1637  std::ostream& verboseOut = verbose ? out : blackHole;
1638 
1639  verboseOut << "Testing upper triangular solves" << endl;
1640  //
1641  // Construct an upper triangular linear system to solve.
1642  //
1643  verboseOut << "-- Generating test matrix" << endl;
1644  const int N = testProblemSize;
1645  mat_type R (N, N);
1646  // Fill the upper triangle of R with random numbers.
1647  for (int j = 0; j < N; ++j) {
1648  for (int i = 0; i <= j; ++i) {
1649  R(i,j) = STS::random ();
1650  }
1651  }
1652  // Compute the Frobenius norm of R for later use.
1653  const magnitude_type R_norm = R.normFrobenius ();
1654  // Fill the right-hand side B with random numbers.
1655  mat_type B (N, 1);
1656  B.random ();
1657  // Compute the Frobenius norm of B for later use.
1658  const magnitude_type B_norm = B.normFrobenius ();
1659 
1660  // Save a copy of the original upper triangular system.
1661  mat_type R_copy (Teuchos::Copy, R, N, N);
1662  mat_type B_copy (Teuchos::Copy, B, N, 1);
1663 
1664  // Solution vector.
1665  mat_type X (N, 1);
1666 
1667  // Solve RX = B.
1668  verboseOut << "-- Solving RX=B" << endl;
1669  // We're ignoring the return values for now.
1670  (void) solveUpperTriangularSystem (LEFT_SIDE, X, R, B, robustness);
1671  // Test the residual error.
1672  mat_type Resid (N, 1);
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)
1679  << endl;
1680 
1681  // Restore R and B.
1682  R.assign (R_copy);
1683  B.assign (B_copy);
1684 
1685  //
1686  // Set up a right-side test problem: YR = B^*.
1687  //
1688  mat_type Y (1, N);
1689  mat_type B_star (1, N);
1690  ops.conjugateTranspose (B_star, B);
1691  mat_type B_star_copy (1, N);
1692  B_star_copy.assign (B_star);
1693  // Solve YR = B^*.
1694  verboseOut << "-- Solving YR=B^*" << endl;
1695  // We're ignoring the return values for now.
1696  (void) solveUpperTriangularSystem (RIGHT_SIDE, Y, R, B_star, robustness);
1697  // Test the residual error.
1698  mat_type Resid2 (1, N);
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)
1704  << endl;
1705 
1706  // FIXME (mfh 14 Oct 2011) The test always "passes" for now;
1707  // you have to inspect the output in order to see whether it
1708  // succeeded. We really should fix the above to use the
1709  // infinity-norm bounds in Higham's book for triangular
1710  // solves. That would automate the test.
1711  return true;
1712  }
1713 
1714  private:
1716  std::ostream& warn_;
1717 
1722  ERobustness defaultRobustness_;
1723 
1724  private:
1736  int
1737  solveLeastSquaresUsingSVD (mat_type& A, mat_type& X)
1738  {
1739  using Teuchos::Array;
1741 
1742  if (defaultRobustness_ > ROBUSTNESS_SOME) {
1743  int count = ops.infNaNCount (A);
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") << ".");
1748  count = ops.infNaNCount (X);
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") << ".");
1753  }
1754  const int N = std::min (A.numRows(), A.numCols());
1755  lapack_type lapack;
1756 
1757  // Rank of A; to be computed by _GELSS and returned.
1758  int rank = N;
1759 
1760  // Use Scalar's machine precision for the rank tolerance,
1761  // not magnitude_type's machine precision.
1762  const magnitude_type rankTolerance = STS::eps();
1763 
1764  // Array of singular values.
1765  Array<magnitude_type> singularValues (N);
1766 
1767  // Extra workspace. This is only used by _GELSS if Scalar is
1768  // complex. Teuchos::LAPACK presents a unified interface to
1769  // _GELSS that always includes the RWORK argument, even though
1770  // LAPACK's SGELSS and DGELSS don't have the RWORK argument.
1771  // We always allocate at least one entry so that &rwork[0]
1772  // makes sense.
1773  Array<magnitude_type> rwork (1);
1774  if (STS::isComplex) {
1775  rwork.resize (std::max (1, 5 * N));
1776  }
1777  //
1778  // Workspace query
1779  //
1780  Scalar lworkScalar = STS::one(); // To be set by workspace query
1781  int info = 0;
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.");
1793  // Allocate workspace. Size > 0 means &work[0] makes sense.
1794  Array<Scalar> work (std::max (1, lwork));
1795  // Solve the least-squares problem.
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.");
1802  return rank;
1803  }
1804 
1811  void
1812  solveGivens (mat_type& y,
1813  mat_type& R,
1814  const mat_type& z,
1815  const int curCol)
1816  {
1817  const int numRows = curCol + 2;
1818 
1819  // Now that we have the updated R factor of H, and the updated
1820  // right-hand side z, solve the least-squares problem by
1821  // solving the upper triangular linear system Ry=z for y.
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());
1825 
1827  R_view, z_view, defaultRobustness_);
1828  }
1829 
1831  void
1832  makeRandomProblem (mat_type& H, mat_type& z)
1833  {
1834  // In GMRES, z always starts out with only the first entry
1835  // being nonzero. That entry always has nonnegative real part
1836  // and zero imaginary part, since it is the initial residual
1837  // norm.
1838  H.random ();
1839  // Zero out the entries below the subdiagonal of H, so that it
1840  // is upper Hessenberg.
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();
1844  }
1845  }
1846  // Initialize z, the right-hand side of the least-squares
1847  // problem. Make the first entry of z nonzero.
1848  {
1849  // It's still possible that a random number will come up
1850  // zero after 1000 trials, but unlikely. Nevertheless, it's
1851  // still important not to allow an infinite loop, for
1852  // example if the pseudorandom number generator is broken
1853  // and always returns zero.
1854  const int numTrials = 1000;
1855  magnitude_type z_init = STM::zero();
1856  for (int trial = 0; trial < numTrials && z_init == STM::zero(); ++trial) {
1857  z_init = STM::random();
1858  }
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.");
1865  const magnitude_type z_first = (z_init < 0) ? -z_init : z_init;
1866 
1867  // NOTE I'm assuming here that "scalar_type = magnitude_type"
1868  // assignments make sense.
1869  z(0,0) = z_first;
1870  }
1871  }
1872 
1876  void
1877  computeGivensRotation (const Scalar& x,
1878  const Scalar& y,
1879  Scalar& theCosine,
1880  Scalar& theSine,
1881  Scalar& result)
1882  {
1883  // _LARTG, an LAPACK aux routine, is slower but more accurate
1884  // than the BLAS' _ROTG.
1885  const bool useLartg = false;
1886 
1887  if (useLartg) {
1888  lapack_type lapack;
1889  // _LARTG doesn't clobber its input arguments x and y.
1890  lapack.LARTG (x, y, &theCosine, &theSine, &result);
1891  } else {
1892  // _ROTG clobbers its first two arguments. x is overwritten
1893  // with the result of applying the Givens rotation: [x; y] ->
1894  // [x (on output); 0]. y is overwritten with the "fast"
1895  // Givens transform (see Golub and Van Loan, 3rd ed.).
1896  Scalar x_temp = x;
1897  Scalar y_temp = y;
1898  blas_type blas;
1899  blas.ROTG (&x_temp, &y_temp, &theCosine, &theSine);
1900  result = x_temp;
1901  }
1902  }
1903 
1905  void
1906  singularValues (const mat_type& A,
1908  {
1909  using Teuchos::Array;
1910  using Teuchos::ArrayView;
1911 
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 "
1920  "matrix A.");
1921 
1922  // Compute the condition number of the matrix A, using a singular
1923  // value decomposition (SVD). LAPACK's SVD routine overwrites the
1924  // input matrix, so make a copy.
1925  mat_type A_copy (numRows, numCols);
1926  A_copy.assign (A);
1927 
1928  // Workspace query.
1929  lapack_type lapack;
1930  int info = 0;
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);
1937 
1938  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1939  "LAPACK _GESVD workspace query failed with INFO = "
1940  << 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.");
1945  // Make sure that the workspace array always has positive
1946  // length, so that &work[0] makes sense.
1947  Teuchos::Array<Scalar> work (std::max (1, lwork));
1948 
1949  // Compute the singular values of A.
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 << ".");
1956  }
1957 
1965  std::pair<magnitude_type, magnitude_type>
1966  extremeSingularValues (const mat_type& A)
1967  {
1968  using Teuchos::Array;
1969 
1970  const int numRows = A.numRows();
1971  const int numCols = A.numCols();
1972 
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]);
1976  }
1977 
1988  magnitude_type
1989  leastSquaresConditionNumber (const mat_type& A,
1990  const mat_type& b,
1991  const magnitude_type& residualNorm)
1992  {
1993  // Extreme singular values of A.
1994  const std::pair<magnitude_type, magnitude_type> sigmaMaxMin =
1995  extremeSingularValues (A);
1996 
1997  // Our solvers currently assume that H has full rank. If the
1998  // test matrix doesn't have full rank, we stop right away.
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 "
2002  "zero.");
2003  // 2-norm condition number of A. We checked above that the
2004  // denominator is nonzero.
2005  const magnitude_type A_cond = sigmaMaxMin.first / sigmaMaxMin.second;
2006 
2007  // "Theta" in the variable names below refers to the angle between
2008  // the vectors b and A*x, where x is the computed solution. It
2009  // measures whether the residual norm is large (near ||b||) or
2010  // small (near 0).
2011  const magnitude_type sinTheta = residualNorm / b.normFrobenius();
2012 
2013  // \sin^2 \theta + \cos^2 \theta = 1.
2014  //
2015  // The range of sine is [-1,1], so squaring it won't overflow.
2016  // We still have to check whether sinTheta > 1, though. This
2017  // is impossible in exact arithmetic, assuming that the
2018  // least-squares solver worked (b-A*0 = b and x minimizes
2019  // ||b-A*x||_2, so ||b-A*0||_2 >= ||b-A*x||_2). However, it
2020  // might just be possible in floating-point arithmetic. We're
2021  // just looking for an estimate, so if sinTheta > 1, we cap it
2022  // at 1.
2023  const magnitude_type cosTheta = (sinTheta > STM::one()) ?
2024  STM::zero() : STM::squareroot (1 - sinTheta * sinTheta);
2025 
2026  // This may result in Inf, if cosTheta is zero. That's OK; in
2027  // that case, the condition number of the (full-rank)
2028  // least-squares problem is rightfully infinite.
2029  const magnitude_type tanTheta = sinTheta / cosTheta;
2030 
2031  // Condition number for the full-rank least-squares problem.
2032  return 2 * A_cond / cosTheta + tanTheta * A_cond * A_cond;
2033  }
2034 
2036  magnitude_type
2037  leastSquaresResidualNorm (const mat_type& A,
2038  const mat_type& x,
2039  const mat_type& b)
2040  {
2041  mat_type r (b.numRows(), b.numCols());
2042 
2043  // r := b - A*x
2044  r.assign (b);
2045  LocalDenseMatrixOps<Scalar> ops;
2046  ops.matMatMult (STS::one(), r, -STS::one(), A, x);
2047  return r.normFrobenius ();
2048  }
2049 
2054  magnitude_type
2055  solutionError (const mat_type& x_approx,
2056  const mat_type& x_exact)
2057  {
2058  const int numRows = x_exact.numRows();
2059  const int numCols = x_exact.numCols();
2060 
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);
2065  }
2066  }
2067  const magnitude_type scalingFactor = x_exact.normFrobenius();
2068 
2069  // If x_exact has zero norm, just use the absolute difference.
2070  return x_diff.normFrobenius() /
2071  (scalingFactor == STM::zero() ? STM::one() : scalingFactor);
2072  }
2073 
2120  magnitude_type
2121  updateColumnGivens (const mat_type& H,
2122  mat_type& R,
2123  mat_type& y,
2124  mat_type& z,
2127  const int curCol)
2128  {
2129  using std::cerr;
2130  using std::endl;
2131 
2132  const int numRows = curCol + 2; // curCol is zero-based
2133  const int LDR = R.stride();
2134  const bool extraDebug = false;
2135 
2136  if (extraDebug) {
2137  cerr << "updateColumnGivens: curCol = " << curCol << endl;
2138  }
2139 
2140  // View of H( 1:curCol+1, curCol ) (in Matlab notation, if
2141  // curCol were a one-based index, as it would be in Matlab but
2142  // is not here).
2143  const mat_type H_col (Teuchos::View, H, numRows, 1, 0, curCol);
2144 
2145  // View of R( 1:curCol+1, curCol ) (again, in Matlab notation,
2146  // if curCol were a one-based index).
2147  mat_type R_col (Teuchos::View, R, numRows, 1, 0, curCol);
2148 
2149  // 1. Copy the current column from H into R, where it will be
2150  // modified.
2151  R_col.assign (H_col);
2152 
2153  if (extraDebug) {
2154  printMatrix<Scalar> (cerr, "R_col before ", R_col);
2155  }
2156 
2157  // 2. Apply all the previous Givens rotations, if any, to the
2158  // current column of the matrix.
2159  blas_type blas;
2160  for (int j = 0; j < curCol; ++j) {
2161  // BLAS::ROT really should take "const Scalar*" for these
2162  // arguments, but it wants a "Scalar*" instead, alas.
2163  Scalar theCosine = theCosines[j];
2164  Scalar theSine = theSines[j];
2165 
2166  if (extraDebug) {
2167  cerr << " j = " << j << ": (cos,sin) = "
2168  << theCosines[j] << "," << theSines[j] << endl;
2169  }
2170  blas.ROT (1, &R_col(j,0), LDR, &R_col(j+1,0), LDR,
2171  &theCosine, &theSine);
2172  }
2173  if (extraDebug && curCol > 0) {
2174  printMatrix<Scalar> (cerr, "R_col after applying previous "
2175  "Givens rotations", R_col);
2176  }
2177 
2178  // 3. Calculate new Givens rotation for R(curCol, curCol),
2179  // R(curCol+1, curCol).
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;
2185 
2186  if (extraDebug) {
2187  cerr << " New cos,sin = " << theCosine << "," << theSine << endl;
2188  }
2189 
2190  // 4. _Apply_ the new Givens rotation. We don't need to
2191  // invoke _ROT here, because computeGivensRotation()
2192  // already gives us the result: [x; y] -> [result; 0].
2193  R_col(curCol, 0) = result;
2194  R_col(curCol+1, 0) = STS::zero();
2195 
2196  if (extraDebug) {
2197  printMatrix<Scalar> (cerr, "R_col after applying current "
2198  "Givens rotation", R_col);
2199  }
2200 
2201  // 5. Apply the resulting Givens rotation to z (the right-hand
2202  // side of the projected least-squares problem).
2203  //
2204  // We prefer overgeneralization to undergeneralization by assuming
2205  // here that z may have more than one column.
2206  const int LDZ = z.stride();
2207  blas.ROT (z.numCols(),
2208  &z(curCol,0), LDZ, &z(curCol+1,0), LDZ,
2209  &theCosine, &theSine);
2210 
2211  if (extraDebug) {
2212  //mat_type R_after (Teuchos::View, R, numRows, numRows-1);
2213  //printMatrix<Scalar> (cerr, "R_after", R_after);
2214  //mat_type z_after (Teuchos::View, z, numRows, z.numCols());
2215  printMatrix<Scalar> (cerr, "z_after", z);
2216  }
2217 
2218  // The last entry of z is the nonzero part of the residual of the
2219  // least-squares problem. Its magnitude gives the residual 2-norm
2220  // of the least-squares problem.
2221  return STS::magnitude( z(numRows-1, 0) );
2222  }
2223 
2257  magnitude_type
2258  solveLapack (const mat_type& H,
2259  mat_type& R,
2260  mat_type& y,
2261  mat_type& z,
2262  const int curCol)
2263  {
2264  const int numRows = curCol + 2;
2265  const int numCols = curCol + 1;
2266  const int LDR = R.stride();
2267 
2268  // Copy H( 1:curCol+1, 1:curCol ) into R( 1:curCol+1, 1:curCol ).
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);
2272 
2273  // The LAPACK least-squares solver overwrites the right-hand side
2274  // vector with the solution, so first copy z into y.
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);
2278 
2279  // Workspace query for the least-squares routine.
2280  int info = 0;
2281  Scalar lworkScalar = STS::zero();
2282  lapack_type lapack;
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" : "") << ".");
2293  std::logic_error,
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.");
2299  std::logic_error,
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.");
2304  // Cast workspace from Scalar to int. Scalar may be complex,
2305  // hence the request for the real part. Don't ask for the
2306  // magnitude, since computing the magnitude may overflow due
2307  // to squaring and square root to int. Hopefully LAPACK
2308  // doesn't ever overflow int this way.
2309  const int lwork = std::max (1, static_cast<int> (STS::real (lworkScalar)));
2310 
2311  // Allocate workspace for solving the least-squares problem.
2312  Teuchos::Array<Scalar> work (lwork);
2313 
2314  // Solve the least-squares problem. The ?: operator prevents
2315  // accessing the first element of the work array, if it has
2316  // length zero.
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),
2321  lwork, &info);
2322 
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" : "") << ".");
2329  // Extract the projected least-squares problem's residual error.
2330  // It's the magnitude of the last entry of y_view on output from
2331  // LAPACK's least-squares solver.
2332  return STS::magnitude( y_view(numRows-1, 0) );
2333  }
2334 
2345  magnitude_type
2346  updateColumnsGivens (const mat_type& H,
2347  mat_type& R,
2348  mat_type& y,
2349  mat_type& z,
2352  const int startCol,
2353  const int endCol)
2354  {
2355  TEUCHOS_TEST_FOR_EXCEPTION(startCol > endCol, std::invalid_argument,
2356  "updateColumnGivens: startCol = " << startCol
2357  << " > endCol = " << endCol << ".");
2358  magnitude_type lastResult = STM::zero();
2359  // [startCol, endCol] is an inclusive range.
2360  for (int curCol = startCol; curCol <= endCol; ++curCol) {
2361  lastResult = updateColumnGivens (H, R, y, z, theCosines, theSines, curCol);
2362  }
2363  return lastResult;
2364  }
2365 
2378  magnitude_type
2379  updateColumnsGivensBlock (const mat_type& H,
2380  mat_type& R,
2381  mat_type& y,
2382  mat_type& z,
2385  const int startCol,
2386  const int endCol)
2387  {
2388  const int numRows = endCol + 2;
2389  const int numColsToUpdate = endCol - startCol + 1;
2390  const int LDR = R.stride();
2391 
2392  // 1. Copy columns [startCol, endCol] from H into R, where they
2393  // will be modified.
2394  {
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);
2398  }
2399 
2400  // 2. Apply all the previous Givens rotations, if any, to
2401  // columns [startCol, endCol] of the matrix. (Remember
2402  // that we're using a left-looking QR factorization
2403  // approach; we haven't yet touched those columns.)
2404  blas_type blas;
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]);
2409  }
2410 
2411  // 3. Update each column in turn of columns [startCol, endCol].
2412  for (int curCol = startCol; curCol < endCol; ++curCol) {
2413  // a. Apply the Givens rotations computed in previous
2414  // iterations of this loop to the current column of R.
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]);
2418  }
2419  // b. Calculate new Givens rotation for R(curCol, curCol),
2420  // R(curCol+1, curCol).
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;
2426 
2427  // c. _Apply_ the new Givens rotation. We don't need to
2428  // invoke _ROT here, because computeGivensRotation()
2429  // already gives us the result: [x; y] -> [result; 0].
2430  R(curCol+1, curCol) = result;
2431  R(curCol+1, curCol) = STS::zero();
2432 
2433  // d. Apply the resulting Givens rotation to z (the right-hand
2434  // side of the projected least-squares problem).
2435  //
2436  // We prefer overgeneralization to undergeneralization by
2437  // assuming here that z may have more than one column.
2438  const int LDZ = z.stride();
2439  blas.ROT (z.numCols(),
2440  &z(curCol,0), LDZ, &z(curCol+1,0), LDZ,
2441  &theCosine, &theSine);
2442  }
2443 
2444  // The last entry of z is the nonzero part of the residual of the
2445  // least-squares problem. Its magnitude gives the residual 2-norm
2446  // of the least-squares problem.
2447  return STS::magnitude( z(numRows-1, 0) );
2448  }
2449  }; // class ProjectedLeastSquaresSolver
2450  } // namespace details
2451 } // namespace Belos
2452 
2453 #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 Thu Mar 28 2024 09:24:29 for Belos by doxygen 1.8.5