trilinos/teuchos/Teuchos__SerialQRDenseSolver_8hpp_source.html

// @HEADER

// *****************************************************************************

//                    Teuchos: Common Tools Package

//

// Copyright 2004 NTESS and the Teuchos contributors.

// SPDX-License-Identifier: BSD-3-Clause

// *****************************************************************************

// @HEADER


#ifndef _TEUCHOS_SERIALQRDENSESOLVER_HPP_

#define _TEUCHOS_SERIALQRDENSESOLVER_HPP_


#include "Teuchos_CompObject.hpp"

#include "Teuchos_BLAS.hpp"

#include "Teuchos_LAPACK.hpp"

#include "Teuchos_RCP.hpp"

#include "Teuchos_ConfigDefs.hpp"

#include "Teuchos_SerialDenseMatrix.hpp"

#include "Teuchos_SerialDenseSolver.hpp"

#include "Teuchos_ScalarTraits.hpp"


#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

#include "Eigen/Dense"

#endif


namespace Teuchos {


  template<typename OrdinalType, typename ScalarType>


  class SerialQRDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,

                              public LAPACK<OrdinalType, ScalarType>

  {


  public:


    typedef typename ScalarTraits<ScalarType>::magnitudeType MagnitudeType;

#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

    typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;

    typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;

    typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;

    typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;

    typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;

    typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;

    typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;

    typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;

    typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;

    typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;

    typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;

    typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;

    typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;

#endif


    SerialQRDenseSolver();


    virtual ~SerialQRDenseSolver();


    int setMatrix(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);


    int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,

                   const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);


    void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}


    void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::CONJ_TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}


    void solveWithTransposeFlag(Teuchos::ETransp trans) {TRANS_ = trans; transpose_ = (trans != Teuchos::NO_TRANS) ? true : false; }


    int factor();


    int solve();


    int computeEquilibrateScaling();


    int equilibrateMatrix();


    int equilibrateRHS();


    int unequilibrateLHS();


    int formQ();


    int formR();


    int multiplyQ (ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C);


    int solveR (ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C);


    bool transpose() {return(transpose_);}


    bool factored() {return(factored_);}


    bool equilibratedA() {return(equilibratedA_);}


    bool equilibratedB() {return(equilibratedB_);}


    bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}


    bool solved() {return(solved_);}


    bool formedQ() {return(formedQ_);}


    bool formedR() {return(formedR_);}


    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getMatrix()  const {return(Matrix_);}


    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix()  const {return(Factor_);}


    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getQ()  const {return(FactorQ_);}


    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getR()  const {return(FactorR_);}


    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}


    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}


    OrdinalType numRows()  const {return(M_);}


    OrdinalType numCols()  const {return(N_);}


    std::vector<ScalarType> tau()  const {return(TAU_);}


    MagnitudeType ANORM()  const {return(ANORM_);}


    void Print(std::ostream& os) const;

  protected:


    void allocateWORK() { LWORK_ = 4*N_; WORK_.resize( LWORK_ ); return;}

    void resetMatrix();

    void resetVectors();


    bool equilibrate_;

    bool shouldEquilibrate_;

    bool equilibratedA_;

    bool equilibratedB_;

    OrdinalType equilibrationModeA_;

    OrdinalType equilibrationModeB_;

    bool transpose_;

    bool factored_;

    bool solved_;

    bool matrixIsZero_;

    bool formedQ_;

    bool formedR_;


    Teuchos::ETransp TRANS_;


    OrdinalType M_;

    OrdinalType N_;

    OrdinalType Min_MN_;

    OrdinalType LDA_;

    OrdinalType LDAF_;

    OrdinalType LDQ_;

    OrdinalType LDR_;

    OrdinalType INFO_;

    OrdinalType LWORK_;


    std::vector<ScalarType> TAU_;


    MagnitudeType ANORM_;

    MagnitudeType BNORM_;


    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Matrix_;

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > TMP_;

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Factor_;

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > FactorQ_;

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > FactorR_;


    ScalarType * A_;

    ScalarType * AF_;

    ScalarType * Q_;

    ScalarType * R_;

    std::vector<ScalarType> WORK_;

#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

    Eigen::HouseholderQR<EigenMatrix> qr_;

#endif


  private:

    // SerialQRDenseSolver copy constructor (put here because we don't want user access)


    SerialQRDenseSolver(const SerialQRDenseSolver<OrdinalType, ScalarType>& Source);

    SerialQRDenseSolver & operator=(const SerialQRDenseSolver<OrdinalType, ScalarType>& Source);


  };


  // Helper traits to distinguish work arrays for real and complex-valued datatypes.

  using namespace details;


//=============================================================================


template<typename OrdinalType, typename ScalarType>


SerialQRDenseSolver<OrdinalType,ScalarType>::SerialQRDenseSolver()

  : CompObject(),

    equilibrate_(false),

    shouldEquilibrate_(false),

    equilibratedA_(false),

    equilibratedB_(false),

    equilibrationModeA_(0),

    equilibrationModeB_(0),

    transpose_(false),

    factored_(false),

    solved_(false),

    matrixIsZero_(false),

    formedQ_(false),

    formedR_(false),

    TRANS_(Teuchos::NO_TRANS),

    M_(0),

    N_(0),

    Min_MN_(0),

    LDA_(0),

    LDAF_(0),

    LDQ_(0),

    LDR_(0),

    INFO_(0),

    LWORK_(0),

    ANORM_(ScalarTraits<MagnitudeType>::zero()),

    BNORM_(ScalarTraits<MagnitudeType>::zero()),

    A_(0),

    AF_(0),

    Q_(0),

    R_(0)

{

  resetMatrix();

}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


SerialQRDenseSolver<OrdinalType,ScalarType>::~SerialQRDenseSolver()

{}


//=============================================================================


template<typename OrdinalType, typename ScalarType>

void SerialQRDenseSolver<OrdinalType,ScalarType>::resetVectors()

{

  LHS_ = Teuchos::null;

  TMP_ = Teuchos::null;

  RHS_ = Teuchos::null;

  solved_ = false;

  equilibratedB_ = false;

}

//=============================================================================


template<typename OrdinalType, typename ScalarType>

void SerialQRDenseSolver<OrdinalType,ScalarType>::resetMatrix()

{

  resetVectors();

  equilibratedA_ = false;

  equilibrationModeA_ = 0;

  equilibrationModeB_ = 0;

  factored_ = false;

  matrixIsZero_ = false;

  formedQ_ = false;

  formedR_ = false;

  M_ = 0;

  N_ = 0;

  Min_MN_ = 0;

  LDA_ = 0;

  LDAF_ = 0;

  LDQ_ = 0;

  LDR_ = 0;

  ANORM_ = -ScalarTraits<MagnitudeType>::one();

  BNORM_ = -ScalarTraits<MagnitudeType>::one();

  A_ = 0;

  AF_ = 0;

  Q_ = 0;

  R_ = 0;

  INFO_ = 0;

  LWORK_ = 0;

}

//=============================================================================


template<typename OrdinalType, typename ScalarType>


int SerialQRDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)

{

  TEUCHOS_TEST_FOR_EXCEPTION(A->numRows() < A->numCols(), std::invalid_argument,

                     "SerialQRDenseSolver<T>::setMatrix: the matrix A must have A.numRows() >= A.numCols()!");


  resetMatrix();

  Matrix_ = A;

  Factor_ = A;

  FactorQ_ = A;

  FactorR_ = A;

  M_ = A->numRows();

  N_ = A->numCols();

  Min_MN_ = TEUCHOS_MIN(M_,N_);

  LDA_ = A->stride();

  LDAF_ = LDA_;

  LDQ_ = LDA_;

  LDR_ = N_;

  A_ = A->values();

  AF_ = A->values();

  Q_ = A->values();

  R_ = A->values();


  return(0);

}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


int SerialQRDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,

                                                           const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)

{

  // Check that these new vectors are consistent

  TEUCHOS_TEST_FOR_EXCEPTION(B->numCols() != X->numCols(), std::invalid_argument,

                     "SerialQRDenseSolver<T>::setVectors: X and B have different numbers of columns!");

  TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,

                     "SerialQRDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");

  TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,

                     "SerialQRDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");

  TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,

                     "SerialQRDenseSolver<T>::setVectors: B has an invalid stride!");

  TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,

                     "SerialQRDenseSolver<T>::setVectors: X has an invalid stride!");


  resetVectors();

  LHS_ = X;

  RHS_ = B;


  return(0);

}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


int SerialQRDenseSolver<OrdinalType,ScalarType>::factor() {


  if (factored()) return(0);


  // Equilibrate matrix if necessary

  int ierr = 0;

  if (equilibrate_) ierr = equilibrateMatrix();

  if (ierr!=0) return(ierr);


  allocateWORK();

  if ((int)TAU_.size()!=Min_MN_) TAU_.resize( Min_MN_ );


  INFO_ = 0;


  // Factor

#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

  EigenMatrixMap aMap( AF_, M_, N_, EigenOuterStride(LDAF_) );

  EigenScalarArray tau;

  EigenScalarArrayMap tauMap(&TAU_[0],TEUCHOS_MIN(M_,N_));

  qr_.compute(aMap);

  tau = qr_.hCoeffs();

  for (OrdinalType i=0; i<tauMap.innerSize(); i++) {

    tauMap(i) = tau(i);

  }

  EigenMatrix qrMat = qr_.matrixQR();

  for (OrdinalType j=0; j<aMap.outerSize(); j++) {

    for (OrdinalType i=0; i<aMap.innerSize(); i++) {

      aMap(i,j) = qrMat(i,j);

    }

  }

#else

  this->GEQRF(M_, N_, AF_, LDAF_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);

#endif


  factored_ = true;


  return(INFO_);

}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


int SerialQRDenseSolver<OrdinalType,ScalarType>::solve() {


  // Check if the matrix is zero

  if (matrixIsZero_) {

    LHS_->putScalar(ScalarTraits<ScalarType>::zero());

    return(0);

  }


  // Equilibrate RHS if necessary

  int ierr = 0;

  if (equilibrate_) {

    ierr = equilibrateRHS();

  }

  if (ierr != 0) return(ierr);


  TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,

                     "SerialQRDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");

  TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,

                     "SerialQRDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");

  if ( TRANS_ == Teuchos::NO_TRANS ) {

    TEUCHOS_TEST_FOR_EXCEPTION( LHS_->numRows() != N_, std::invalid_argument,

                     "SerialQRDenseSolver<T>::solve: No transpose: must have LHS_->numRows() = N_");

  } else {

    TEUCHOS_TEST_FOR_EXCEPTION( LHS_->numRows() != M_, std::invalid_argument,

                     "SerialQRDenseSolver<T>::solve: Transpose: must have LHS_->numRows() = M_");

  }


  if (shouldEquilibrate() && !equilibratedA_)

    std::cout << "WARNING!  SerialQRDenseSolver<T>::solve: System should be equilibrated!" << std::endl;


  // Matrix must be factored

  if (!factored()) factor();


  TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,

                     std::logic_error, "SerialQRDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");


  TMP_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(M_, RHS_->numCols()) );

  for (OrdinalType j=0; j<RHS_->numCols(); j++) {

    for (OrdinalType i=0; i<RHS_->numRows(); i++) {

      (*TMP_)(i,j) = (*RHS_)(i,j);

    }

  }


  INFO_ = 0;


  // Solve

  if ( TRANS_ == Teuchos::NO_TRANS ) {


    // Solve A lhs = rhs

    this->multiplyQ( Teuchos::CONJ_TRANS, *TMP_ );

    this->solveR( Teuchos::NO_TRANS, *TMP_ );


  } else {


    // Solve A**H lhs = rhs

    this->solveR( Teuchos::CONJ_TRANS, *TMP_ );

    for (OrdinalType j = 0; j < RHS_->numCols(); j++ ) {

      for (OrdinalType i = N_; i < M_; i++ ) {

        (*TMP_)(i, j) = ScalarTraits<ScalarType>::zero();

      }

    }

    this->multiplyQ( Teuchos::NO_TRANS, *TMP_ );


  }

  for (OrdinalType j = 0; j < LHS_->numCols(); j++ ) {

    for (OrdinalType i = 0; i < LHS_->numRows(); i++ ) {

      (*LHS_)(i, j) = (*TMP_)(i,j);

    }

  }


  solved_ = true;


  // Unequilibrate LHS if necessary

  if (equilibrate_) {

    ierr = unequilibrateLHS();

  }

  if (ierr != 0) {

    return ierr;

  }


  return INFO_;

}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


int SerialQRDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling()

{

  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();

  MagnitudeType prec = ScalarTraits<ScalarType>::prec();

  ScalarType sZero = ScalarTraits<ScalarType>::zero();

  ScalarType sOne  = ScalarTraits<ScalarType>::one();

  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);


  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);

  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);


  // Compute maximum absolute value of matrix entries

  OrdinalType i, j;

  MagnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());

  for (j = 0; j < N_; j++) {

    for (i = 0; i < M_; i++) {

      anorm = TEUCHOS_MAX( anorm, ScalarTraits<ScalarType>::magnitude((*Matrix_)(i,j)) );

    }

  }


  ANORM_ = anorm;


  if (ANORM_ > mZero && ANORM_ < smlnum) {

    // Scale matrix norm up to smlnum

    shouldEquilibrate_ = true;

  } else if (ANORM_ > bignum) {

    // Scale matrix norm down to bignum

    shouldEquilibrate_ = true;

  } else if (ANORM_ == mZero) {

    // Matrix all zero. Return zero solution

    matrixIsZero_ = true;

  }


  return(0);

}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


int SerialQRDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()

{

  if (equilibratedA_) return(0);


  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();

  MagnitudeType prec = ScalarTraits<ScalarType>::prec();

  ScalarType sZero = ScalarTraits<ScalarType>::zero();

  ScalarType sOne  = ScalarTraits<ScalarType>::one();

  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);


  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);

  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);

  OrdinalType BW = 0;


  // Compute maximum absolute value of matrix entries

  OrdinalType i, j;

  MagnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());

  for (j = 0; j < N_; j++) {

    for (i = 0; i < M_; i++) {

      anorm = TEUCHOS_MAX( anorm, ScalarTraits<ScalarType>::magnitude((*Matrix_)(i,j)) );

    }

  }


  ANORM_ = anorm;

  int ierr1 = 0;

  if (ANORM_ > mZero && ANORM_ < smlnum) {

    // Scale matrix norm up to smlnum

    this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, smlnum, M_, N_, A_, LDA_, &INFO_);

    equilibrationModeA_ = 1;

  } else if (ANORM_ > bignum) {

    // Scale matrix norm down to bignum

    this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, bignum, M_, N_, A_, LDA_, &INFO_);

    equilibrationModeA_ = 2;

  } else if (ANORM_ == mZero) {

    // Matrix all zero. Return zero solution

    matrixIsZero_ = true;

  }


  equilibratedA_ = true;


  return(ierr1);

}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


int SerialQRDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()

{

  if (equilibratedB_) return(0);


  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();

  MagnitudeType prec = ScalarTraits<ScalarType>::prec();

  ScalarType sZero = ScalarTraits<ScalarType>::zero();

  ScalarType sOne  = ScalarTraits<ScalarType>::one();

  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);


  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);

  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);

  OrdinalType BW = 0;


  // Compute maximum absolute value of rhs entries

  OrdinalType i, j;

  MagnitudeType bnorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());

  for (j = 0; j <RHS_->numCols(); j++) {

    for (i = 0; i < RHS_->numRows(); i++) {

      bnorm = TEUCHOS_MAX( bnorm, ScalarTraits<ScalarType>::magnitude((*RHS_)(i,j)) );

    }

  }


  BNORM_ = bnorm;


  int ierr1 = 0;

  if (BNORM_ > mZero && BNORM_ < smlnum) {

    // Scale RHS norm up to smlnum

    this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, BNORM_, smlnum, RHS_->numRows(), RHS_->numCols(),

                RHS_->values(), RHS_->stride(), &INFO_);

    equilibrationModeB_ = 1;

  } else if (BNORM_ > bignum) {

    // Scale RHS norm down to bignum

    this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, BNORM_, bignum, RHS_->numRows(), RHS_->numCols(),

                RHS_->values(), RHS_->stride(), &INFO_);

    equilibrationModeB_ = 2;

  }


  equilibratedB_ = true;


  return(ierr1);

}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


int SerialQRDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()

{

  if (!equilibratedB_) return(0);


  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();

  MagnitudeType prec = ScalarTraits<ScalarType>::prec();

  ScalarType sZero = ScalarTraits<ScalarType>::zero();

  ScalarType sOne  = ScalarTraits<ScalarType>::one();

  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);

  (void) mZero; // Silence "unused variable" compiler warning.


  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);

  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);

  OrdinalType BW = 0;


  int ierr1 = 0;

  if (equilibrationModeA_ == 1) {

    this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, smlnum, LHS_->numRows(), LHS_->numCols(),

                LHS_->values(), LHS_->stride(), &INFO_);

  } else if (equilibrationModeA_ == 2) {

    this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, bignum, LHS_->numRows(), LHS_->numCols(),

                LHS_->values(), LHS_->stride(), &INFO_);

  }

  if (equilibrationModeB_ == 1) {

    this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, smlnum, BNORM_, LHS_->numRows(), LHS_->numCols(),

                LHS_->values(), LHS_->stride(), &INFO_);

  } else if (equilibrationModeB_ == 2) {

    this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, bignum, BNORM_, LHS_->numRows(), LHS_->numCols(),

                LHS_->values(), LHS_->stride(), &INFO_);

  }


  return(ierr1);

}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


int SerialQRDenseSolver<OrdinalType,ScalarType>::formQ() {


  // Matrix must be factored first

  if (!factored()) factor();


  // Store Q separately from factored matrix

  if (AF_ == Q_) {

    FactorQ_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(*Factor_) );

    Q_ = FactorQ_->values();

    LDQ_ = FactorQ_->stride();

  }


  INFO_ = 0;


  // Form Q

#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

  EigenMatrixMap qMap( Q_, M_, N_, EigenOuterStride(LDQ_) );

  EigenMatrix qMat = qr_.householderQ();

  for (OrdinalType j=0; j<qMap.outerSize(); j++) {

    for (OrdinalType i=0; i<qMap.innerSize(); i++) {

      qMap(i,j) = qMat(i,j);

    }

  }

#else

  this->UNGQR(M_, N_, N_, Q_, LDQ_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);

#endif


  if (INFO_!=0) return(INFO_);


  formedQ_ = true;


  return(INFO_);

}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


int SerialQRDenseSolver<OrdinalType,ScalarType>::formR() {


  // Matrix must be factored first

  if (!factored()) factor();


  // Store R separately from factored matrix

  if (AF_ == R_) {

    FactorR_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(N_, N_) );

    R_ = FactorR_->values();

    LDR_ = FactorR_->stride();

  }


  // Form R

  for (OrdinalType j=0; j<N_; j++) {

    for (OrdinalType i=0; i<=j; i++) {

      (*FactorR_)(i,j) = (*Factor_)(i,j);

    }

  }


  formedR_ = true;


  return(0);

}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


int  SerialQRDenseSolver<OrdinalType, ScalarType>::multiplyQ(ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C)

{


  // Check that the matrices are consistent.

  TEUCHOS_TEST_FOR_EXCEPTION(C.numRows()!=M_, std::invalid_argument,

                     "SerialQRDenseSolver<T>::multiplyQ: C.numRows() != M_!");

  TEUCHOS_TEST_FOR_EXCEPTION(C.values()==0, std::invalid_argument,

                     "SerialQRDenseSolver<T>::multiplyQ: C is an empty SerialDenseMatrix<T>!");


  // Matrix must be factored

  if (!factored()) factor();


  INFO_ = 0;


  // Multiply

#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

  EigenMatrixMap cMap( C.values(), C.numRows(), C.numCols(), EigenOuterStride(C.stride()) );

  if ( transq == Teuchos::NO_TRANS ) {

    // C = Q * C

    cMap = qr_.householderQ() * cMap;

  } else {

    // C = Q**H * C

    cMap = qr_.householderQ().adjoint() * cMap;

    for (OrdinalType j = 0; j < C.numCols(); j++) {

      for (OrdinalType i = N_; i < C.numRows(); i++ ) {

        cMap(i, j) = ScalarTraits<ScalarType>::zero();

      }

    }

  }

#else

  Teuchos::ETransp NO_TRANS = Teuchos::NO_TRANS;

  Teuchos::ETransp TRANS = (Teuchos::ScalarTraits<ScalarType>::isComplex) ? Teuchos::CONJ_TRANS : Teuchos::TRANS;

  Teuchos::ESide SIDE = Teuchos::LEFT_SIDE;


  if ( transq == Teuchos::NO_TRANS ) {


    // C = Q * C

    this->UNMQR(ESideChar[SIDE], ETranspChar[NO_TRANS], M_, C.numCols(), N_, AF_, LDAF_,

                &TAU_[0], C.values(), C.stride(), &WORK_[0], LWORK_, &INFO_);


  } else {


    // C = Q**H * C

    this->UNMQR(ESideChar[SIDE], ETranspChar[TRANS], M_, C.numCols(), N_, AF_, LDAF_,

                &TAU_[0], C.values(), C.stride(), &WORK_[0], LWORK_, &INFO_);


    for (OrdinalType j = 0; j < C.numCols(); j++) {

      for (OrdinalType i = N_; i < C.numRows(); i++ ) {

        C(i, j) = ScalarTraits<ScalarType>::zero();

      }

    }

  }

#endif


  return(INFO_);


}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


int  SerialQRDenseSolver<OrdinalType, ScalarType>::solveR(ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C)

{


  // Check that the matrices are consistent.

  TEUCHOS_TEST_FOR_EXCEPTION(C.numRows()<N_ || C.numRows()>M_, std::invalid_argument,

                     "SerialQRDenseSolver<T>::solveR: must have N_ < C.numRows() < M_!");

  TEUCHOS_TEST_FOR_EXCEPTION(C.values()==0, std::invalid_argument,

                     "SerialQRDenseSolver<T>::solveR: C is an empty SerialDenseMatrix<T>!");


  // Matrix must be factored

  if (!factored()) factor();


  INFO_ = 0;


  // Solve

#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

  EigenMatrixMap cMap( C.values(), N_, C.numCols(), EigenOuterStride(C.stride()) );

  // Check for singularity first like TRTRS

  for (OrdinalType j=0; j<N_; j++) {

    if ((qr_.matrixQR())(j,j) == ScalarTraits<ScalarType>::zero()) {

      INFO_ = j+1;

      return(INFO_);

    }

  }

  if ( transr == Teuchos::NO_TRANS ) {

    // C = R \ C

    qr_.matrixQR().topRows(N_).template triangularView<Eigen::Upper>().solveInPlace(cMap);

  } else {

    // C = R**H \ C

    qr_.matrixQR().topRows(N_).template triangularView<Eigen::Upper>().adjoint().solveInPlace(cMap);

  }

#else

  Teuchos::ETransp NO_TRANS = Teuchos::NO_TRANS;

  Teuchos::ETransp TRANS = (Teuchos::ScalarTraits<ScalarType>::isComplex) ? Teuchos::CONJ_TRANS : Teuchos::TRANS;

  Teuchos::EUplo UPLO = Teuchos::UPPER_TRI;

  Teuchos::EDiag DIAG = Teuchos::NON_UNIT_DIAG;


  ScalarType* RMAT = (formedR_) ? R_ : AF_;

  OrdinalType LDRMAT = (formedR_) ? LDR_ : LDAF_;


  if ( transr == Teuchos::NO_TRANS ) {


    // C = R \ C

    this->TRTRS(EUploChar[UPLO], ETranspChar[NO_TRANS], EDiagChar[DIAG], N_, C.numCols(),

                RMAT, LDRMAT, C.values(), C.stride(), &INFO_);


  } else {


    // C = R**H \ C

    this->TRTRS(EUploChar[UPLO], ETranspChar[TRANS], EDiagChar[DIAG], N_, C.numCols(),

                RMAT, LDRMAT, C.values(), C.stride(), &INFO_);


  }

#endif


  return(INFO_);


}


//=============================================================================


template<typename OrdinalType, typename ScalarType>


void SerialQRDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {


  if (Matrix_!=Teuchos::null) os << "Solver Matrix"          << std::endl << *Matrix_ << std::endl;

  if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;

  if (Q_!=Teuchos::null) os << "Solver Factor Q" << std::endl << *Q_ << std::endl;

  if (LHS_   !=Teuchos::null) os << "Solver LHS"             << std::endl << *LHS_    << std::endl;

  if (RHS_   !=Teuchos::null) os << "Solver RHS"             << std::endl << *RHS_    << std::endl;


}


} // namespace Teuchos


#endif /* _TEUCHOS_SERIALQRDENSESOLVER_HPP_ */

Teuchos_BLAS.hpp
Templated interface class to BLAS routines.

Teuchos_CompObject.hpp
Object for storing data and providing functionality that is common to all computational classes.

Teuchos_ConfigDefs.hpp
Teuchos header file which uses auto-configuration information to include necessary C++ headers.

Teuchos_LAPACK.hpp
Templated interface class to LAPACK routines.

Teuchos_RCP.hpp
Reference-counted pointer class and non-member templated function implementations.

Teuchos_ScalarTraits.hpp
Defines basic traits for the scalar field type.

Teuchos_SerialDenseMatrix.hpp
Templated serial dense matrix class.

Teuchos_SerialDenseSolver.hpp
Templated class for solving dense linear problems.

Teuchos::BLAS::BLAS
BLAS(void)
Default constructor.
Definition Teuchos_BLAS.hpp:222

Teuchos::CompObject::CompObject
CompObject()
Default constructor.
Definition Teuchos_CompObject.cpp:20

Teuchos::LAPACK::TRTRS
void TRTRS(const char &UPLO, const char &TRANS, const char &DIAG, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
Solves a triangular linear system of the form A*X=B or A**T*X=B, where A is a triangular matrix.
Definition Teuchos_LAPACK.cpp:238

Teuchos::LAPACK::UNGQR
void UNGQR(const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, ScalarType *A, const OrdinalType &lda, const ScalarType *TAU, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
Compute explicit QR factor from QR factorization (GEQRF) (complex case).
Definition Teuchos_LAPACK.cpp:358

Teuchos::LAPACK::UNMQR
void UNMQR(const char &SIDE, const char &TRANS, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const ScalarType *A, const OrdinalType &lda, const ScalarType *TAU, ScalarType *C, const OrdinalType &ldc, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
Apply Householder reflectors (complex case).
Definition Teuchos_LAPACK.cpp:339

Teuchos::LAPACK::LASCL
void LASCL(const char &TYPE, const OrdinalType &kl, const OrdinalType &ku, const MagnitudeType cfrom, const MagnitudeType cto, const OrdinalType &m, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
Multiplies the m by n matrix A by the real scalar cto/cfrom.
Definition Teuchos_LAPACK.cpp:157

Teuchos::LAPACK::LAPACK
LAPACK(void)
Default Constructor.
Definition Teuchos_LAPACK.hpp:74

Teuchos::LAPACK::GEQRF
void GEQRF(const OrdinalType &m, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, ScalarType *TAU, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
Computes a QR factorization of a general m by n matrix A.
Definition Teuchos_LAPACK.cpp:141

Teuchos::Object::Object
Object(int tracebackModeIn=-1)
Default Constructor.
Definition Teuchos_Object.cpp:20

Teuchos::RCP
Smart reference counting pointer class for automatic garbage collection.
Definition Teuchos_RCPDecl.hpp:397

Teuchos::SerialDenseMatrix
This class creates and provides basic support for dense rectangular matrix of templated type.
Definition Teuchos_SerialDenseMatrix.hpp:38

Teuchos::SerialDenseMatrix::stride
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
Definition Teuchos_SerialDenseMatrix.hpp:329

Teuchos::SerialDenseMatrix::numRows
OrdinalType numRows() const
Returns the row dimension of this matrix.
Definition Teuchos_SerialDenseMatrix.hpp:323

Teuchos::SerialDenseMatrix::values
ScalarType * values() const
Data array access method.
Definition Teuchos_SerialDenseMatrix.hpp:228

Teuchos::SerialDenseMatrix::numCols
OrdinalType numCols() const
Returns the column dimension of this matrix.
Definition Teuchos_SerialDenseMatrix.hpp:326

Teuchos::SerialQRDenseSolver::factorWithEquilibration
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
Definition Teuchos_SerialQRDenseSolver.hpp:154

Teuchos::SerialQRDenseSolver::solveR
int solveR(ETransp transr, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Solve input matrix on the left with the upper triangular matrix R or its adjoint.
Definition Teuchos_SerialQRDenseSolver.hpp:922

Teuchos::SerialQRDenseSolver::solved
bool solved()
Returns true if the current set of vectors has been solved.
Definition Teuchos_SerialQRDenseSolver.hpp:250

Teuchos::SerialQRDenseSolver::equilibratedB
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
Definition Teuchos_SerialQRDenseSolver.hpp:244

Teuchos::SerialQRDenseSolver::getLHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
Definition Teuchos_SerialQRDenseSolver.hpp:276

Teuchos::SerialQRDenseSolver::formedR
bool formedR()
Returns true if R has been formed explicitly.
Definition Teuchos_SerialQRDenseSolver.hpp:256

Teuchos::SerialQRDenseSolver::formR
int formR()
Explicitly forms the upper triangular matrix R.
Definition Teuchos_SerialQRDenseSolver.hpp:834

Teuchos::SerialQRDenseSolver::shouldEquilibrate
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
Definition Teuchos_SerialQRDenseSolver.hpp:247

Teuchos::SerialQRDenseSolver::setVectors
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
Sets the pointers for left and right hand side vector(s).
Definition Teuchos_SerialQRDenseSolver.hpp:477

Teuchos::SerialQRDenseSolver::getMatrix
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
Definition Teuchos_SerialQRDenseSolver.hpp:264

Teuchos::SerialQRDenseSolver::solveWithTranspose
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the adjoint of this matrix,...
Definition Teuchos_SerialQRDenseSolver.hpp:157

Teuchos::SerialQRDenseSolver::transpose
bool transpose()
Returns true if adjoint of this matrix has and will be used.
Definition Teuchos_SerialQRDenseSolver.hpp:235

Teuchos::SerialQRDenseSolver::multiplyQ
int multiplyQ(ETransp transq, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Left multiply the input matrix by the unitary matrix Q or its adjoint.
Definition Teuchos_SerialQRDenseSolver.hpp:861

Teuchos::SerialQRDenseSolver::getQ
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
Definition Teuchos_SerialQRDenseSolver.hpp:270

Teuchos::SerialQRDenseSolver::ANORM
MagnitudeType ANORM() const
Returns the absolute value of the largest element of this matrix (returns -1 if not yet computed).
Definition Teuchos_SerialQRDenseSolver.hpp:291

Teuchos::SerialQRDenseSolver::numCols
OrdinalType numCols() const
Returns column dimension of system.
Definition Teuchos_SerialQRDenseSolver.hpp:285

Teuchos::SerialQRDenseSolver::formedQ
bool formedQ()
Returns true if Q has been formed explicitly.
Definition Teuchos_SerialQRDenseSolver.hpp:253

Teuchos::SerialQRDenseSolver::tau
std::vector< ScalarType > tau() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
Definition Teuchos_SerialQRDenseSolver.hpp:288

Teuchos::SerialQRDenseSolver::numRows
OrdinalType numRows() const
Returns row dimension of system.
Definition Teuchos_SerialQRDenseSolver.hpp:282

Teuchos::SerialQRDenseSolver::~SerialQRDenseSolver
virtual ~SerialQRDenseSolver()
SerialQRDenseSolver destructor.
Definition Teuchos_SerialQRDenseSolver.hpp:404

Teuchos::SerialQRDenseSolver::factor
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
Definition Teuchos_SerialQRDenseSolver.hpp:501

Teuchos::SerialQRDenseSolver::getFactoredMatrix
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
Definition Teuchos_SerialQRDenseSolver.hpp:267

Teuchos::SerialQRDenseSolver::solve
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
Definition Teuchos_SerialQRDenseSolver.hpp:543

Teuchos::SerialQRDenseSolver::SerialQRDenseSolver
SerialQRDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
Definition Teuchos_SerialQRDenseSolver.hpp:368

Teuchos::SerialQRDenseSolver::solveWithTransposeFlag
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
Definition Teuchos_SerialQRDenseSolver.hpp:160

Teuchos::SerialQRDenseSolver::factored
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
Definition Teuchos_SerialQRDenseSolver.hpp:238

Teuchos::SerialQRDenseSolver::Print
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
Definition Teuchos_SerialQRDenseSolver.hpp:984

Teuchos::SerialQRDenseSolver::getR
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getR() const
Returns pointer to R (assuming factorization has been performed).
Definition Teuchos_SerialQRDenseSolver.hpp:273

Teuchos::SerialQRDenseSolver::setMatrix
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
Definition Teuchos_SerialQRDenseSolver.hpp:450

Teuchos::SerialQRDenseSolver::formQ
int formQ()
Explicitly forms the unitary matrix Q.
Definition Teuchos_SerialQRDenseSolver.hpp:797

Teuchos::SerialQRDenseSolver::equilibrateMatrix
int equilibrateMatrix()
Equilibrates the this matrix.
Definition Teuchos_SerialQRDenseSolver.hpp:668

Teuchos::SerialQRDenseSolver::getRHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
Definition Teuchos_SerialQRDenseSolver.hpp:279

Teuchos::SerialQRDenseSolver::equilibrateRHS
int equilibrateRHS()
Equilibrates the current RHS.
Definition Teuchos_SerialQRDenseSolver.hpp:714

Teuchos::SerialQRDenseSolver::equilibratedA
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
Definition Teuchos_SerialQRDenseSolver.hpp:241

Teuchos::SerialQRDenseSolver::computeEquilibrateScaling
int computeEquilibrateScaling()
Determines if this matrix should be scaled.
Definition Teuchos_SerialQRDenseSolver.hpp:629

Teuchos::SerialQRDenseSolver::unequilibrateLHS
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
Definition Teuchos_SerialQRDenseSolver.hpp:760

TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
Definition Teuchos_TestForException.hpp:146

Teuchos
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos,...

Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
Definition Teuchos_RCPDecl.hpp:1234

Teuchos::ETransp
ETransp
Definition Teuchos_BLAS_types.hpp:59

Teuchos::CONJ_TRANS
@ CONJ_TRANS
Definition Teuchos_BLAS_types.hpp:62

Teuchos::NO_TRANS
@ NO_TRANS
Definition Teuchos_BLAS_types.hpp:60

Teuchos::TRANS
@ TRANS
Definition Teuchos_BLAS_types.hpp:61

Teuchos::ESide
ESide
Definition Teuchos_BLAS_types.hpp:54

Teuchos::LEFT_SIDE
@ LEFT_SIDE
Definition Teuchos_BLAS_types.hpp:55

Teuchos::EDiag
EDiag
Definition Teuchos_BLAS_types.hpp:71

Teuchos::NON_UNIT_DIAG
@ NON_UNIT_DIAG
Definition Teuchos_BLAS_types.hpp:73

Teuchos::EUplo
EUplo
Definition Teuchos_BLAS_types.hpp:65

Teuchos::UPPER_TRI
@ UPPER_TRI
Definition Teuchos_BLAS_types.hpp:66

Teuchos::FULL
@ FULL
Definition Teuchos_BLAS_types.hpp:77

Teuchos::ScalarTraits
This structure defines some basic traits for a scalar field type.
Definition Teuchos_ScalarTraitsDecl.hpp:59

Teuchos::ScalarTraits::magnitude
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
Definition Teuchos_ScalarTraitsDecl.hpp:100

Teuchos::ScalarTraits::one
static T one()
Returns representation of one for this scalar type.
Definition Teuchos_ScalarTraitsDecl.hpp:104

Teuchos::ScalarTraits::magnitudeType
T magnitudeType
Mandatory typedef for result of magnitude.
Definition Teuchos_ScalarTraitsDecl.hpp:61

Teuchos::ScalarTraits::zero
static T zero()
Returns representation of zero for this scalar type.
Definition Teuchos_ScalarTraitsDecl.hpp:102

Teuchos::ScalarTraits::isComplex
static const bool isComplex
Determines if scalar type is std::complex.
Definition Teuchos_ScalarTraitsDecl.hpp:69

Teuchos::ScalarTraits::sfmin
static magnitudeType sfmin()
Returns safe minimum (sfmin), such that 1/sfmin does not overflow.
Definition Teuchos_ScalarTraitsDecl.hpp:82

Teuchos::ScalarTraits::prec
static magnitudeType prec()
Returns eps*base.
Definition Teuchos_ScalarTraitsDecl.hpp:86