Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
52#include "Teuchos_SerialDenseMatrix.hpp"
53#include "Teuchos_StandardParameterEntryValidators.hpp"
54
58
59namespace 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,
79 const Teuchos::SerialDenseMatrix<int, Scalar>& A)
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,
115 const Teuchos::SerialDenseMatrix<int, Scalar>& A,
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;
166 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
167
175 Teuchos::SerialDenseMatrix<int,Scalar> H;
176
188 Teuchos::SerialDenseMatrix<int,Scalar> R;
189
200 Teuchos::SerialDenseMatrix<int,Scalar> y;
201
210 Teuchos::SerialDenseMatrix<int,Scalar> z;
211
221 Teuchos::Array<Scalar> theCosines;
222
227 Teuchos::Array<Scalar> theSines;
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
265 reset (const typename Teuchos::ScalarTraits<Scalar>::magnitudeType beta)
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
292 reallocateAndReset (const typename Teuchos::ScalarTraits<Scalar>::magnitudeType beta,
293 const int maxNumIterations)
294 {
295 typedef Teuchos::ScalarTraits<Scalar> STS;
296 typedef Teuchos::ScalarTraits<magnitude_type> STM;
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;
356 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
359 typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
360
361 private:
362 typedef Teuchos::ScalarTraits<scalar_type> STS;
363 typedef Teuchos::ScalarTraits<magnitude_type> STM;
364 typedef Teuchos::BLAS<int, scalar_type> blas_type;
365 typedef Teuchos::LAPACK<int, scalar_type> lapack_type;
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
414 partition (Teuchos::RCP<mat_type>& A_11,
415 Teuchos::RCP<mat_type>& A_21,
416 Teuchos::RCP<mat_type>& A_12,
417 Teuchos::RCP<mat_type>& A_22,
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
484 TEUCHOS_TEST_FOR_EXCEPTION(
485 B.numRows() != numRows || B.numCols() != numCols,
486 std::invalid_argument,
487 "matAdd: The input matrices A and B have incompatible dimensions. "
488 "A is " << numRows << " x " << numCols << ", but B is " <<
489 B.numRows () << " x " << B.numCols () << ".");
490 if (numRows == 0 || numCols == 0) {
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
511 TEUCHOS_TEST_FOR_EXCEPTION(
512 B.numRows() != numRows || B.numCols() != numCols,
513 std::invalid_argument,
514 "matSub: The input matrices A and B have incompatible dimensions. "
515 "A is " << numRows << " x " << numCols << ", but B is " <<
516 B.numRows () << " x " << B.numCols () << ".");
517 if (numRows == 0 || numCols == 0) {
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 {
539 TEUCHOS_TEST_FOR_EXCEPTION(B.numCols() != R.numRows(),
540 std::invalid_argument,
541 "rightUpperTriSolve: R and B have incompatible "
542 "dimensions. B has " << B.numCols() << " columns, "
543 "but R has " << R.numRows() << " rows.");
544 blas_type blas;
545 blas.TRSM (Teuchos::RIGHT_SIDE, Teuchos::UPPER_TRI,
546 Teuchos::NO_TRANS, Teuchos::NON_UNIT_DIAG,
547 R.numCols(), B.numCols(),
548 STS::one(), R.values(), R.stride(),
549 B.values(), B.stride());
550 }
551
557 void
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
566 TEUCHOS_TEST_FOR_EXCEPTION(A.numCols() != B.numRows(),
567 std::invalid_argument,
568 "matMatMult: The input matrices A and B have "
569 "incompatible dimensions. A is " << A.numRows()
570 << " x " << A.numCols() << ", but B is "
571 << B.numRows() << " x " << B.numCols() << ".");
572 TEUCHOS_TEST_FOR_EXCEPTION(A.numRows() != C.numRows(),
573 std::invalid_argument,
574 "matMatMult: The input matrix A and the output "
575 "matrix C have incompatible dimensions. A has "
576 << A.numRows() << " rows, but C has " << C.numRows()
577 << " rows.");
578 TEUCHOS_TEST_FOR_EXCEPTION(B.numCols() != C.numCols(),
579 std::invalid_argument,
580 "matMatMult: The input matrix B and the output "
581 "matrix C have incompatible dimensions. B has "
582 << B.numCols() << " columns, but C has "
583 << C.numCols() << " columns.");
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> >
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> >
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 =
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 =
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 =
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
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
874 inline Teuchos::RCP<Teuchos::ParameterEntryValidator>
876 {
877 using Teuchos::stringToIntegralParameterEntryValidator;
878
879 Teuchos::Array<std::string> strs (3);
883 Teuchos::Array<std::string> docs (3);
884 docs[0] = "Use the BLAS' triangular solve. This may result in Inf or "
885 "NaN output if the triangular matrix is rank deficient.";
886 docs[1] = "Robustness somewhere between \"None\" and \"Lots\".";
887 docs[2] = "Solve the triangular system in a least-squares sense, using "
888 "an SVD-based algorithm. This will always succeed, though the "
889 "solution may not make sense for GMRES.";
890 Teuchos::Array<ERobustness> ints (3);
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;
976 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
979 typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
980
981 private:
982 typedef Teuchos::ScalarTraits<scalar_type> STS;
983 typedef Teuchos::ScalarTraits<magnitude_type> STM;
984 typedef Teuchos::BLAS<int, scalar_type> blas_type;
985 typedef Teuchos::LAPACK<int, scalar_type> lapack_type;
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>
1077 solveUpperTriangularSystem (Teuchos::ESide side,
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>
1124 solveUpperTriangularSystem (Teuchos::ESide side,
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 {
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 }
1207 TEUCHOS_TEST_FOR_EXCEPTION(robustness < ROBUSTNESS_NONE ||
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.
1240 blas.TRSM(side, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
1241 Teuchos::NON_UNIT_DIAG, X.numRows(), X.numCols(),
1242 STS::one(), R.values(), R.stride(),
1243 X.values(), X.stride());
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.
1251 blas.TRSM(side, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
1252 Teuchos::NON_UNIT_DIAG, X.numRows(), X.numCols(),
1253 STS::one(), R.values(), R.stride(),
1254 X.values(), X.stride());
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.
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;
1634 typedef Teuchos::SerialDenseMatrix<int, scalar_type> mat_type;
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
1723
1724 private:
1736 int
1738 {
1739 using Teuchos::Array;
1741
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
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
1826 (void) solveUpperTriangularSystem (Teuchos::LEFT_SIDE, y_view,
1827 R_view, z_view, defaultRobustness_);
1828 }
1829
1831 void
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
1907 Teuchos::ArrayView<magnitude_type> sigmas)
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>
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
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 =
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
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);
2046 ops.matMatMult (STS::one(), r, -STS::one(), A, x);
2047 return r.normFrobenius ();
2048 }
2049
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
2122 mat_type& R,
2123 mat_type& y,
2124 mat_type& z,
2125 Teuchos::ArrayView<scalar_type> theCosines,
2126 Teuchos::ArrayView<scalar_type> theSines,
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
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" : "") << ".");
2292 TEUCHOS_TEST_FOR_EXCEPTION(STS::real(lworkScalar) < STM::zero(),
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.");
2298 TEUCHOS_TEST_FOR_EXCEPTION(STS::isComplex && STS::imag(lworkScalar) != STM::zero(),
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
2347 mat_type& R,
2348 mat_type& y,
2349 mat_type& z,
2350 Teuchos::ArrayView<scalar_type> theCosines,
2351 Teuchos::ArrayView<scalar_type> theSines,
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
2380 mat_type& R,
2381 mat_type& y,
2382 mat_type& z,
2383 Teuchos::ArrayView<scalar_type> theCosines,
2384 Teuchos::ArrayView<scalar_type> theSines,
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
Belos header file which uses auto-configuration information to include necessary C++ headers.
Collection of types and exceptions used within the Belos solvers.
Low-level operations on non-distributed dense matrices.
void rightUpperTriSolve(mat_type &B, const mat_type &R) const
In Matlab notation: B = B / R, where R is upper triangular.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
The type of a dense matrix (or vector) of scalar_type.
void conjugateTranspose(mat_type &A_star, const mat_type &A) const
A_star := (conjugate) transpose of A.
void matSub(mat_type &A, const mat_type &B) const
A := A - B.
void zeroOutStrictLowerTriangle(mat_type &A) const
Zero out everything below the diagonal of A.
int infNaNCount(const mat_type &A, const bool upperTriangular=false) const
Return the number of Inf or NaN entries in the matrix A.
std::pair< bool, std::pair< magnitude_type, magnitude_type > > isUpperTriangular(const mat_type &A) const
Is the matrix A upper triangular / trapezoidal?
void ensureUpperHessenberg(const mat_type &A, const char *const matrixName, const magnitude_type relativeTolerance) const
Throw an exception if A is not "approximately" upper Hessenberg.
void conjugateTransposeOfUpperTriangular(mat_type &L, const mat_type &R) const
L := (conjugate) transpose of R (upper triangular).
void matScale(mat_type &A, const scalar_type &alpha) const
A := alpha * A.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
The type of the magnitude of a scalar_type value.
void partition(Teuchos::RCP< mat_type > &A_11, Teuchos::RCP< mat_type > &A_21, Teuchos::RCP< mat_type > &A_12, Teuchos::RCP< mat_type > &A_22, mat_type &A, const int numRows1, const int numRows2, const int numCols1, const int numCols2)
A -> [A_11, A_21, A_12, A_22].
void matAdd(mat_type &A, const mat_type &B) const
A := A + B.
Teuchos::ScalarTraits< magnitude_type > STM
void axpy(mat_type &Y, const scalar_type &alpha, const mat_type &X) const
Y := Y + alpha * X.
void matMatMult(const scalar_type &beta, mat_type &C, const scalar_type &alpha, const mat_type &A, const mat_type &B) const
C := beta*C + alpha*A*B.
Scalar scalar_type
The template parameter of this class.
std::pair< bool, std::pair< magnitude_type, magnitude_type > > isUpperHessenberg(const mat_type &A) const
Is the matrix A upper Hessenberg?
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.
void ensureUpperHessenberg(const mat_type &A, const char *const matrixName) const
Throw an exception if A is not (strictly) upper Hessenberg.
void ensureUpperTriangular(const mat_type &A, const char *const matrixName) const
Throw an exception if A is not upper triangular / trapezoidal.
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.
"Container" for the GMRES projected least-squares problem.
Teuchos::SerialDenseMatrix< int, Scalar > y
Current solution of the projected least-squares problem.
void reallocateAndReset(const typename Teuchos::ScalarTraits< Scalar >::magnitudeType beta, const int maxNumIterations)
(Re)allocate and reset the projected least-squares problem.
Teuchos::SerialDenseMatrix< int, Scalar > R
Upper triangular factor from the QR factorization of H.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
The type of the magnitude of scalar_type values.
ProjectedLeastSquaresProblem(const int maxNumIterations)
Constructor.
Teuchos::SerialDenseMatrix< int, Scalar > H
The upper Hessenberg matrix from GMRES.
void reset(const typename Teuchos::ScalarTraits< Scalar >::magnitudeType beta)
Reset the projected least-squares problem.
Scalar scalar_type
The type of the entries in the projected least-squares problem.
Teuchos::SerialDenseMatrix< int, Scalar > z
Current right-hand side of the projected least-squares problem.
Teuchos::Array< Scalar > theCosines
Array of cosines from the computed Givens rotations.
Teuchos::Array< Scalar > theSines
Array of sines from the computed Givens rotations.
Methods for solving GMRES' projected least-squares problem.
magnitude_type solveLapack(const mat_type &H, mat_type &R, mat_type &y, mat_type &z, const int curCol)
Solve the least-squares problem using LAPACK.
bool testGivensRotations(std::ostream &out)
Test Givens rotations.
void computeGivensRotation(const Scalar &x, const Scalar &y, Scalar &theCosine, Scalar &theSine, Scalar &result)
Compute the Givens rotation corresponding to [x; y].
magnitude_type updateColumnsGivensBlock(const mat_type &H, mat_type &R, mat_type &y, mat_type &z, Teuchos::ArrayView< scalar_type > theCosines, Teuchos::ArrayView< scalar_type > theSines, const int startCol, const int endCol)
Update columns [startCol,endCol] of the projected least-squares problem.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
The type of a dense matrix (or vector) of scalar_type.
bool testUpdateColumn(std::ostream &out, const int numCols, const bool testBlockGivens=false, const bool extraVerbose=false)
Test update and solve using Givens rotations.
magnitude_type leastSquaresResidualNorm(const mat_type &A, const mat_type &x, const mat_type &b)
(Frobenius norm if b has more than one column).
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.
void solve(ProjectedLeastSquaresProblem< Scalar > &problem, const int curCol)
Solve the projected least-squares problem.
void singularValues(const mat_type &A, Teuchos::ArrayView< magnitude_type > sigmas)
Compute the singular values of A. Store them in the sigmas array.
magnitude_type updateColumn(ProjectedLeastSquaresProblem< Scalar > &problem, const int curCol)
Update column curCol of the projected least-squares problem.
magnitude_type updateColumnsGivens(const mat_type &H, mat_type &R, mat_type &y, mat_type &z, Teuchos::ArrayView< scalar_type > theCosines, Teuchos::ArrayView< scalar_type > theSines, const int startCol, const int endCol)
Update columns [startCol,endCol] of the projected least-squares problem.
Scalar scalar_type
The template parameter of this class.
ERobustness defaultRobustness_
Default robustness level, for things like triangular solves.
magnitude_type updateColumns(ProjectedLeastSquaresProblem< Scalar > &problem, const int startCol, const int endCol)
Update columns [startCol,endCol] of the projected least-squares problem.
std::pair< magnitude_type, magnitude_type > extremeSingularValues(const mat_type &A)
The (largest, smallest) singular values of the given matrix.
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).
ProjectedLeastSquaresSolver(std::ostream &warnStream, const ERobustness defaultRobustness=ROBUSTNESS_NONE)
Constructor.
std::pair< int, bool > solveUpperTriangularSystemInPlace(Teuchos::ESide side, mat_type &X, const mat_type &R)
Solve square upper triangular linear system(s) in place.
std::ostream & warn_
Stream to which to output warnings.
bool testTriangularSolves(std::ostream &out, const int testProblemSize, const ERobustness robustness, const bool verbose=false)
Test upper triangular solves.
std::pair< int, bool > solveUpperTriangularSystem(Teuchos::ESide side, mat_type &X, const mat_type &R, const mat_type &B)
Solve the given square upper triangular linear system(s).
magnitude_type updateColumnGivens(const mat_type &H, mat_type &R, mat_type &y, mat_type &z, Teuchos::ArrayView< scalar_type > theCosines, Teuchos::ArrayView< scalar_type > theSines, const int curCol)
Update current column using Givens rotations.
magnitude_type leastSquaresConditionNumber(const mat_type &A, const mat_type &b, const magnitude_type &residualNorm)
Normwise 2-norm condition number of the least-squares problem.
magnitude_type solutionError(const mat_type &x_approx, const mat_type &x_exact)
||x_approx - x_exact||_2 // ||x_exact||_2.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
The type of the magnitude of a scalar_type value.
int solveLeastSquaresUsingSVD(mat_type &A, mat_type &X)
Solve using the SVD.
void solveGivens(mat_type &y, mat_type &R, const mat_type &z, const int curCol)
Solve the projected least-squares problem, assuming Givens rotations updates.
void makeRandomProblem(mat_type &H, mat_type &z)
Make a random projected least-squares problem.
ERobustness
Robustness level of projected least-squares solver operations.
ERobustness robustnessStringToEnum(const std::string &x)
Convert the given robustness string value to an ERobustness enum.
std::string robustnessEnumToString(const ERobustness x)
Convert the given ERobustness enum value to a string.
Teuchos::RCP< Teuchos::ParameterEntryValidator > robustnessValidator()
Make a ParameterList validator for ERobustness.
Namespace containing implementation details of Belos solvers.