Teuchos_SerialSpdDenseSolver.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                    Teuchos: Common Tools Package
//                 Copyright (2004) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
// @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_SerialDenseMatrix.hpp"

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;
    


    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<ScalarType> FERR()  const {return(FERR_);}
    
    std::vector<ScalarType> BERR()  const {return(BERR_);}
    
    std::vector<ScalarType> 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<ScalarType> FERR_;
    std::vector<ScalarType> BERR_;
    std::vector<ScalarType> WORK_;
    std::vector<ScalarType> R_;
    
  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);
    
  };  

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

template<typename OrdinalType, typename ScalarType>
SerialSpdDenseSolver<OrdinalType,ScalarType>::SerialSpdDenseSolver()
  : CompObject(),
    Object("Teuchos::SerialSpdDenseSolver"),
    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<ScalarType>::zero()),
    RCOND_(ScalarTraits<ScalarType>::zero()),
    SCOND_(ScalarTraits<ScalarType>::zero()),
    AMAX_(ScalarTraits<ScalarType>::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.
  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!");
  TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
         "SerialSpdDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
  TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
         "SerialSpdDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
  TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
         "SerialSpdDenseSolver<T>::setVectors: B has an invalid stride!");
  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

  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;
  this->POTRF(Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
  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();
    equilibratedB_ = true;
  }
  if (ierr != 0) return(ierr);  // Can't equilibrate B, so return.

  TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) , 
                     std::logic_error, "SerialSpdDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
  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!");
  TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument, 
                     "SerialSpdDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");

  if (shouldEquilibrate() && !equilibratedA_)
    std::cout << "WARNING!  SerialSpdDenseSolver<T>::solve: System should be equilibrated!" << std::endl;

  if (inverted()) {

    TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument, 
                        "SerialSpdDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");

    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
    
    if (RHS_->values()!=LHS_->values()) {
       (*LHS_) = (*RHS_); // Copy B to X if needed
    }
    INFO_ = 0;
    this->POTRS(Matrix_->UPLO(), numRowCols_, RHS_->numCols(), AF_, LDAF_, LHS_->values(), LHS_->stride(), &INFO_);
    if (INFO_!=0) return(INFO_);
    solved_ = true;

  }
  int ierr1=0;
  if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
  if (ierr1!=0) 
    return(ierr1);
  else
    return(ierr);
  
  if (equilibrate_) ierr1 = unequilibrateLHS();
  return(ierr1);
}
//=============================================================================

template<typename OrdinalType, typename ScalarType>
int SerialSpdDenseSolver<OrdinalType,ScalarType>::applyRefinement()
{
  TEST_FOR_EXCEPTION(!solved(), std::logic_error,
         "SerialSpdDenseSolver<T>::applyRefinement: Must have an existing solution!");
  TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
         "SerialSpdDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");

  OrdinalType NRHS = RHS_->numCols();
  FERR_.resize( NRHS );
  BERR_.resize( NRHS );
  allocateWORK();
  allocateIWORK();
  
  INFO_ = 0;
  this->PORFS(Matrix_->UPLO(), numRowCols_, NRHS, A_, LDA_, AF_, LDAF_,
        RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(), 
              &FERR_[0], &BERR_[0], &WORK_[0], &IWORK_[0], &INFO_);
  
  solutionErrorsEstimated_ = true;
  reciprocalConditionEstimated_ = true;
  solutionRefined_ = true;
  
  return(INFO_);

}

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

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 = RHS_->stride(), NRHS = RHS_->numCols();
  ScalarType * X = RHS_->values();
  ScalarType * ptr;
  for (j=0; j<NRHS; 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()
{
  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_);
}

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

template<typename OrdinalType, typename ScalarType>
int SerialSpdDenseSolver<OrdinalType,ScalarType>::reciprocalConditionEstimate(MagnitudeType & Value)
{
  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;
  this->POCON( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &IWORK_[0], &INFO_);
  reciprocalConditionEstimated_ = true;
  Value = RCOND_;

  return(INFO_);
}
//=============================================================================

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 */