43#ifndef TEUCHOS_SERIALSPDDENSESOLVER_H
44#define TEUCHOS_SERIALSPDDENSESOLVER_H
59#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
132 template<
typename OrdinalType,
typename ScalarType>
134 public LAPACK<OrdinalType, ScalarType>
139#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
140 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor>
EigenMatrix;
141 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1>
EigenVector;
144 typedef typename Eigen::Map<EigenVector,0,EigenInnerStride >
EigenVectorMap;
146 typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride >
EigenMatrixMap;
148 typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType>
EigenPermutationMatrix;
151 typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1>
EigenScalarArray;
316 MagnitudeType
ANORM()
const {
return(ANORM_);}
319 MagnitudeType
RCOND()
const {
return(RCOND_);}
324 MagnitudeType
SCOND() {
return(SCOND_);};
327 MagnitudeType
AMAX()
const {
return(AMAX_);}
330 std::vector<MagnitudeType>
FERR()
const {
return(FERR_);}
333 std::vector<MagnitudeType>
BERR()
const {
return(BERR_);}
336 std::vector<MagnitudeType>
R()
const {
return(R_);}
343 void Print(std::ostream& os)
const;
348 void allocateWORK() { LWORK_ = 4*numRowCols_; WORK_.resize( LWORK_ );
return;}
349 void allocateIWORK() { IWORK_.resize( numRowCols_ );
return;}
354 bool shouldEquilibrate_;
359 bool estimateSolutionErrors_;
360 bool solutionErrorsEstimated_;
363 bool reciprocalConditionEstimated_;
364 bool refineSolution_;
365 bool solutionRefined_;
367 OrdinalType numRowCols_;
373 std::vector<int> IWORK_;
375 MagnitudeType ANORM_;
376 MagnitudeType RCOND_;
377 MagnitudeType SCOND_;
380 RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > Matrix_;
381 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
382 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
383 RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > Factor_;
387 std::vector<MagnitudeType> FERR_;
388 std::vector<MagnitudeType> BERR_;
389 std::vector<ScalarType> WORK_;
390 std::vector<MagnitudeType> R_;
391#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
392 Eigen::LLT<EigenMatrix,Eigen::Upper> lltu_;
393 Eigen::LLT<EigenMatrix,Eigen::Lower> lltl_;
409template<
typename OrdinalType,
typename ScalarType>
413 shouldEquilibrate_(
false),
414 equilibratedA_(
false),
415 equilibratedB_(
false),
418 estimateSolutionErrors_(
false),
419 solutionErrorsEstimated_(
false),
422 reciprocalConditionEstimated_(
false),
423 refineSolution_(
false),
424 solutionRefined_(
false),
441template<
typename OrdinalType,
typename ScalarType>
447template<
typename OrdinalType,
typename ScalarType>
450 LHS_ = Teuchos::null;
451 RHS_ = Teuchos::null;
452 reciprocalConditionEstimated_ =
false;
453 solutionRefined_ =
false;
455 solutionErrorsEstimated_ =
false;
456 equilibratedB_ =
false;
460template<
typename OrdinalType,
typename ScalarType>
461void SerialSpdDenseSolver<OrdinalType,ScalarType>::resetMatrix()
464 equilibratedA_ =
false;
482template<
typename OrdinalType,
typename ScalarType>
488 numRowCols_ =
A->numRows();
497template<
typename OrdinalType,
typename ScalarType>
503 "SerialSpdDenseSolver<T>::setVectors: X and B are not the same size!");
505 "SerialSpdDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
507 "SerialSpdDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
509 "SerialSpdDenseSolver<T>::setVectors: B has an invalid stride!");
511 "SerialSpdDenseSolver<T>::setVectors: X has an invalid stride!");
520template<
typename OrdinalType,
typename ScalarType>
523 estimateSolutionErrors_ =
flag;
526 refineSolution_ = refineSolution_ ||
flag;
530template<
typename OrdinalType,
typename ScalarType>
533 if (factored())
return(0);
536 "SerialSpdDenseSolver<T>::factor: Cannot factor an inverted matrix!");
538 ANORM_ = Matrix_->normOne();
545 if (refineSolution_ ) {
547 AF_ = Factor_->values();
548 LDAF_ = Factor_->stride();
552 if (equilibrate_)
ierr = equilibrateMatrix();
558#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
561 if (Matrix_->upper())
567 this->POTRF(Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
577template<
typename OrdinalType,
typename ScalarType>
591 ierr = equilibrateRHS();
596 "SerialSpdDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
598 "SerialSpdDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
603 "SerialSpdDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
605 std::logic_error,
"SerialSpdDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
609 numRowCols_, 1.0, AF_, LDAF_, RHS_->values(), RHS_->stride(), 0.0,
610 LHS_->values(), LHS_->stride());
611 if (INFO_!=0)
return(INFO_);
616 if (!factored()) factor();
619 std::logic_error,
"SerialSpdDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
621 if (RHS_->values()!=LHS_->values()) {
626#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
630 if (Matrix_->upper())
636 this->POTRS(Matrix_->UPLO(), numRowCols_, RHS_->numCols(), AF_, LDAF_, LHS_->values(), LHS_->stride(), &INFO_);
639 if (INFO_!=0)
return(INFO_);
645 if (shouldEquilibrate() && !equilibratedA_)
646 std::cout <<
"WARNING! SerialSpdDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
649 if (refineSolution_ && !inverted())
ierr1 = applyRefinement();
653 if (equilibrate_)
ierr1 = unequilibrateLHS();
658template<
typename OrdinalType,
typename ScalarType>
662 "SerialSpdDenseSolver<T>::applyRefinement: Must have an existing solution!");
664 "SerialSpdDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
667#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
672 FERR_.resize(
NRHS );
673 BERR_.resize(
NRHS );
678 std::vector<typename details::lapack_traits<ScalarType>::iwork_type>
PORFS_WORK( numRowCols_ );
679 this->PORFS(Matrix_->UPLO(), numRowCols_,
NRHS, A_, LDA_, AF_, LDAF_,
680 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
681 &FERR_[0], &BERR_[0], &WORK_[0], &
PORFS_WORK[0], &INFO_);
683 solutionErrorsEstimated_ =
true;
684 reciprocalConditionEstimated_ =
true;
685 solutionRefined_ =
true;
693template<
typename OrdinalType,
typename ScalarType>
696 if (R_.size()!=0)
return(0);
698 R_.resize( numRowCols_ );
701 this->POEQU (numRowCols_, AF_, LDAF_, &R_[0], &SCOND_, &AMAX_, &INFO_);
705 shouldEquilibrate_ =
true;
712template<
typename OrdinalType,
typename ScalarType>
718 if (equilibratedA_)
return(0);
719 if (R_.size()==0)
ierr = computeEquilibrateScaling();
721 if (Matrix_->upper()) {
724 for (
j=0;
j<numRowCols_;
j++) {
727 for (
i=0;
i<=
j;
i++) {
736 for (
j=0;
j<numRowCols_;
j++) {
738 ptr1 = AF_ +
j*LDAF_;
740 for (
i=0;
i<=
j;
i++) {
752 for (
j=0;
j<numRowCols_;
j++) {
753 ptr = A_ +
j +
j*LDA_;
755 for (
i=
j;
i<numRowCols_;
i++) {
764 for (
j=0;
j<numRowCols_;
j++) {
765 ptr = A_ +
j +
j*LDA_;
768 for (
i=
j;
i<numRowCols_;
i++) {
778 equilibratedA_ =
true;
785template<
typename OrdinalType,
typename ScalarType>
791 if (equilibratedB_)
return(0);
792 if (R_.size()==0)
ierr = computeEquilibrateScaling();
800 for (
i=0;
i<numRowCols_;
i++) {
806 equilibratedB_ =
true;
813template<
typename OrdinalType,
typename ScalarType>
818 if (!equilibratedB_)
return(0);
825 for (
i=0;
i<numRowCols_;
i++) {
836template<
typename OrdinalType,
typename ScalarType>
839#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
843 if (!factored()) factor();
846 this->POTRI( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
849 if (Matrix_->upper())
863template<
typename OrdinalType,
typename ScalarType>
866#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
870 if (reciprocalConditionEstimated()) {
878 if (!factored())
ierr = factor();
886 std::vector<typename details::lapack_traits<ScalarType>::iwork_type>
POCON_WORK( numRowCols_ );
887 this->POCON( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &
POCON_WORK[0], &INFO_);
888 reciprocalConditionEstimated_ =
true;
896template<
typename OrdinalType,
typename ScalarType>
899 if (Matrix_!=Teuchos::null) os <<
"Solver Matrix" << std::endl << *Matrix_ << std::endl;
900 if (Factor_!=Teuchos::null) os <<
"Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
901 if (LHS_ !=Teuchos::null) os <<
"Solver LHS" << std::endl << *LHS_ << std::endl;
902 if (RHS_ !=Teuchos::null) os <<
"Solver RHS" << std::endl << *RHS_ << std::endl;
Templated interface class to BLAS routines.
Object for storing data and providing functionality that is common to all computational classes.
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
Templated interface class to LAPACK routines.
Reference-counted pointer class and non-member templated function implementations.
Templated serial dense matrix class.
Templated class for solving dense linear problems.
Templated serial, dense, symmetric matrix class.
Functionality and data that is common to all computational classes.
The Templated LAPACK Wrapper Class.
Smart reference counting pointer class for automatic garbage collection.
RCP(ENull null_arg=null)
Initialize RCP<T> to NULL.
RCP< T > rcp(const boost::shared_ptr< T > &sptr)
Conversion function that takes in a boost::shared_ptr object and spits out a Teuchos::RCP object.
Ptr< T > ptr() const
Get a safer wrapper raw C++ pointer to the underlying object.
A class for constructing and using Hermitian positive definite dense matrices.
MagnitudeType SCOND()
Ratio of smallest to largest equilibration scale factors for the this matrix (returns -1 if not yet c...
int setMatrix(const RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > &A_in)
Sets the pointers for coefficient matrix.
OrdinalType numCols() const
Returns column dimension of system.
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
int applyRefinement()
Apply Iterative Refinement.
bool transpose()
Returns true if transpose of this matrix has and will be used.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
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).
int factor()
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
bool solutionRefined()
Returns true if the current set of vectors has been refined.
virtual ~SerialSpdDenseSolver()
SerialSpdDenseSolver destructor.
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
bool solved()
Returns true if the current set of vectors has been solved.
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement.
OrdinalType numRows() const
Returns row dimension of system.
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed).
bool factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
SerialSpdDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
int equilibrateMatrix()
Equilibrates the this matrix.
int invert()
Inverts the this matrix.
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
bool inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the this matrix.
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed).
int equilibrateRHS()
Equilibrates the current RHS.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos,...
Teuchos implementation details.
This structure defines some basic traits for a scalar field type.
static T one()
Returns representation of one for this scalar type.
T magnitudeType
Mandatory typedef for result of magnitude.