BelosProjectedLeastSquaresSolver.hpp

Go to the documentation of this file.

//@HEADER
// ************************************************************************
//
//                 Belos: Block Linear Solvers Package
//                  Copyright 2004 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER

#ifndef __Belos_ProjectedLeastSquaresSolver_hpp
#define __Belos_ProjectedLeastSquaresSolver_hpp

#include "BelosConfigDefs.hpp"
#include "BelosTypes.hpp"
#include "Teuchos_Array.hpp"
#include "Teuchos_BLAS.hpp"
#include "Teuchos_LAPACK.hpp"
#include "Teuchos_oblackholestream.hpp"
#include "Teuchos_ScalarTraits.hpp"
#include "Teuchos_SerialDenseMatrix.hpp"
#include "Teuchos_StandardParameterEntryValidators.hpp"


namespace Belos {

  namespace details {

    // Anonymous namespace restricts contents to file scope.
    namespace {
      template<class Scalar>
      void 
      printMatrix (std::ostream& out, 
       const std::string& name,
       const Teuchos::SerialDenseMatrix<int, Scalar>& A)
      {
  using std::endl;

        const int numRows = A.numRows();
  const int numCols = A.numCols();

  out << name << " = " << endl << '[';
  if (numCols == 1) {
    // Compact form for column vectors; valid Matlab.
    for (int i = 0; i < numRows; ++i) {
      out << A(i,0);
      if (i < numRows-1) {
        out << "; ";
      } 
    }
  } else {
    for (int i = 0; i < numRows; ++i) {
      for (int j = 0; j < numCols; ++j) {
        out << A(i,j);
        if (j < numCols-1) {
    out << ", ";
        } else if (i < numRows-1) {
    out << ';' << endl;
        }
      }
    }
  }
  out << ']' << endl;
      }

      template<class Scalar>
      void 
      print (std::ostream& out, 
       const Teuchos::SerialDenseMatrix<int, Scalar>& A,
       const std::string& linePrefix)
      {
  using std::endl;

        const int numRows = A.numRows();
  const int numCols = A.numCols();

  out << linePrefix << '[';
  if (numCols == 1) {
    // Compact form for column vectors; valid Matlab.
    for (int i = 0; i < numRows; ++i) {
      out << A(i,0);
      if (i < numRows-1) {
        out << "; ";
      } 
    }
  } else {
    for (int i = 0; i < numRows; ++i) {
      for (int j = 0; j < numCols; ++j) {
        if (numRows > 1) {
    out << linePrefix << "  ";
        }
        out << A(i,j);
        if (j < numCols-1) {
    out << ", ";
        } else if (i < numRows-1) {
    out << ';' << endl;
        }
      }
    }
  }
  out << linePrefix << ']' << endl;
      }
    } // namespace (anonymous)

    template<class Scalar>
    class ProjectedLeastSquaresProblem {
    public:
      typedef Scalar scalar_type;
      typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;

      Teuchos::SerialDenseMatrix<int,Scalar> H;

      Teuchos::SerialDenseMatrix<int,Scalar> R;

      Teuchos::SerialDenseMatrix<int,Scalar> y;

      Teuchos::SerialDenseMatrix<int,Scalar> z;

      Teuchos::Array<Scalar> theCosines;

      Teuchos::Array<Scalar> theSines;
      
      ProjectedLeastSquaresProblem (const int maxNumIterations) :
  H (maxNumIterations+1, maxNumIterations),
  R (maxNumIterations+1, maxNumIterations),
  y (maxNumIterations+1, 1),
  z (maxNumIterations+1, 1),
  theCosines (maxNumIterations+1),
  theSines (maxNumIterations+1)
      {}

      void 
      reset (const typename Teuchos::ScalarTraits<Scalar>::magnitudeType beta)
      {
  typedef Teuchos::ScalarTraits<Scalar> STS;

  // Zero out the right-hand side of the least-squares problem.
  z.putScalar (STS::zero());

  // Promote the initial residual norm from a magnitude type to
  // a scalar type, so we can assign it to the first entry of z.
  const Scalar initialResidualNorm (beta);
  z(0,0) = initialResidualNorm;
      }

      void 
      reallocateAndReset (const typename Teuchos::ScalarTraits<Scalar>::magnitudeType beta,
        const int maxNumIterations) 
      {
  typedef Teuchos::ScalarTraits<Scalar> STS;
  typedef Teuchos::ScalarTraits<magnitude_type> STM;

  TEUCHOS_TEST_FOR_EXCEPTION(beta < STM::zero(), std::invalid_argument,
         "ProjectedLeastSquaresProblem::reset: initial "
         "residual beta = " << beta << " < 0.");
  TEUCHOS_TEST_FOR_EXCEPTION(maxNumIterations <= 0, std::invalid_argument,
         "ProjectedLeastSquaresProblem::reset: maximum number "
         "of iterations " << maxNumIterations << " <= 0.");

  if (H.numRows() < maxNumIterations+1 || H.numCols() < maxNumIterations) {
    const int errcode = H.reshape (maxNumIterations+1, maxNumIterations);
    TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error, 
           "Failed to reshape H into a " << (maxNumIterations+1) 
           << " x " << maxNumIterations << " matrix.");
  }
  (void) H.putScalar (STS::zero());

  if (R.numRows() < maxNumIterations+1 || R.numCols() < maxNumIterations) {
    const int errcode = R.reshape (maxNumIterations+1, maxNumIterations);
    TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error, 
           "Failed to reshape R into a " << (maxNumIterations+1) 
           << " x " << maxNumIterations << " matrix.");
  }
  (void) R.putScalar (STS::zero());

  if (y.numRows() < maxNumIterations+1 || y.numCols() < 1) {
    const int errcode = y.reshape (maxNumIterations+1, 1);
    TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error, 
           "Failed to reshape y into a " << (maxNumIterations+1) 
           << " x " << 1 << " matrix.");
  }
  (void) y.putScalar (STS::zero());

  if (z.numRows() < maxNumIterations+1 || z.numCols() < 1) {
    const int errcode = z.reshape (maxNumIterations+1, 1);
    TEUCHOS_TEST_FOR_EXCEPTION(errcode != 0, std::runtime_error, 
           "Failed to reshape z into a " << (maxNumIterations+1) 
           << " x " << 1 << " matrix.");
  }
  reset (beta);
      }

    };


    template<class Scalar>
    class LocalDenseMatrixOps {
    public:
      typedef Scalar scalar_type;
      typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
      typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;

    private:
      typedef Teuchos::ScalarTraits<scalar_type> STS;
      typedef Teuchos::ScalarTraits<magnitude_type> STM;
      typedef Teuchos::BLAS<int, scalar_type> blas_type;
      typedef Teuchos::LAPACK<int, scalar_type> lapack_type;
      
    public:
      void
      conjugateTranspose (mat_type& A_star, const mat_type& A) const
      {
  for (int i = 0; i < A.numRows(); ++i) {
    for (int j = 0; j < A.numCols(); ++j) {
      A_star(j,i) = STS::conjugate (A(i,j));
    }
  }
      }

      void
      conjugateTransposeOfUpperTriangular (mat_type& L, const mat_type& R) const
      {
  const int N = R.numCols();

  for (int j = 0; j < N; ++j) {
    for (int i = 0; i <= j; ++i) {
      L(j,i) = STS::conjugate (R(i,j));
    }
  }
      }

      void
      zeroOutStrictLowerTriangle (mat_type& A) const
      {
  const int N = std::min (A.numRows(), A.numCols());

  for (int j = 0; j < N; ++j) {
    for (int i = j+1; i < A.numRows(); ++i) {
      A(i,j) = STS::zero();
    }
  }
      }

      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)
      {
  using Teuchos::rcp;
  using Teuchos::View;

  A_11 = rcp (new mat_type (View, A, numRows1, numCols1, 0, 0));
  A_21 = rcp (new mat_type (View, A, numRows2, numCols1, numRows1, 0));
  A_12 = rcp (new mat_type (View, A, numRows1, numCols2, 0, numCols1));
  A_22 = rcp (new mat_type (View, A, numRows2, numCols2, numRows1, numCols1));
      }

      void
      matScale (mat_type& A, const scalar_type& alpha) const
      {
  const int LDA = A.stride();
  const int numRows = A.numRows();
  const int numCols = A.numCols();

  if (numRows == 0 || numCols == 0) {
    return;
  } else {
    for (int j = 0; j < numCols; ++j) {
      scalar_type* const A_j = &A(0,j);

      for (int i = 0; i < numRows; ++i) {
        A_j[i] *= alpha;
      }
    }
  }
      }

      void
      axpy (mat_type& Y, 
      const scalar_type& alpha, 
      const mat_type& X) const
      {
  const int numRows = Y.numRows();
  const int numCols = Y.numCols();

  TEUCHOS_TEST_FOR_EXCEPTION(numRows != X.numRows() || numCols != X.numCols(), 
         std::invalid_argument, "Dimensions of X and Y don't "
         "match.  X is " << X.numRows() << " x " << X.numCols() 
         << ", and Y is " << numRows << " x " << numCols << ".");
  for (int j = 0; j < numCols; ++j) {
    for (int i = 0; i < numRows; ++i) {
      Y(i,j) += alpha * X(i,j);
    }
  }
      }

      void 
      matAdd (mat_type& A, const mat_type& B) const
      {
  const int LDA = A.stride();
  const int LDB = B.stride();
  const int numRows = A.numRows();
  const int numCols = A.numCols();

  TEUCHOS_TEST_FOR_EXCEPTION(B.numRows() != numRows || B.numCols() != numCols,
         std::invalid_argument,
         "matAdd: The input matrices A and B have "
         "incompatible dimensions.  A is " << numRows 
         << " x " << numCols << ", but B is " << B.numRows()
         << " x " << B.numCols() << ".");
  if (numRows == 0 || numCols == 0) {
    return;
  } else {
    for (int j = 0; j < numCols; ++j) {
      scalar_type* const A_j = &A(0,j);
      const scalar_type* const B_j = &B(0,j);

      for (int i = 0; i < numRows; ++i) {
        A_j[i] += B_j[i];
      }
    }
  }
      }

      void 
      matSub (mat_type& A, const mat_type& B) const
      {
  const int LDA = A.stride();
  const int LDB = B.stride();
  const int numRows = A.numRows();
  const int numCols = A.numCols();

  TEUCHOS_TEST_FOR_EXCEPTION(B.numRows() != numRows || B.numCols() != numCols,
         std::invalid_argument,
         "matSub: The input matrices A and B have "
         "incompatible dimensions.  A is " << numRows 
         << " x " << numCols << ", but B is " << B.numRows()
         << " x " << B.numCols() << ".");
  if (numRows == 0 || numCols == 0) {
    return;
  } else {
    for (int j = 0; j < numCols; ++j) {
      scalar_type* const A_j = &A(0,j);
      const scalar_type* const B_j = &B(0,j);

      for (int i = 0; i < numRows; ++i) {
        A_j[i] -= B_j[i];
      }
    }
  }
      }

      void
      rightUpperTriSolve (mat_type& B,
        const mat_type& R) const
      {
  TEUCHOS_TEST_FOR_EXCEPTION(B.numCols() != R.numRows(),
         std::invalid_argument,
         "rightUpperTriSolve: R and B have incompatible "
         "dimensions.  B has " << B.numCols() << " columns, "
         "but R has " << R.numRows() << " rows.");
  blas_type blas;
  blas.TRSM (Teuchos::RIGHT_SIDE, Teuchos::UPPER_TRI, 
       Teuchos::NO_TRANS, Teuchos::NON_UNIT_DIAG, 
       R.numCols(), B.numCols(), 
       STS::one(), R.values(), R.stride(), 
       B.values(), B.stride());
      }

      void
      matMatMult (const scalar_type& beta,
      mat_type& C, 
      const scalar_type& alpha,
      const mat_type& A, 
      const mat_type& B) const
      {
  using Teuchos::NO_TRANS;

  TEUCHOS_TEST_FOR_EXCEPTION(A.numCols() != B.numRows(),
         std::invalid_argument,
         "matMatMult: The input matrices A and B have "
         "incompatible dimensions.  A is " << A.numRows() 
         << " x " << A.numCols() << ", but B is " 
         << B.numRows() << " x " << B.numCols() << ".");
  TEUCHOS_TEST_FOR_EXCEPTION(A.numRows() != C.numRows(),
         std::invalid_argument,
         "matMatMult: The input matrix A and the output "
         "matrix C have incompatible dimensions.  A has " 
         << A.numRows() << " rows, but C has " << C.numRows()
         << " rows.");
  TEUCHOS_TEST_FOR_EXCEPTION(B.numCols() != C.numCols(),
         std::invalid_argument,
         "matMatMult: The input matrix B and the output "
         "matrix C have incompatible dimensions.  B has " 
         << B.numCols() << " columns, but C has " 
         << C.numCols() << " columns.");
  blas_type blas;
  blas.GEMM (NO_TRANS, NO_TRANS, C.numRows(), C.numCols(), A.numCols(), 
       alpha, A.values(), A.stride(), B.values(), B.stride(),
       beta, C.values(), C.stride());
      }

      int 
      infNaNCount (const mat_type& A, const bool upperTriangular=false) const
      {
  int count = 0;
  for (int j = 0; j < A.numCols(); ++j) {
    if (upperTriangular) {
      for (int i = 0; i <= j && i < A.numRows(); ++i) {
        if (STS::isnaninf (A(i,j))) {
    ++count;
        }
      }
    } else {
      for (int i = 0; i < A.numRows(); ++i) {
        if (STS::isnaninf (A(i,j))) {
    ++count;
        }
      }
    }
  }
  return count;
      }

      std::pair<bool, std::pair<magnitude_type, magnitude_type> >
      isUpperTriangular (const mat_type& A) const
      {
  magnitude_type lowerTri = STM::zero();
  magnitude_type upperTri = STM::zero();
  int count = 0;

  for (int j = 0; j < A.numCols(); ++j) {
    // Compute the Frobenius norm of the upper triangle /
    // trapezoid of A.  The second clause of the loop upper
    // bound is for matrices with fewer rows than columns.
    for (int i = 0; i <= j && i < A.numRows(); ++i) {
      const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
      upperTri += A_ij_mag * A_ij_mag;
    }
    // Scan the strict lower triangle / trapezoid of A.
    for (int i = j+1; i < A.numRows(); ++i) {
      const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
      lowerTri += A_ij_mag * A_ij_mag;
      if (A_ij_mag != STM::zero()) {
        ++count;
      }
    }
  }
  return std::make_pair (count == 0, std::make_pair (lowerTri, upperTri));
      }


      std::pair<bool, std::pair<magnitude_type, magnitude_type> >
      isUpperHessenberg (const mat_type& A) const
      {
  magnitude_type lower = STM::zero();
  magnitude_type upper = STM::zero();
  int count = 0;

  for (int j = 0; j < A.numCols(); ++j) {
    // Compute the Frobenius norm of the upper Hessenberg part
    // of A.  The second clause of the loop upper bound is for
    // matrices with fewer rows than columns.
    for (int i = 0; i <= j+1 && i < A.numRows(); ++i) {
      const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
      upper += A_ij_mag * A_ij_mag;
    }
    // Scan the strict lower part of A.
    for (int i = j+2; i < A.numRows(); ++i) {
      const magnitude_type A_ij_mag = STS::magnitude (A(i,j));
      lower += A_ij_mag * A_ij_mag;
      if (A_ij_mag != STM::zero()) {
        ++count;
      }
    }
  }
  return std::make_pair (count == 0, std::make_pair (lower, upper));
      }

      void
      ensureUpperTriangular (const mat_type& A,
           const char* const matrixName) const
      {
  std::pair<bool, std::pair<magnitude_type, magnitude_type> > result = 
    isUpperTriangular (A);

  TEUCHOS_TEST_FOR_EXCEPTION(! result.first, std::invalid_argument,
         "The " << A.numRows() << " x " << A.numCols() 
         << " matrix " << matrixName << " is not upper "
         "triangular.  ||tril(A)||_F = " 
         << result.second.first << " and ||A||_F = " 
         << result.second.second << ".");
      }

      void
      ensureUpperHessenberg (const mat_type& A,
           const char* const matrixName) const
      {
  std::pair<bool, std::pair<magnitude_type, magnitude_type> > result = 
    isUpperHessenberg (A);

  TEUCHOS_TEST_FOR_EXCEPTION(! result.first, std::invalid_argument,
         "The " << A.numRows() << " x " << A.numCols() 
         << " matrix " << matrixName << " is not upper "
         "triangular.  ||tril(A(2:end, :))||_F = " 
         << result.second.first << " and ||A||_F = " 
         << result.second.second << ".");
      }

      void
      ensureUpperHessenberg (const mat_type& A,
           const char* const matrixName,
           const magnitude_type relativeTolerance) const
      {
  std::pair<bool, std::pair<magnitude_type, magnitude_type> > result = 
    isUpperHessenberg (A);

  if (result.first) {
    // Mollified relative departure from upper Hessenberg.
    const magnitude_type err = (result.second.second == STM::zero() ? 
              result.second.first : 
              result.second.first / result.second.second);
    TEUCHOS_TEST_FOR_EXCEPTION(err > relativeTolerance, std::invalid_argument,
           "The " << A.numRows() << " x " << A.numCols() 
           << " matrix " << matrixName << " is not upper "
           "triangular.  ||tril(A(2:end, :))||_F "
           << (result.second.second == STM::zero() ? "" : " / ||A||_F")
           << " = " << err << " > " << relativeTolerance << ".");
  }
      }

      void
      ensureMinimumDimensions (const mat_type& A, 
             const char* const matrixName,
             const int minNumRows, 
             const int minNumCols) const
      {
  TEUCHOS_TEST_FOR_EXCEPTION(A.numRows() < minNumRows || A.numCols() < minNumCols,
         std::invalid_argument,
         "The matrix " << matrixName << " is " << A.numRows() 
         << " x " << A.numCols() << ", and therefore does not "
         "satisfy the minimum dimensions " << minNumRows 
         << " x " << minNumCols << ".");
      }

      void
      ensureEqualDimensions (const mat_type& A, 
           const char* const matrixName,
           const int numRows, 
           const int numCols) const
      {
  TEUCHOS_TEST_FOR_EXCEPTION(A.numRows() != numRows || A.numCols() != numCols,
         std::invalid_argument,
         "The matrix " << matrixName << " is supposed to be " 
         << numRows << " x " << numCols << ", but is " 
         << A.numRows() << " x " << A.numCols() << " instead.");
      }

    };

    enum ERobustness {
      ROBUSTNESS_NONE,
      ROBUSTNESS_SOME,
      ROBUSTNESS_LOTS,
      ROBUSTNESS_INVALID
    };

    std::string
    robustnessEnumToString (const ERobustness x) 
    {
      const char* strings[] = {"None", "Some", "Lots"};
      TEUCHOS_TEST_FOR_EXCEPTION(x < ROBUSTNESS_NONE || x >= ROBUSTNESS_INVALID,
       std::invalid_argument,
       "Invalid enum value " << x << ".");
      return std::string (strings[x]);
    }

    ERobustness
    robustnessStringToEnum (const std::string& x)
    {
      const char* strings[] = {"None", "Some", "Lots"};
      for (int r = 0; r < static_cast<int> (ROBUSTNESS_INVALID); ++r) {
  if (x == strings[r]) {
    return static_cast<ERobustness> (r);
  }
      }
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
       "Invalid robustness string " << x << ".");
    }

    Teuchos::RCP<Teuchos::ParameterEntryValidator>
    robustnessValidator ()
    {
      using Teuchos::stringToIntegralParameterEntryValidator;

      Teuchos::Array<std::string> strs (3);
      strs[0] = robustnessEnumToString (ROBUSTNESS_NONE);
      strs[1] = robustnessEnumToString (ROBUSTNESS_SOME);
      strs[2] = robustnessEnumToString (ROBUSTNESS_LOTS);
      Teuchos::Array<std::string> docs (3);
      docs[0] = "Use the BLAS' triangular solve.  This may result in Inf or "
  "NaN output if the triangular matrix is rank deficient.";
      docs[1] = "Robustness somewhere between \"None\" and \"Lots\".";
      docs[2] = "Solve the triangular system in a least-squares sense, using "
  "an SVD-based algorithm.  This will always succeed, though the "
  "solution may not make sense for GMRES.";
      Teuchos::Array<ERobustness> ints (3);
      ints[0] = ROBUSTNESS_NONE;
      ints[1] = ROBUSTNESS_SOME;
      ints[2] = ROBUSTNESS_LOTS;
      const std::string pname ("Robustness of Projected Least-Squares Solve");

      return stringToIntegralParameterEntryValidator<ERobustness> (strs, docs,
                   ints, pname);
    }
    
    template<class Scalar>
    class ProjectedLeastSquaresSolver {
    public:
      typedef Scalar scalar_type;
      typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
      typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;

    private:
      typedef Teuchos::ScalarTraits<scalar_type> STS;
      typedef Teuchos::ScalarTraits<magnitude_type> STM;
      typedef Teuchos::BLAS<int, scalar_type> blas_type;
      typedef Teuchos::LAPACK<int, scalar_type> lapack_type;

    public:
      ProjectedLeastSquaresSolver (std::ostream& warnStream,
           const ERobustness defaultRobustness=ROBUSTNESS_NONE) :
  warn_ (warnStream),
  defaultRobustness_ (defaultRobustness) 
      {}

      magnitude_type
      updateColumn (ProjectedLeastSquaresProblem<Scalar>& problem,
        const int curCol) 
      {
  return updateColumnGivens (problem.H, problem.R, problem.y, problem.z,
           problem.theCosines, problem.theSines, curCol);
      }

      magnitude_type
      updateColumns (ProjectedLeastSquaresProblem<Scalar>& problem,
         const int startCol,
         const int endCol) 
      {
  return updateColumnsGivens (problem.H, problem.R, problem.y, problem.z,
            problem.theCosines, problem.theSines, 
            startCol, endCol);
      }

      void
      solve (ProjectedLeastSquaresProblem<Scalar>& problem,
       const int curCol) 
      {
  solveGivens (problem.y, problem.R, problem.z, curCol);
      }

      std::pair<int, bool>
      solveUpperTriangularSystem (Teuchos::ESide side,
          mat_type& X,
          const mat_type& R,
          const mat_type& B)
      {
  return solveUpperTriangularSystem (side, X, R, B, defaultRobustness_);
      }

      std::pair<int, bool>
      solveUpperTriangularSystem (Teuchos::ESide side,
          mat_type& X,
          const mat_type& R,
          const mat_type& B,
          const ERobustness robustness)
      {
  TEUCHOS_TEST_FOR_EXCEPTION(X.numRows() != B.numRows(), std::invalid_argument,
         "The output X and right-hand side B have different "
         "numbers of rows.  X has " << X.numRows() << " rows"
         ", and B has " << B.numRows() << " rows.");
  // If B has more columns than X, we ignore the remaining
  // columns of B when solving the upper triangular system.  If
  // B has _fewer_ columns than X, we can't solve for all the
  // columns of X, so we throw an exception.
  TEUCHOS_TEST_FOR_EXCEPTION(X.numCols() > B.numCols(), std::invalid_argument,
         "The output X has more columns than the "
         "right-hand side B.  X has " << X.numCols() 
         << " columns and B has " << B.numCols() 
         << " columns.");
  // See above explaining the number of columns in B_view.
  mat_type B_view (Teuchos::View, B, B.numRows(), X.numCols());

  // Both the BLAS' _TRSM and LAPACK's _LATRS overwrite the
  // right-hand side with the solution, so first copy B_view
  // into X.
  X.assign (B_view);

  // Solve the upper triangular system.
  return solveUpperTriangularSystemInPlace (side, X, R, robustness);
      }

      std::pair<int, bool>
      solveUpperTriangularSystemInPlace (Teuchos::ESide side,
           mat_type& X,
           const mat_type& R)
      {
  return solveUpperTriangularSystemInPlace (side, X, R, defaultRobustness_);
      }

      std::pair<int, bool>
      solveUpperTriangularSystemInPlace (Teuchos::ESide side,
           mat_type& X,
           const mat_type& R,
           const ERobustness robustness)
      {
  using Teuchos::Array;
  using Teuchos::Copy;
  using Teuchos::LEFT_SIDE;
  using Teuchos::RIGHT_SIDE;
  LocalDenseMatrixOps<Scalar> ops;

  const int M = R.numRows();
  const int N = R.numCols();
  TEUCHOS_TEST_FOR_EXCEPTION(M < N, std::invalid_argument, 
         "The input matrix R has fewer columns than rows.  "
         "R is " << M << " x " << N << ".");
  // Ignore any additional rows of R by working with a square view.
  mat_type R_view (Teuchos::View, R, N, N);

  if (side == LEFT_SIDE) {
    TEUCHOS_TEST_FOR_EXCEPTION(X.numRows() < N, std::invalid_argument,
           "The input/output matrix X has only "
           << X.numRows() << " rows, but needs at least " 
           << N << " rows to match the matrix for a "
           "left-side solve R \\ X.");
  } else if (side == RIGHT_SIDE) {
    TEUCHOS_TEST_FOR_EXCEPTION(X.numCols() < N, std::invalid_argument,
           "The input/output matrix X has only "
           << X.numCols() << " columns, but needs at least " 
           << N << " columns to match the matrix for a "
           "right-side solve X / R.");
  }
  TEUCHOS_TEST_FOR_EXCEPTION(robustness < ROBUSTNESS_NONE || 
         robustness >= ROBUSTNESS_INVALID, 
         std::invalid_argument,
         "Invalid robustness value " << robustness << ".");

  // In robust mode, scan the matrix and right-hand side(s) for
  // Infs and NaNs.  Only look at the upper triangle of the
  // matrix.
  if (robustness > ROBUSTNESS_NONE) {
    int count = ops.infNaNCount (R_view, true);
    TEUCHOS_TEST_FOR_EXCEPTION(count > 0, std::runtime_error,
           "There " << (count != 1 ? "are" : "is")
           << " " << count << " Inf or NaN entr"
           << (count != 1 ? "ies" : "y") 
           << " in the upper triangle of R.");
    count = ops.infNaNCount (X, false);
    TEUCHOS_TEST_FOR_EXCEPTION(count > 0, std::runtime_error,
           "There " << (count != 1 ? "are" : "is")
           << " " << count << " Inf or NaN entr"
           << (count != 1 ? "ies" : "y") << " in the "
           "right-hand side(s) X.");
  }

  // Pair of values to return from this method.
  int rank = N;
  bool foundRankDeficiency = false;

  // Solve for X.
  blas_type blas;

  if (robustness == ROBUSTNESS_NONE) {
    // Fast triangular solve using the BLAS' _TRSM.  This does
    // no checking for rank deficiency.
    blas.TRSM(side, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
        Teuchos::NON_UNIT_DIAG, X.numRows(), X.numCols(), 
        STS::one(), R.values(), R.stride(), 
        X.values(), X.stride());
  } else if (robustness < ROBUSTNESS_INVALID) {
    // Save a copy of X, since X contains the right-hand side on
    // input.
    mat_type B (Copy, X, X.numRows(), X.numCols());

    // Fast triangular solve using the BLAS' _TRSM.  This does
    // no checking for rank deficiency.
    blas.TRSM(side, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
        Teuchos::NON_UNIT_DIAG, X.numRows(), X.numCols(), 
        STS::one(), R.values(), R.stride(), 
        X.values(), X.stride());

    // Check for Infs or NaNs in X.  If there are any, then
    // assume that TRSM failed, and use a more robust algorithm.
    if (ops.infNaNCount (X, false) != 0) {

      warn_ << "Upper triangular solve: Found Infs and/or NaNs in the "
        "solution after using the fast algorithm.  Retrying using a more "
        "robust algorithm." << std::endl;

      // Restore X from the copy.
      X.assign (B);

      // Find the minimum-norm solution to the least-squares
      // problem $\min_x \|RX - B\|_2$, using the singular value
      // decomposition (SVD).
      LocalDenseMatrixOps<Scalar> ops;
      if (side == LEFT_SIDE) {
        // _GELSS overwrites its matrix input, so make a copy.
        mat_type R_copy (Teuchos::Copy, R_view, N, N);

        // Zero out the lower triangle of R_copy, since the
        // mat_type constructor copies all the entries, not just
        // the upper triangle.  _GELSS will read all the entries
        // of the input matrix.
        ops.zeroOutStrictLowerTriangle (R_copy);

        // Solve the least-squares problem.
        rank = solveLeastSquaresUsingSVD (R_copy, X);
      } else {
        // If solving with R on the right-hand side, the interface
        // requires that instead of solving $\min \|XR - B\|_2$,
        // we have to solve $\min \|R^* X^* - B^*\|_2$.  We
        // compute (conjugate) transposes in newly allocated
        // temporary matrices X_star resp. R_star.  (B is already
        // in X and _GELSS overwrites its input vector X with the
        // solution.)
        mat_type X_star (X.numCols(), X.numRows());
        ops.conjugateTranspose (X_star, X);
        mat_type R_star (N, N); // Filled with zeros automatically.
        ops.conjugateTransposeOfUpperTriangular (R_star, R);

        // Solve the least-squares problem.
        rank = solveLeastSquaresUsingSVD (R_star, X_star);

        // Copy the transpose of X_star back into X.
        ops.conjugateTranspose (X, X_star);
      }
      if (rank < N) {
        foundRankDeficiency = true;
      }
    }
  } else {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, 
           "Should never get here!  Invalid robustness value " 
           << robustness << ".  Please report this bug to the "
           "Belos developers.");
  }
  return std::make_pair (rank, foundRankDeficiency);
      }


    public:
      bool
      testGivensRotations (std::ostream& out)
      {
  using std::endl;

  out << "Testing Givens rotations:" << endl;
  Scalar x = STS::random();
  Scalar y = STS::random();
  out << "  x = " << x << ", y = " << y << endl;

  Scalar theCosine, theSine, result;
  blas_type blas;
  computeGivensRotation (x, y, theCosine, theSine, result);
  out << "-- After computing rotation:" << endl;
  out << "---- cos,sin = " << theCosine << "," << theSine << endl;
  out << "---- x = " << x << ", y = " << y 
      << ", result = " << result << endl;

  blas.ROT (1, &x, 1, &y, 1, &theCosine, &theSine);
  out << "-- After applying rotation:" << endl;
  out << "---- cos,sin = " << theCosine << "," << theSine << endl;
  out << "---- x = " << x << ", y = " << y << endl;

  // Allow only a tiny bit of wiggle room for zeroing-out of y.
  if (STS::magnitude(y) > 2*STS::eps())
    return false;
  else 
    return true;
      }

      bool
      testUpdateColumn (std::ostream& out, 
      const int numCols,
      const bool testBlockGivens=false,
      const bool extraVerbose=false) 
      {
  using Teuchos::Array;
  using std::endl;
  
  TEUCHOS_TEST_FOR_EXCEPTION(numCols <= 0, std::invalid_argument,
         "numCols = " << numCols << " <= 0.");
  const int numRows = numCols + 1;

  mat_type H (numRows, numCols);
  mat_type z (numRows, 1);

  mat_type R_givens (numRows, numCols);
  mat_type y_givens (numRows, 1);
  mat_type z_givens (numRows, 1);
  Array<Scalar> theCosines (numCols);
  Array<Scalar> theSines (numCols);

  mat_type R_blockGivens (numRows, numCols);
  mat_type y_blockGivens (numRows, 1);
  mat_type z_blockGivens (numRows, 1);
  Array<Scalar> blockCosines (numCols);
  Array<Scalar> blockSines (numCols);
  const int panelWidth = std::min (3, numCols);

  mat_type R_lapack (numRows, numCols);
  mat_type y_lapack (numRows, 1);
  mat_type z_lapack (numRows, 1);

  // Make a random least-squares problem.  
  makeRandomProblem (H, z);
  if (extraVerbose) {
    printMatrix<Scalar> (out, "H", H);
    printMatrix<Scalar> (out, "z", z);
  }

  // Set up the right-hand side copies for each of the methods.
  // Each method is free to overwrite its given right-hand side.
  z_givens.assign (z);
  if (testBlockGivens) {
    z_blockGivens.assign (z);
  }
  z_lapack.assign (z);

  //
  // Imitate how one would update the least-squares problem in a
  // typical GMRES implementation, for each updating method.
  //
  // Update using Givens rotations, one at a time.
  magnitude_type residualNormGivens = STM::zero();
  for (int curCol = 0; curCol < numCols; ++curCol) {
    residualNormGivens = updateColumnGivens (H, R_givens, y_givens, z_givens, 
               theCosines, theSines, curCol);
  }
  solveGivens (y_givens, R_givens, z_givens, numCols-1);

  // Update using the "panel left-looking" Givens approach, with
  // the given panel width.
  magnitude_type residualNormBlockGivens = STM::zero();
  if (testBlockGivens) {
    const bool testBlocksAtATime = true;
    if (testBlocksAtATime) {
      // Blocks of columns at a time.
      for (int startCol = 0; startCol < numCols; startCol += panelWidth) {
        int endCol = std::min (startCol + panelWidth - 1, numCols - 1);
        residualNormBlockGivens = 
    updateColumnsGivens (H, R_blockGivens, y_blockGivens, z_blockGivens, 
             blockCosines, blockSines, startCol, endCol);
      }
    } else {
      // One column at a time.  This is good as a sanity check
      // to make sure updateColumnsGivens() with a single column
      // does the same thing as updateColumnGivens().
      for (int startCol = 0; startCol < numCols; ++startCol) {
        residualNormBlockGivens = 
    updateColumnsGivens (H, R_blockGivens, y_blockGivens, z_blockGivens, 
             blockCosines, blockSines, startCol, startCol);
      }
    }
    // The panel version of Givens should compute the same
    // cosines and sines as the non-panel version, and should
    // update the right-hand side z in the same way.  Thus, we
    // should be able to use the same triangular solver.
    solveGivens (y_blockGivens, R_blockGivens, z_blockGivens, numCols-1);
  }

  // Solve using LAPACK's least-squares solver.
  const magnitude_type residualNormLapack = 
    solveLapack (H, R_lapack, y_lapack, z_lapack, numCols-1);

  // Compute the condition number of the least-squares problem.
  // This requires a residual, so use the residual from the
  // LAPACK method.  All that the method needs for an accurate
  // residual norm is forward stability.
  const magnitude_type leastSquaresCondNum = 
    leastSquaresConditionNumber (H, z, residualNormLapack);

  // Compute the relative least-squares solution error for both
  // Givens methods.  We assume that the LAPACK solution is
  // "exact" and compare against the Givens rotations solution.
  // This is taking liberties with the definition of condition
  // number, but it's the best we can do, since we don't know
  // the exact solution and don't have an extended-precision
  // solver.

  // The solution lives only in y[0 .. numCols-1].
  mat_type y_givens_view (Teuchos::View, y_givens, numCols, 1);
  mat_type y_blockGivens_view (Teuchos::View, y_blockGivens, numCols, 1);
  mat_type y_lapack_view (Teuchos::View, y_lapack, numCols, 1);

  const magnitude_type givensSolutionError = 
    solutionError (y_givens_view, y_lapack_view);
  const magnitude_type blockGivensSolutionError = testBlockGivens ?
    solutionError (y_blockGivens_view, y_lapack_view) : 
    STM::zero();

  // If printing out the matrices, copy out the upper triangular
  // factors for printing.  (Both methods are free to leave data
  // below the lower triangle.)
  if (extraVerbose) {
    mat_type R_factorFromGivens (numCols, numCols);
    mat_type R_factorFromBlockGivens (numCols, numCols);
    mat_type R_factorFromLapack (numCols, numCols);

    for (int j = 0; j < numCols; ++j) {
      for (int i = 0; i <= j; ++i) {
        R_factorFromGivens(i,j) = R_givens(i,j);
        if (testBlockGivens) {
    R_factorFromBlockGivens(i,j) = R_blockGivens(i,j);
        }
        R_factorFromLapack(i,j) = R_lapack(i,j);
      }
    }

    printMatrix<Scalar> (out, "R_givens", R_factorFromGivens);
    printMatrix<Scalar> (out, "y_givens", y_givens_view);
    printMatrix<Scalar> (out, "z_givens", z_givens);
      
    if (testBlockGivens) {
      printMatrix<Scalar> (out, "R_blockGivens", R_factorFromBlockGivens);
      printMatrix<Scalar> (out, "y_blockGivens", y_blockGivens_view);
      printMatrix<Scalar> (out, "z_blockGivens", z_blockGivens);
    }
      
    printMatrix<Scalar> (out, "R_lapack", R_factorFromLapack);
    printMatrix<Scalar> (out, "y_lapack", y_lapack_view);
    printMatrix<Scalar> (out, "z_lapack", z_lapack);
  }

  // Compute the (Frobenius) norm of the original matrix H.
  const magnitude_type H_norm = H.normFrobenius();

  out << "||H||_F = " << H_norm << endl;

  out << "||H y_givens - z||_2 / ||H||_F = " 
      << leastSquaresResidualNorm (H, y_givens_view, z) / H_norm << endl;
  if (testBlockGivens) {
    out << "||H y_blockGivens - z||_2 / ||H||_F = " 
        << leastSquaresResidualNorm (H, y_blockGivens_view, z) / H_norm << endl;
  }
  out << "||H y_lapack - z||_2 / ||H||_F = " 
      << leastSquaresResidualNorm (H, y_lapack_view, z) / H_norm << endl;

  out << "||y_givens - y_lapack||_2 / ||y_lapack||_2 = " 
      << givensSolutionError << endl;
  if (testBlockGivens) {
    out << "||y_blockGivens - y_lapack||_2 / ||y_lapack||_2 = " 
        << blockGivensSolutionError << endl;
  }

  out << "Least-squares condition number = " 
      << leastSquaresCondNum << endl;

  // Now for the controversial part of the test: judging whether
  // we succeeded.  This includes the problem's condition
  // number, which is a measure of the maximum perturbation in
  // the solution for which we can still say that the solution
  // is valid.  We include a little wiggle room by including a
  // factor proportional to the square root of the number of
  // floating-point operations that influence the last entry
  // (the conventional Wilkinsonian heuristic), times 10 for
  // good measure.
  //
  // (The square root looks like it has something to do with an
  // average-case probabilistic argument, but doesn't really.
  // What's an "average problem"?)
  const magnitude_type wiggleFactor = 
    10 * STM::squareroot( numRows*numCols );
  const magnitude_type solutionErrorBoundFactor = 
    wiggleFactor * leastSquaresCondNum;
  const magnitude_type solutionErrorBound = 
    solutionErrorBoundFactor * STS::eps();
  out << "Solution error bound: " << solutionErrorBoundFactor 
      << " * eps = " << solutionErrorBound << endl;

  // Remember that NaN is not greater than, not less than, and
  // not equal to any other number, including itself.  Some
  // compilers will rudely optimize away the "x != x" test.
  if (STM::isnaninf (solutionErrorBound)) {
    // Hm, the solution error bound is Inf or NaN.  This
    // probably means that the test problem was generated
    // incorrectly.  We should return false in this case.
    return false; 
  } else { // solution error bound is finite.
    if (STM::isnaninf (givensSolutionError)) {
      return false; 
    } else if (givensSolutionError > solutionErrorBound) {
      return false;
    } else if (testBlockGivens) {
      if (STM::isnaninf (blockGivensSolutionError)) {
        return false; 
      } else if (blockGivensSolutionError > solutionErrorBound) {
        return false;
      } else { // Givens and Block Givens tests succeeded.
        return true;
      }
    } else { // Not testing block Givens; Givens test succeeded.
      return true;
    }
  }
      }

      bool 
      testTriangularSolves (std::ostream& out,
          const int testProblemSize,
          const ERobustness robustness,
          const bool verbose=false)
      {
  using Teuchos::LEFT_SIDE;
  using Teuchos::RIGHT_SIDE;
  using std::endl;
  typedef Teuchos::SerialDenseMatrix<int, scalar_type> mat_type;

  Teuchos::oblackholestream blackHole;
  std::ostream& verboseOut = verbose ? out : blackHole;

  verboseOut << "Testing upper triangular solves" << endl;
  //
  // Construct an upper triangular linear system to solve.
  //
  verboseOut << "-- Generating test matrix" << endl;
  const int N = testProblemSize;
  mat_type R (N, N);
  // Fill the upper triangle of R with random numbers.
  for (int j = 0; j < N; ++j) {
    for (int i = 0; i <= j; ++i) {
      R(i,j) = STS::random ();
    }
  }
  // Compute the Frobenius norm of R for later use.
  const magnitude_type R_norm = R.normFrobenius ();
  // Fill the right-hand side B with random numbers.
  mat_type B (N, 1);
  B.random ();
  // Compute the Frobenius norm of B for later use.
  const magnitude_type B_norm = B.normFrobenius ();

  // Save a copy of the original upper triangular system.
  mat_type R_copy (Teuchos::Copy, R, N, N);
  mat_type B_copy (Teuchos::Copy, B, N, 1);

  // Solution vector.
  mat_type X (N, 1);

  // Solve RX = B.
  verboseOut << "-- Solving RX=B" << endl;
  // We're ignoring the return values for now.
  (void) solveUpperTriangularSystem (LEFT_SIDE, X, R, B, robustness);
  // Test the residual error.
  mat_type Resid (N, 1);
  Resid.assign (B_copy);
  Belos::details::LocalDenseMatrixOps<scalar_type> ops;
  ops.matMatMult (STS::one(), Resid, -STS::one(), R_copy, X);
  verboseOut << "---- ||R*X - B||_F = " << Resid.normFrobenius() << endl;
  verboseOut << "---- ||R||_F ||X||_F + ||B||_F = " 
       << (R_norm * X.normFrobenius() + B_norm)
       << endl;

  // Restore R and B.
  R.assign (R_copy);
  B.assign (B_copy);

  // 
  // Set up a right-side test problem: YR = B^*.
  //
  mat_type Y (1, N);
  mat_type B_star (1, N);
  ops.conjugateTranspose (B_star, B);
  mat_type B_star_copy (1, N);
  B_star_copy.assign (B_star);
  // Solve YR = B^*.
  verboseOut << "-- Solving YR=B^*" << endl;
  // We're ignoring the return values for now.
  (void) solveUpperTriangularSystem (RIGHT_SIDE, Y, R, B_star, robustness);
  // Test the residual error.
  mat_type Resid2 (1, N);
  Resid2.assign (B_star_copy);
  ops.matMatMult (STS::one(), Resid2, -STS::one(), Y, R_copy);
  verboseOut << "---- ||Y*R - B^*||_F = " << Resid2.normFrobenius() << endl;
  verboseOut << "---- ||Y||_F ||R||_F + ||B^*||_F = " 
       << (Y.normFrobenius() * R_norm + B_norm)
       << endl;

  // FIXME (mfh 14 Oct 2011) The test always "passes" for now;
  // you have to inspect the output in order to see whether it
  // succeeded.  We really should fix the above to use the
  // infinity-norm bounds in Higham's book for triangular
  // solves.  That would automate the test.
  return true;
      }

    private:
      std::ostream& warn_;

      ERobustness defaultRobustness_;

    private:
      int 
      solveLeastSquaresUsingSVD (mat_type& A, mat_type& X)
      {
  using Teuchos::Array;
  LocalDenseMatrixOps<Scalar> ops;

  if (defaultRobustness_ > ROBUSTNESS_SOME) {
    int count = ops.infNaNCount (A);
    TEUCHOS_TEST_FOR_EXCEPTION(count != 0, std::invalid_argument, 
           "solveLeastSquaresUsingSVD: The input matrix A "
           "contains " << count << "Inf and/or NaN entr" 
           << (count != 1 ? "ies" : "y") << ".");
    count = ops.infNaNCount (X);
    TEUCHOS_TEST_FOR_EXCEPTION(count != 0, std::invalid_argument, 
           "solveLeastSquaresUsingSVD: The input matrix X "
           "contains " << count << "Inf and/or NaN entr" 
           << (count != 1 ? "ies" : "y") << ".");
  }
  const int N = std::min (A.numRows(), A.numCols());
  lapack_type lapack;

  // Rank of A; to be computed by _GELSS and returned.
  int rank = N;

  // Use Scalar's machine precision for the rank tolerance,
  // not magnitude_type's machine precision.
  const magnitude_type rankTolerance = STS::eps();

  // Array of singular values.
  Array<magnitude_type> singularValues (N);

  // Extra workspace.  This is only used by _GELSS if Scalar is
  // complex.  Teuchos::LAPACK presents a unified interface to
  // _GELSS that always includes the RWORK argument, even though
  // LAPACK's SGELSS and DGELSS don't have the RWORK argument.
  // We always allocate at least one entry so that &rwork[0]
  // makes sense.
  Array<magnitude_type> rwork (1);
  if (STS::isComplex) {
    rwork.resize (std::max (1, 5 * N));
  }
  //
  // Workspace query
  //
  Scalar lworkScalar = STS::one(); // To be set by workspace query
  int info = 0;
  lapack.GELSS (A.numRows(), A.numCols(), X.numCols(), 
          A.values(), A.stride(), X.values(), X.stride(),
          &singularValues[0], rankTolerance, &rank,
          &lworkScalar, -1, &rwork[0], &info);
  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
         "_GELSS workspace query returned INFO = " 
         << info << " != 0.");
  const int lwork = static_cast<int> (STS::real (lworkScalar));
  TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
         "_GELSS workspace query returned LWORK = " 
         << lwork << " < 0.");
  // Allocate workspace.  Size > 0 means &work[0] makes sense.
  Array<Scalar> work (std::max (1, lwork));
  // Solve the least-squares problem.
  lapack.GELSS (A.numRows(), A.numCols(), X.numCols(), 
          A.values(), A.stride(), X.values(), X.stride(),
          &singularValues[0], rankTolerance, &rank,
          &work[0], lwork, &rwork[0], &info);
  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
         "_GELSS returned INFO = " << info << " != 0.");
  return rank;
      }

      void
      solveGivens (mat_type& y, 
       mat_type& R, 
       const mat_type& z, 
       const int curCol)
      {
  const int numRows = curCol + 2;

  // Now that we have the updated R factor of H, and the updated
  // right-hand side z, solve the least-squares problem by
  // solving the upper triangular linear system Ry=z for y.
  const mat_type R_view (Teuchos::View, R, numRows-1, numRows-1);
  const mat_type z_view (Teuchos::View, z, numRows-1, z.numCols());
  mat_type y_view (Teuchos::View, y, numRows-1, y.numCols());

  (void) solveUpperTriangularSystem (Teuchos::LEFT_SIDE, y_view, 
             R_view, z_view, defaultRobustness_);
      }

      void
      makeRandomProblem (mat_type& H, mat_type& z) 
      {
  // In GMRES, z always starts out with only the first entry
  // being nonzero.  That entry always has nonnegative real part
  // and zero imaginary part, since it is the initial residual
  // norm.
  H.random ();
  // Zero out the entries below the subdiagonal of H, so that it
  // is upper Hessenberg.
  for (int j = 0; j < H.numCols(); ++j) {
    for (int i = j+2; i < H.numRows(); ++i) {
      H(i,j) = STS::zero();
    }
  }
  // Initialize z, the right-hand side of the least-squares
  // problem.  Make the first entry of z nonzero.
  {
    // It's still possible that a random number will come up
    // zero after 1000 trials, but unlikely.  Nevertheless, it's
    // still important not to allow an infinite loop, for
    // example if the pseudorandom number generator is broken
    // and always returns zero.
    const int numTrials = 1000;
    magnitude_type z_init = STM::zero();
    for (int trial = 0; trial < numTrials && z_init == STM::zero(); ++trial) {
      z_init = STM::random();
    }
    TEUCHOS_TEST_FOR_EXCEPTION(z_init == STM::zero(), std::runtime_error,
           "After " << numTrials << " trial" 
           << (numTrials != 1 ? "s" : "") 
           << ", we were unable to generate a nonzero pseudo"
           "random real number.  This most likely indicates a "
           "broken pseudorandom number generator.");
    const magnitude_type z_first = (z_init < 0) ? -z_init : z_init;

    // NOTE I'm assuming here that "scalar_type = magnitude_type"
    // assignments make sense.
    z(0,0) = z_first;
  }
      }

      void
      computeGivensRotation (const Scalar& x, 
           const Scalar& y, 
           Scalar& theCosine, 
           Scalar& theSine,
           Scalar& result) 
      {
  // _LARTG, an LAPACK aux routine, is slower but more accurate
  // than the BLAS' _ROTG.
  const bool useLartg = false;

  if (useLartg) {
    lapack_type lapack;
    // _LARTG doesn't clobber its input arguments x and y.
    lapack.LARTG (x, y, &theCosine, &theSine, &result);
  } else {
    // _ROTG clobbers its first two arguments.  x is overwritten
    // with the result of applying the Givens rotation: [x; y] ->
    // [x (on output); 0].  y is overwritten with the "fast"
    // Givens transform (see Golub and Van Loan, 3rd ed.).
    Scalar x_temp = x;
    Scalar y_temp = y;
    blas_type blas;
    blas.ROTG (&x_temp, &y_temp, &theCosine, &theSine);
    result = x_temp;
  }
      }

      void
      singularValues (const mat_type& A, 
          Teuchos::ArrayView<magnitude_type> sigmas) 
      {
  using Teuchos::Array;
  using Teuchos::ArrayView;

  const int numRows = A.numRows();
  const int numCols = A.numCols();
  TEUCHOS_TEST_FOR_EXCEPTION(sigmas.size() < std::min (numRows, numCols),
         std::invalid_argument,
         "The sigmas array is only of length " << sigmas.size()
         << ", but must be of length at least "
         << std::min (numRows, numCols)
         << " in order to hold all the singular values of the "
         "matrix A.");

  // Compute the condition number of the matrix A, using a singular
  // value decomposition (SVD).  LAPACK's SVD routine overwrites the
  // input matrix, so make a copy.
  mat_type A_copy (numRows, numCols);
  A_copy.assign (A);

  // Workspace query.
  lapack_type lapack;
  int info = 0;
  Scalar lworkScalar = STS::zero();
  Array<magnitude_type> rwork (std::max (std::min (numRows, numCols) - 1, 1));
  lapack.GESVD ('N', 'N', numRows, numCols, 
          A_copy.values(), A_copy.stride(), &sigmas[0], 
          (Scalar*) NULL, 1, (Scalar*) NULL, 1, 
          &lworkScalar, -1, &rwork[0], &info);

  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
         "LAPACK _GESVD workspace query failed with INFO = " 
         << info << ".");
  const int lwork = static_cast<int> (STS::real (lworkScalar));
  TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
         "LAPACK _GESVD workspace query returned LWORK = " 
         << lwork << " < 0.");
  // Make sure that the workspace array always has positive
  // length, so that &work[0] makes sense.
  Teuchos::Array<Scalar> work (std::max (1, lwork));

  // Compute the singular values of A.
  lapack.GESVD ('N', 'N', numRows, numCols, 
          A_copy.values(), A_copy.stride(), &sigmas[0], 
          (Scalar*) NULL, 1, (Scalar*) NULL, 1, 
          &work[0], lwork, &rwork[0], &info);
  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
         "LAPACK _GESVD failed with INFO = " << info << ".");
      }

      std::pair<magnitude_type, magnitude_type>
      extremeSingularValues (const mat_type& A) 
      {
  using Teuchos::Array;

  const int numRows = A.numRows();
  const int numCols = A.numCols();

  Array<magnitude_type> sigmas (std::min (numRows, numCols));
  singularValues (A, sigmas);
  return std::make_pair (sigmas[0], sigmas[std::min(numRows, numCols) - 1]);
      }
      
      magnitude_type 
      leastSquaresConditionNumber (const mat_type& A,
           const mat_type& b,
           const magnitude_type& residualNorm) 
      {
  // Extreme singular values of A.
  const std::pair<magnitude_type, magnitude_type> sigmaMaxMin = 
    extremeSingularValues (A);

  // Our solvers currently assume that H has full rank.  If the
  // test matrix doesn't have full rank, we stop right away.
  TEUCHOS_TEST_FOR_EXCEPTION(sigmaMaxMin.second == STM::zero(), std::runtime_error,
         "The test matrix is rank deficient; LAPACK's _GESVD "
         "routine reports that its smallest singular value is "
         "zero.");
  // 2-norm condition number of A.  We checked above that the
  // denominator is nonzero.
  const magnitude_type A_cond = sigmaMaxMin.first / sigmaMaxMin.second;

  // "Theta" in the variable names below refers to the angle between
  // the vectors b and A*x, where x is the computed solution.  It
  // measures whether the residual norm is large (near ||b||) or
  // small (near 0).
  const magnitude_type sinTheta = residualNorm / b.normFrobenius();

  // \sin^2 \theta + \cos^2 \theta = 1.
  //
  // The range of sine is [-1,1], so squaring it won't overflow.
  // We still have to check whether sinTheta > 1, though.  This
  // is impossible in exact arithmetic, assuming that the
  // least-squares solver worked (b-A*0 = b and x minimizes
  // ||b-A*x||_2, so ||b-A*0||_2 >= ||b-A*x||_2).  However, it
  // might just be possible in floating-point arithmetic.  We're
  // just looking for an estimate, so if sinTheta > 1, we cap it
  // at 1.
  const magnitude_type cosTheta = (sinTheta > STM::one()) ? 
    STM::zero() : STM::squareroot (1 - sinTheta * sinTheta);

  // This may result in Inf, if cosTheta is zero.  That's OK; in
  // that case, the condition number of the (full-rank)
  // least-squares problem is rightfully infinite.
  const magnitude_type tanTheta = sinTheta / cosTheta;

  // Condition number for the full-rank least-squares problem.
  return 2 * A_cond / cosTheta + tanTheta * A_cond * A_cond;
      }
      
      magnitude_type
      leastSquaresResidualNorm (const mat_type& A,
        const mat_type& x,
        const mat_type& b) 
      {
  mat_type r (b.numRows(), b.numCols());

  // r := b - A*x
  r.assign (b);
  LocalDenseMatrixOps<Scalar> ops;
  ops.matMatMult (STS::one(), r, -STS::one(), A, x);
  return r.normFrobenius ();
      }

      magnitude_type
      solutionError (const mat_type& x_approx,
         const mat_type& x_exact) 
      {
  const int numRows = x_exact.numRows();
  const int numCols = x_exact.numCols();

  mat_type x_diff (numRows, numCols);
  for (int j = 0; j < numCols; ++j) {
    for (int i = 0; i < numRows; ++i) {
      x_diff(i,j) = x_exact(i,j) - x_approx(i,j);
    }
  }
  const magnitude_type scalingFactor = x_exact.normFrobenius();

  // If x_exact has zero norm, just use the absolute difference.
  return x_diff.normFrobenius() / 
    (scalingFactor == STM::zero() ? STM::one() : scalingFactor);
      }

      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) 
      {
  using std::cerr;
  using std::endl;

  const int numRows = curCol + 2; // curCol is zero-based
  const int LDR = R.stride();
  const bool extraDebug = false;

  if (extraDebug) {
    cerr << "updateColumnGivens: curCol = " << curCol << endl;
  }

  // View of H( 1:curCol+1, curCol ) (in Matlab notation, if
  // curCol were a one-based index, as it would be in Matlab but
  // is not here).
  const mat_type H_col (Teuchos::View, H, numRows, 1, 0, curCol);

  // View of R( 1:curCol+1, curCol ) (again, in Matlab notation,
  // if curCol were a one-based index).
  mat_type R_col (Teuchos::View, R, numRows, 1, 0, curCol);

  // 1. Copy the current column from H into R, where it will be
  //    modified.
  R_col.assign (H_col);

  if (extraDebug) {
    printMatrix<Scalar> (cerr, "R_col before ", R_col);
  }

  // 2. Apply all the previous Givens rotations, if any, to the
  //    current column of the matrix.
  blas_type blas;
  for (int j = 0; j < curCol; ++j) {
    // BLAS::ROT really should take "const Scalar*" for these
    // arguments, but it wants a "Scalar*" instead, alas.
    Scalar theCosine = theCosines[j];
    Scalar theSine = theSines[j];
    
    if (extraDebug) {
      cerr << "  j = " << j << ": (cos,sin) = " 
     << theCosines[j] << "," << theSines[j] << endl;
    }
    blas.ROT (1, &R_col(j,0), LDR, &R_col(j+1,0), LDR,
        &theCosine, &theSine);
  }
  if (extraDebug && curCol > 0) {
    printMatrix<Scalar> (cerr, "R_col after applying previous "
             "Givens rotations", R_col);
  }

  // 3. Calculate new Givens rotation for R(curCol, curCol),
  //    R(curCol+1, curCol).
  Scalar theCosine, theSine, result;
  computeGivensRotation (R_col(curCol,0), R_col(curCol+1,0), 
             theCosine, theSine, result);
  theCosines[curCol] = theCosine;
  theSines[curCol] = theSine;

  if (extraDebug) {
    cerr << "  New cos,sin = " << theCosine << "," << theSine << endl;
  }

  // 4. _Apply_ the new Givens rotation.  We don't need to
  //    invoke _ROT here, because computeGivensRotation()
  //    already gives us the result: [x; y] -> [result; 0].
  R_col(curCol, 0) = result;
  R_col(curCol+1, 0) = STS::zero();

  if (extraDebug) {
    printMatrix<Scalar> (cerr, "R_col after applying current "
             "Givens rotation", R_col);
  }

  // 5. Apply the resulting Givens rotation to z (the right-hand
  //    side of the projected least-squares problem).
  //
  // We prefer overgeneralization to undergeneralization by assuming
  // here that z may have more than one column.
  const int LDZ = z.stride();
  blas.ROT (z.numCols(), 
      &z(curCol,0), LDZ, &z(curCol+1,0), LDZ,
      &theCosine, &theSine);

  if (extraDebug) {
    //mat_type R_after (Teuchos::View, R, numRows, numRows-1);
    //printMatrix<Scalar> (cerr, "R_after", R_after);
    //mat_type z_after (Teuchos::View, z, numRows, z.numCols());
    printMatrix<Scalar> (cerr, "z_after", z);
  }

  // The last entry of z is the nonzero part of the residual of the
  // least-squares problem.  Its magnitude gives the residual 2-norm
  // of the least-squares problem.
  return STS::magnitude( z(numRows-1, 0) );
      }

      magnitude_type
      solveLapack (const mat_type& H, 
       mat_type& R, 
       mat_type& y, 
       mat_type& z,
       const int curCol) 
      {
  const int numRows = curCol + 2;
  const int numCols = curCol + 1;
  const int LDR = R.stride();

  // Copy H( 1:curCol+1, 1:curCol ) into R( 1:curCol+1, 1:curCol ).
  const mat_type H_view (Teuchos::View, H, numRows, numCols);
  mat_type R_view (Teuchos::View, R, numRows, numCols);
  R_view.assign (H_view);

  // The LAPACK least-squares solver overwrites the right-hand side
  // vector with the solution, so first copy z into y.
  mat_type y_view (Teuchos::View, y, numRows, y.numCols());
  mat_type z_view (Teuchos::View, z, numRows, y.numCols());
  y_view.assign (z_view);

  // Workspace query for the least-squares routine.
  int info = 0;
  Scalar lworkScalar = STS::zero();
  lapack_type lapack;
  lapack.GELS ('N', numRows, numCols, y_view.numCols(),
         NULL, LDR, NULL, y_view.stride(), 
         &lworkScalar, -1, &info);
  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
         "LAPACK _GELS workspace query failed with INFO = " 
         << info << ", for a " << numRows << " x " << numCols
         << " matrix with " << y_view.numCols() 
         << " right hand side"
         << ((y_view.numCols() != 1) ? "s" : "") << ".");
  TEUCHOS_TEST_FOR_EXCEPTION(STS::real(lworkScalar) < STM::zero(),
         std::logic_error,
         "LAPACK _GELS workspace query returned an LWORK with "
         "negative real part: LWORK = " << lworkScalar 
         << ".  That should never happen.  Please report this "
         "to the Belos developers.");
  TEUCHOS_TEST_FOR_EXCEPTION(STS::isComplex && STS::imag(lworkScalar) != STM::zero(),
         std::logic_error,
         "LAPACK _GELS workspace query returned an LWORK with "
         "nonzero imaginary part: LWORK = " << lworkScalar 
         << ".  That should never happen.  Please report this "
         "to the Belos developers.");
  // Cast workspace from Scalar to int.  Scalar may be complex,
  // hence the request for the real part.  Don't ask for the
  // magnitude, since computing the magnitude may overflow due
  // to squaring and square root to int.  Hopefully LAPACK
  // doesn't ever overflow int this way.
  const int lwork = std::max (1, static_cast<int> (STS::real (lworkScalar)));

  // Allocate workspace for solving the least-squares problem.
  Teuchos::Array<Scalar> work (lwork);

  // Solve the least-squares problem.  The ?: operator prevents
  // accessing the first element of the work array, if it has
  // length zero.
  lapack.GELS ('N', numRows, numCols, y_view.numCols(),
         R_view.values(), R_view.stride(),
         y_view.values(), y_view.stride(),
         (lwork > 0 ? &work[0] : (Scalar*) NULL), 
         lwork, &info);

  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
         "Solving projected least-squares problem with LAPACK "
         "_GELS failed with INFO = " << info << ", for a " 
         << numRows << " x " << numCols << " matrix with " 
         << y_view.numCols() << " right hand side"
         << (y_view.numCols() != 1 ? "s" : "") << ".");
  // Extract the projected least-squares problem's residual error.
  // It's the magnitude of the last entry of y_view on output from
  // LAPACK's least-squares solver.  
  return STS::magnitude( y_view(numRows-1, 0) );  
      }

      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) 
      {
  TEUCHOS_TEST_FOR_EXCEPTION(startCol > endCol, std::invalid_argument,
         "updateColumnGivens: startCol = " << startCol 
         << " > endCol = " << endCol << ".");
  magnitude_type lastResult = STM::zero();
  // [startCol, endCol] is an inclusive range.
  for (int curCol = startCol; curCol <= endCol; ++curCol) {
    lastResult = updateColumnGivens (H, R, y, z, theCosines, theSines, curCol);
  }
  return lastResult;
      }

      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) 
      {
  const int numRows = endCol + 2;
  const int numColsToUpdate = endCol - startCol + 1;
  const int LDR = R.stride();

  // 1. Copy columns [startCol, endCol] from H into R, where they
  //    will be modified.
  {
    const mat_type H_view (Teuchos::View, H, numRows, numColsToUpdate, 0, startCol);
    mat_type R_view (Teuchos::View, R, numRows, numColsToUpdate, 0, startCol);
    R_view.assign (H_view);
  }

  // 2. Apply all the previous Givens rotations, if any, to
  //    columns [startCol, endCol] of the matrix.  (Remember
  //    that we're using a left-looking QR factorization
  //    approach; we haven't yet touched those columns.)
  blas_type blas;
  for (int j = 0; j < startCol; ++j) {
    blas.ROT (numColsToUpdate,
        &R(j, startCol), LDR, &R(j+1, startCol), LDR, 
        &theCosines[j], &theSines[j]);
  }

  // 3. Update each column in turn of columns [startCol, endCol].
  for (int curCol = startCol; curCol < endCol; ++curCol) {
    // a. Apply the Givens rotations computed in previous
    //    iterations of this loop to the current column of R.
    for (int j = startCol; j < curCol; ++j) {
      blas.ROT (1, &R(j, curCol), LDR, &R(j+1, curCol), LDR, 
          &theCosines[j], &theSines[j]);
    }
    // b. Calculate new Givens rotation for R(curCol, curCol),
    //    R(curCol+1, curCol).
    Scalar theCosine, theSine, result;
    computeGivensRotation (R(curCol, curCol), R(curCol+1, curCol), 
         theCosine, theSine, result);
    theCosines[curCol] = theCosine;
    theSines[curCol] = theSine;

    // c. _Apply_ the new Givens rotation.  We don't need to
    //    invoke _ROT here, because computeGivensRotation()
    //    already gives us the result: [x; y] -> [result; 0].
    R(curCol+1, curCol) = result;
    R(curCol+1, curCol) = STS::zero();

    // d. Apply the resulting Givens rotation to z (the right-hand
    //    side of the projected least-squares problem).  
    //
    // We prefer overgeneralization to undergeneralization by
    // assuming here that z may have more than one column.
    const int LDZ = z.stride();
    blas.ROT (z.numCols(), 
        &z(curCol,0), LDZ, &z(curCol+1,0), LDZ,
        &theCosine, &theSine);
  }

  // The last entry of z is the nonzero part of the residual of the
  // least-squares problem.  Its magnitude gives the residual 2-norm
  // of the least-squares problem.
  return STS::magnitude( z(numRows-1, 0) );
      }
    }; // class ProjectedLeastSquaresSolver
  } // namespace details
} // namespace Belos

#endif // __Belos_ProjectedLeastSquaresSolver_hpp