trilinos/teuchos/Teuchos__SerialSpdDenseSolver_8hpp_source.html

// @HEADER

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

//                    Teuchos: Common Tools Package

//

// Copyright 2004 NTESS and the Teuchos contributors.

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

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

// @HEADER


#ifndef TEUCHOS_SERIALSPDDENSESOLVER_H

#define TEUCHOS_SERIALSPDDENSESOLVER_H


#include "Teuchos_CompObject.hpp"

#include "Teuchos_BLAS.hpp"

#include "Teuchos_LAPACK.hpp"

#include "Teuchos_RCP.hpp"

#include "Teuchos_ConfigDefs.hpp"

#include "Teuchos_SerialSymDenseMatrix.hpp"

#include "Teuchos_SerialDenseSolver.hpp"

#include "Teuchos_SerialDenseMatrix.hpp"


#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

#include "Eigen/Dense"

#endif


namespace Teuchos {


  template<typename OrdinalType, typename ScalarType>


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

                            public LAPACK<OrdinalType, ScalarType>

  {

  public:


    typedef typename Teuchos::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


    SerialSpdDenseSolver();


    virtual ~SerialSpdDenseSolver();


    int setMatrix(const RCP<SerialSymDenseMatrix<OrdinalType,ScalarType> >& A_in);


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

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


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


    void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}


    void estimateSolutionErrors(bool flag);


    int factor();


    int solve();


    int invert();


    int computeEquilibrateScaling();


    int equilibrateMatrix();


    int equilibrateRHS();


    int applyRefinement();


    int unequilibrateLHS();


    int reciprocalConditionEstimate(MagnitudeType & Value);


    bool transpose() {return(transpose_);}


    bool factored() {return(factored_);}


    bool equilibratedA() {return(equilibratedA_);}


    bool equilibratedB() {return(equilibratedB_);}


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


    bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}


    bool inverted() {return(inverted_);}


    bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}


    bool solved() {return(solved_);}


    bool solutionRefined() {return(solutionRefined_);}


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


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


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


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


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


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


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


    MagnitudeType RCOND()  const {return(RCOND_);}


    MagnitudeType SCOND() {return(SCOND_);};


    MagnitudeType AMAX()  const {return(AMAX_);}


    std::vector<MagnitudeType> FERR()  const {return(FERR_);}


    std::vector<MagnitudeType> BERR()  const {return(BERR_);}


    std::vector<MagnitudeType> R()  const {return(R_);}


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


  protected:


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

    void allocateIWORK() { IWORK_.resize( numRowCols_ ); return;}

    void resetMatrix();

    void resetVectors();


    bool equilibrate_;

    bool shouldEquilibrate_;

    bool equilibratedA_;

    bool equilibratedB_;

    bool transpose_;

    bool factored_;

    bool estimateSolutionErrors_;

    bool solutionErrorsEstimated_;

    bool solved_;

    bool inverted_;

    bool reciprocalConditionEstimated_;

    bool refineSolution_;

    bool solutionRefined_;


    OrdinalType numRowCols_;

    OrdinalType LDA_;

    OrdinalType LDAF_;

    OrdinalType INFO_;

    OrdinalType LWORK_;


    std::vector<int> IWORK_;


    MagnitudeType ANORM_;

    MagnitudeType RCOND_;

    MagnitudeType SCOND_;

    MagnitudeType AMAX_;


    RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > Matrix_;

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;

    RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > Factor_;


    ScalarType * A_;

    ScalarType * AF_;

    std::vector<MagnitudeType> FERR_;

    std::vector<MagnitudeType> BERR_;

    std::vector<ScalarType> WORK_;

    std::vector<MagnitudeType> R_;

#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

    Eigen::LLT<EigenMatrix,Eigen::Upper> lltu_;

    Eigen::LLT<EigenMatrix,Eigen::Lower> lltl_;

#endif


  private:

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


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

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


  };


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

  using namespace details;


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


template<typename OrdinalType, typename ScalarType>


SerialSpdDenseSolver<OrdinalType,ScalarType>::SerialSpdDenseSolver()

  : CompObject(),

    equilibrate_(false),

    shouldEquilibrate_(false),

    equilibratedA_(false),

    equilibratedB_(false),

    transpose_(false),

    factored_(false),

    estimateSolutionErrors_(false),

    solutionErrorsEstimated_(false),

    solved_(false),

    inverted_(false),

    reciprocalConditionEstimated_(false),

    refineSolution_(false),

    solutionRefined_(false),

    numRowCols_(0),

    LDA_(0),

    LDAF_(0),

    INFO_(0),

    LWORK_(0),

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

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

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

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

    A_(0),

    AF_(0)

{

  resetMatrix();

}


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


template<typename OrdinalType, typename ScalarType>


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

{}


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


template<typename OrdinalType, typename ScalarType>

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

{

  LHS_ = Teuchos::null;

  RHS_ = Teuchos::null;

  reciprocalConditionEstimated_ = false;

  solutionRefined_ = false;

  solved_ = false;

  solutionErrorsEstimated_ = false;

  equilibratedB_ = false;

}

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


template<typename OrdinalType, typename ScalarType>

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

{

  resetVectors();

  equilibratedA_ = false;

  factored_ = false;

  inverted_ = false;

  numRowCols_ = 0;

  LDA_ = 0;

  LDAF_ = 0;

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

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

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

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

  A_ = 0;

  AF_ = 0;

  INFO_ = 0;

  LWORK_ = 0;

  R_.resize(0);

}

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


template<typename OrdinalType, typename ScalarType>


int SerialSpdDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialSymDenseMatrix<OrdinalType,ScalarType> >& A)

{

  resetMatrix();

  Matrix_ = A;

  Factor_ = A;

  numRowCols_ = A->numRows();

  LDA_ = A->stride();

  LDAF_ = LDA_;

  A_ = A->values();

  AF_ = A->values();

  return(0);

}


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


template<typename OrdinalType, typename ScalarType>


int SerialSpdDenseSolver<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->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,

                     "SerialSpdDenseSolver<T>::setVectors: X and B are not the same size!");

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

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

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

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

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

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

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

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


  resetVectors();

  LHS_ = X;

  RHS_ = B;

  return(0);

}


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


template<typename OrdinalType, typename ScalarType>


void SerialSpdDenseSolver<OrdinalType,ScalarType>::estimateSolutionErrors(bool flag)

{

  estimateSolutionErrors_ = flag;


  // If the errors are estimated, this implies that the solution must be refined

  refineSolution_ = refineSolution_ || flag;

}


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


template<typename OrdinalType, typename ScalarType>


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


  if (factored()) return(0); // Already factored


  TEUCHOS_TEST_FOR_EXCEPTION(inverted(), std::logic_error,

                     "SerialSpdDenseSolver<T>::factor: Cannot factor an inverted matrix!");


  ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A


  // If we want to refine the solution, then the factor must

  // be stored separatedly from the original matrix


  if (A_ == AF_)

    if (refineSolution_ ) {

      Factor_ = rcp( new SerialSymDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );

      AF_ = Factor_->values();

      LDAF_ = Factor_->stride();

    }


  int ierr = 0;

  if (equilibrate_) ierr = equilibrateMatrix();


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


  INFO_ = 0;


#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

  EigenMatrixMap aMap( AF_, numRowCols_, numRowCols_, EigenOuterStride(LDAF_) );


  if (Matrix_->upper())

    lltu_.compute(aMap);

  else

    lltl_.compute(aMap);


#else

  this->POTRF(Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);

#endif


  factored_ = true;


  return(INFO_);

}


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


template<typename OrdinalType, typename ScalarType>


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


  // We will call one of four routines depending on what services the user wants and

  // whether or not the matrix has been inverted or factored already.

  //

  // If the matrix has been inverted, use DGEMM to compute solution.

  // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will

  // call the X interface.

  // Otherwise, if the matrix is already factored we will call the TRS interface.

  // Otherwise, if the matrix is unfactored we will call the SV interface.


  int ierr = 0;

  if (equilibrate_) {

    ierr = equilibrateRHS();

  }

  if (ierr != 0) return(ierr);  // Can't equilibrate B, so return.


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

                     "SerialSpdDenseSolver<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,

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


  if (inverted()) {


    TEUCHOS_TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument,

                        "SerialSpdDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");

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

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


    INFO_ = 0;

    this->GEMM(Teuchos::NO_TRANS, Teuchos::NO_TRANS, numRowCols_, RHS_->numCols(),

               numRowCols_, 1.0, AF_, LDAF_, RHS_->values(), RHS_->stride(), 0.0,

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

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

    solved_ = true;

  }

  else {


    if (!factored()) factor(); // Matrix must be factored


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

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


    if (RHS_->values()!=LHS_->values()) {

       (*LHS_) = (*RHS_); // Copy B to X if needed

    }

    INFO_ = 0;


#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

    EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );

    EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );


  if (Matrix_->upper())

    lhsMap = lltu_.solve(rhsMap);

  else

    lhsMap = lltl_.solve(rhsMap);


#else

    this->POTRS(Matrix_->UPLO(), numRowCols_, RHS_->numCols(), AF_, LDAF_, LHS_->values(), LHS_->stride(), &INFO_);

#endif


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

    solved_ = true;


  }


  // Warn users that the system should be equilibrated, but it wasn't.

  if (shouldEquilibrate() && !equilibratedA_)

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


  int ierr1=0;

  if (refineSolution_ && !inverted()) ierr1 = applyRefinement();

  if (ierr1!=0)

    return(ierr1);


  if (equilibrate_) ierr1 = unequilibrateLHS();

  return(ierr1);

}


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


template<typename OrdinalType, typename ScalarType>


int SerialSpdDenseSolver<OrdinalType,ScalarType>::applyRefinement()

{

  TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,

                     "SerialSpdDenseSolver<T>::applyRefinement: Must have an existing solution!");

  TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,

                     "SerialSpdDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");


#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

  // Implement templated PORFS or use Eigen.

  return (-1);

#else

  OrdinalType NRHS = RHS_->numCols();

  FERR_.resize( NRHS );

  BERR_.resize( NRHS );

  allocateWORK();

  allocateIWORK();


  INFO_ = 0;

  std::vector<typename details::lapack_traits<ScalarType>::iwork_type> PORFS_WORK( numRowCols_ );

  this->PORFS(Matrix_->UPLO(), numRowCols_, NRHS, A_, LDA_, AF_, LDAF_,

              RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),

              &FERR_[0], &BERR_[0], &WORK_[0], &PORFS_WORK[0], &INFO_);


  solutionErrorsEstimated_ = true;

  reciprocalConditionEstimated_ = true;

  solutionRefined_ = true;


  return(INFO_);

#endif

}


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


template<typename OrdinalType, typename ScalarType>


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

{

  if (R_.size()!=0) return(0); // Already computed


  R_.resize( numRowCols_ );


  INFO_ = 0;

  this->POEQU (numRowCols_, AF_, LDAF_, &R_[0], &SCOND_, &AMAX_, &INFO_);

  if (SCOND_<0.1*ScalarTraits<MagnitudeType>::one() ||

      AMAX_ < ScalarTraits<ScalarType>::rmin() ||

      AMAX_ > ScalarTraits<ScalarType>::rmax())

    shouldEquilibrate_ = true;


  return(INFO_);

}


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


template<typename OrdinalType, typename ScalarType>


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

{

  OrdinalType i, j;

  int ierr = 0;


  if (equilibratedA_) return(0); // Already done.

  if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R if needed.

  if (ierr!=0) return(ierr);     // If return value is not zero, then can't equilibrate.

  if (Matrix_->upper()) {

    if (A_==AF_) {

      ScalarType * ptr;

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

        ptr = A_ + j*LDA_;

        ScalarType s1 = R_[j];

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

          *ptr = *ptr*s1*R_[i];

          ptr++;

        }

      }

    }

    else {

      ScalarType * ptr;

      ScalarType * ptr1;

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

        ptr = A_ + j*LDA_;

        ptr1 = AF_ + j*LDAF_;

        ScalarType s1 = R_[j];

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

          *ptr = *ptr*s1*R_[i];

          ptr++;

          *ptr1 = *ptr1*s1*R_[i];

          ptr1++;

        }

      }

    }

  }

  else {

    if (A_==AF_) {

      ScalarType * ptr;

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

        ptr = A_ + j + j*LDA_;

        ScalarType s1 = R_[j];

        for (i=j; i<numRowCols_; i++) {

          *ptr = *ptr*s1*R_[i];

          ptr++;

        }

      }

    }

    else {

      ScalarType * ptr;

      ScalarType * ptr1;

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

        ptr = A_ + j + j*LDA_;

        ptr1 = AF_ + j + j*LDAF_;

        ScalarType s1 = R_[j];

        for (i=j; i<numRowCols_; i++) {

          *ptr = *ptr*s1*R_[i];

          ptr++;

          *ptr1 = *ptr1*s1*R_[i];

          ptr1++;

        }

      }

    }

  }


  equilibratedA_ = true;


  return(0);

}


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


template<typename OrdinalType, typename ScalarType>


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

{

  OrdinalType i, j;

  int ierr = 0;


  if (equilibratedB_) return(0); // Already done.

  if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R if needed.

  if (ierr!=0) return(ierr);     // Can't count on R being computed.


  OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();

  ScalarType * B = RHS_->values();

  ScalarType * ptr;

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

    ptr = B + j*LDB;

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

      *ptr = *ptr*R_[i];

      ptr++;

    }

  }


  equilibratedB_ = true;


  return(0);

}


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


template<typename OrdinalType, typename ScalarType>


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

{

  OrdinalType i, j;


  if (!equilibratedB_) return(0); // Nothing to do


  OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();

  ScalarType * X = LHS_->values();

  ScalarType * ptr;

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

    ptr = X + j*LDX;

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

      *ptr = *ptr*R_[i];

      ptr++;

    }

  }


  return(0);

}


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


template<typename OrdinalType, typename ScalarType>


int SerialSpdDenseSolver<OrdinalType,ScalarType>::invert()

{

#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

  // Implement templated inverse or use Eigen.

  return (-1);

#else

  if (!factored()) factor(); // Need matrix factored.


  INFO_ = 0;

  this->POTRI( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);


  // Copy lower/upper triangle to upper/lower triangle to make full inverse.

  if (Matrix_->upper())

    Matrix_->setLower();

  else

    Matrix_->setUpper();


  inverted_ = true;

  factored_ = false;


  return(INFO_);

#endif

}


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


template<typename OrdinalType, typename ScalarType>


int SerialSpdDenseSolver<OrdinalType,ScalarType>::reciprocalConditionEstimate(MagnitudeType & Value)

{

#ifdef HAVE_TEUCHOSNUMERICS_EIGEN

  // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.

  return (-1);

#else

  if (reciprocalConditionEstimated()) {

    Value = RCOND_;

    return(0); // Already computed, just return it.

  }


  if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();


  int ierr = 0;

  if (!factored()) ierr = factor(); // Need matrix factored.

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


  allocateWORK();

  allocateIWORK();


  // We will assume a one-norm condition number

  INFO_ = 0;

  std::vector<typename details::lapack_traits<ScalarType>::iwork_type> POCON_WORK( numRowCols_ );

  this->POCON( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &POCON_WORK[0], &INFO_);

  reciprocalConditionEstimated_ = true;

  Value = RCOND_;


  return(INFO_);

#endif

}


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


template<typename OrdinalType, typename ScalarType>


void SerialSpdDenseSolver<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 (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_SERIALSPDDENSESOLVER_H */

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_SerialDenseMatrix.hpp
Templated serial dense matrix class.

Teuchos_SerialDenseSolver.hpp
Templated class for solving dense linear problems.

Teuchos_SerialSymDenseMatrix.hpp
Templated serial, dense, symmetric matrix class.

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

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

Teuchos::DefaultBLASImpl::GEMM
void GEMM(ETransp transa, ETransp transb, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
General matrix-matrix multiply.
Definition Teuchos_BLAS.hpp:1109

Teuchos::LAPACK::PORFS
void PORFS(const char &UPLO, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *A, const OrdinalType &lda, const ScalarType *AF, const OrdinalType &ldaf, const ScalarType *B, const OrdinalType &ldb, ScalarType *X, const OrdinalType &ldx, ScalarType *FERR, ScalarType *BERR, ScalarType *WORK, OrdinalType *IWORK, OrdinalType *info) const
Improves the computed solution to a system of linear equations when the coefficient matrix is symmetr...
Definition Teuchos_LAPACK.cpp:117

Teuchos::LAPACK::POCON
void POCON(const char &UPLO, const OrdinalType &n, const ScalarType *A, const OrdinalType &lda, const ScalarType &anorm, ScalarType *rcond, ScalarType *WORK, OrdinalType *IWORK, OrdinalType *info) const
Estimates the reciprocal of the condition number (1-norm) of a real symmetric positive definite matri...
Definition Teuchos_LAPACK.cpp:106

Teuchos::LAPACK::POTRS
void POTRS(const char &UPLO, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
Solves a system of linear equations A*X=B, where A is a symmetric positive definite matrix factored b...
Definition Teuchos_LAPACK.cpp:98

Teuchos::LAPACK::POEQU
void POEQU(const OrdinalType &n, const ScalarType *A, const OrdinalType &lda, MagnitudeType *S, MagnitudeType *scond, MagnitudeType *amax, OrdinalType *info) const
Computes row and column scalings intended to equilibrate a symmetric positive definite matrix A and r...
Definition Teuchos_LAPACK.cpp:114

Teuchos::LAPACK::POTRI
void POTRI(const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
Computes the inverse of a real symmetric positive definite matrix A using the Cholesky factorization ...
Definition Teuchos_LAPACK.cpp:102

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

Teuchos::LAPACK::POTRF
void POTRF(const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
Computes Cholesky factorization of a real symmetric positive definite matrix A.
Definition Teuchos_LAPACK.cpp:94

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

Teuchos::Ptr::ptr
Ptr< T > ptr(T *p)
Create a pointer to an object from a raw pointer.
Definition Teuchos_PtrDecl.hpp:296

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::SerialSpdDenseSolver::SCOND
MagnitudeType SCOND()
Ratio of smallest to largest equilibration scale factors for the this matrix (returns -1 if not yet c...
Definition Teuchos_SerialSpdDenseSolver.hpp:292

Teuchos::SerialSpdDenseSolver::setMatrix
int setMatrix(const RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > &A_in)
Sets the pointers for coefficient matrix.
Definition Teuchos_SerialSpdDenseSolver.hpp:451

Teuchos::SerialSpdDenseSolver::numCols
OrdinalType numCols() const
Returns column dimension of system.
Definition Teuchos_SerialSpdDenseSolver.hpp:281

Teuchos::SerialSpdDenseSolver::FERR
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
Definition Teuchos_SerialSpdDenseSolver.hpp:298

Teuchos::SerialSpdDenseSolver::applyRefinement
int applyRefinement()
Apply Iterative Refinement.
Definition Teuchos_SerialSpdDenseSolver.hpp:627

Teuchos::SerialSpdDenseSolver::transpose
bool transpose()
Returns true if transpose of this matrix has and will be used.
Definition Teuchos_SerialSpdDenseSolver.hpp:232

Teuchos::SerialSpdDenseSolver::solve
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
Definition Teuchos_SerialSpdDenseSolver.hpp:546

Teuchos::SerialSpdDenseSolver::R
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
Definition Teuchos_SerialSpdDenseSolver.hpp:304

Teuchos::SerialSpdDenseSolver::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_SerialSpdDenseSolver.hpp:466

Teuchos::SerialSpdDenseSolver::factor
int factor()
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF.
Definition Teuchos_SerialSpdDenseSolver.hpp:499

Teuchos::SerialSpdDenseSolver::getRHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
Definition Teuchos_SerialSpdDenseSolver.hpp:275

Teuchos::SerialSpdDenseSolver::estimateSolutionErrors
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
Definition Teuchos_SerialSpdDenseSolver.hpp:489

Teuchos::SerialSpdDenseSolver::solutionRefined
bool solutionRefined()
Returns true if the current set of vectors has been refined.
Definition Teuchos_SerialSpdDenseSolver.hpp:259

Teuchos::SerialSpdDenseSolver::~SerialSpdDenseSolver
virtual ~SerialSpdDenseSolver()
SerialSpdDenseSolver destructor.
Definition Teuchos_SerialSpdDenseSolver.hpp:410

Teuchos::SerialSpdDenseSolver::getMatrix
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
Definition Teuchos_SerialSpdDenseSolver.hpp:266

Teuchos::SerialSpdDenseSolver::solved
bool solved()
Returns true if the current set of vectors has been solved.
Definition Teuchos_SerialSpdDenseSolver.hpp:256

Teuchos::SerialSpdDenseSolver::reciprocalConditionEstimated
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
Definition Teuchos_SerialSpdDenseSolver.hpp:253

Teuchos::SerialSpdDenseSolver::getFactoredMatrix
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
Definition Teuchos_SerialSpdDenseSolver.hpp:269

Teuchos::SerialSpdDenseSolver::solveToRefinedSolution
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement.
Definition Teuchos_SerialSpdDenseSolver.hpp:156

Teuchos::SerialSpdDenseSolver::numRows
OrdinalType numRows() const
Returns row dimension of system.
Definition Teuchos_SerialSpdDenseSolver.hpp:278

Teuchos::SerialSpdDenseSolver::ANORM
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
Definition Teuchos_SerialSpdDenseSolver.hpp:284

Teuchos::SerialSpdDenseSolver::solutionErrorsEstimated
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
Definition Teuchos_SerialSpdDenseSolver.hpp:247

Teuchos::SerialSpdDenseSolver::AMAX
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed).
Definition Teuchos_SerialSpdDenseSolver.hpp:295

Teuchos::SerialSpdDenseSolver::factored
bool factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
Definition Teuchos_SerialSpdDenseSolver.hpp:235

Teuchos::SerialSpdDenseSolver::Print
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
Definition Teuchos_SerialSpdDenseSolver.hpp:865

Teuchos::SerialSpdDenseSolver::SerialSpdDenseSolver
SerialSpdDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
Definition Teuchos_SerialSpdDenseSolver.hpp:378

Teuchos::SerialSpdDenseSolver::BERR
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
Definition Teuchos_SerialSpdDenseSolver.hpp:301

Teuchos::SerialSpdDenseSolver::unequilibrateLHS
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
Definition Teuchos_SerialSpdDenseSolver.hpp:782

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

Teuchos::SerialSpdDenseSolver::equilibrateMatrix
int equilibrateMatrix()
Equilibrates the this matrix.
Definition Teuchos_SerialSpdDenseSolver.hpp:681

Teuchos::SerialSpdDenseSolver::invert
int invert()
Inverts the this matrix.
Definition Teuchos_SerialSpdDenseSolver.hpp:805

Teuchos::SerialSpdDenseSolver::reciprocalConditionEstimate
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
Definition Teuchos_SerialSpdDenseSolver.hpp:832

Teuchos::SerialSpdDenseSolver::equilibratedB
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
Definition Teuchos_SerialSpdDenseSolver.hpp:241

Teuchos::SerialSpdDenseSolver::getLHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
Definition Teuchos_SerialSpdDenseSolver.hpp:272

Teuchos::SerialSpdDenseSolver::inverted
bool inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
Definition Teuchos_SerialSpdDenseSolver.hpp:250

Teuchos::SerialSpdDenseSolver::computeEquilibrateScaling
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the this matrix.
Definition Teuchos_SerialSpdDenseSolver.hpp:662

Teuchos::SerialSpdDenseSolver::RCOND
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed).
Definition Teuchos_SerialSpdDenseSolver.hpp:287

Teuchos::SerialSpdDenseSolver::equilibrateRHS
int equilibrateRHS()
Equilibrates the current RHS.
Definition Teuchos_SerialSpdDenseSolver.hpp:754

Teuchos::SerialSpdDenseSolver::shouldEquilibrate
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
Definition Teuchos_SerialSpdDenseSolver.hpp:244

Teuchos::SerialSpdDenseSolver::equilibratedA
bool equilibratedA()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
Definition Teuchos_SerialSpdDenseSolver.hpp:238

Teuchos::SerialSymDenseMatrix
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
Definition Teuchos_SerialSymDenseMatrix.hpp:92

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::NO_TRANS
@ NO_TRANS
Definition Teuchos_BLAS_types.hpp:60

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

Teuchos::ScalarTraits::rmax
static magnitudeType rmax()
Overflow theshold - (base^emax)*(1-eps).
Definition Teuchos_ScalarTraitsDecl.hpp:98

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::rmin
static magnitudeType rmin()
Returns the underflow threshold - base^(emin-1).
Definition Teuchos_ScalarTraitsDecl.hpp:94