42#ifndef _TEUCHOS_SERIALDENSESOLVER_HPP_
43#define _TEUCHOS_SERIALDENSESOLVER_HPP_
56#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
136 template<
typename OrdinalType,
typename ScalarType>
138 public LAPACK<OrdinalType, ScalarType>
144#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
145 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor>
EigenMatrix;
146 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1>
EigenVector;
149 typedef typename Eigen::Map<EigenVector,0,EigenInnerStride >
EigenVectorMap;
151 typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride >
EigenMatrixMap;
153 typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType>
EigenPermutationMatrix;
156 typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1>
EigenScalarArray;
335 std::vector<OrdinalType>
IPIV()
const {
return(IPIV_);}
338 MagnitudeType
ANORM()
const {
return(ANORM_);}
341 MagnitudeType
RCOND()
const {
return(RCOND_);}
346 MagnitudeType
ROWCND()
const {
return(ROWCND_);}
351 MagnitudeType
COLCND()
const {
return(COLCND_);}
354 MagnitudeType
AMAX()
const {
return(AMAX_);}
357 std::vector<MagnitudeType>
FERR()
const {
return(FERR_);}
360 std::vector<MagnitudeType>
BERR()
const {
return(BERR_);}
363 std::vector<MagnitudeType>
R()
const {
return(R_);}
366 std::vector<MagnitudeType>
C()
const {
return(C_);}
372 void Print(std::ostream& os)
const;
376 void allocateWORK() { LWORK_ = 4*N_; WORK_.resize( LWORK_ );
return;}
382 bool shouldEquilibrate_;
387 bool estimateSolutionErrors_;
388 bool solutionErrorsEstimated_;
391 bool reciprocalConditionEstimated_;
392 bool refineSolution_;
393 bool solutionRefined_;
405 std::vector<OrdinalType> IPIV_;
407 MagnitudeType ANORM_;
408 MagnitudeType RCOND_;
409 MagnitudeType ROWCND_;
410 MagnitudeType COLCND_;
413 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Matrix_;
414 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
415 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
416 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Factor_;
420 std::vector<MagnitudeType> FERR_;
421 std::vector<MagnitudeType> BERR_;
422 std::vector<ScalarType> WORK_;
423 std::vector<MagnitudeType> R_;
424 std::vector<MagnitudeType> C_;
425#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
426 Eigen::PartialPivLU<EigenMatrix> lu_;
433 SerialDenseSolver & operator=(
const SerialDenseSolver<OrdinalType, ScalarType>& Source);
441 struct lapack_traits {
442 typedef int iwork_type;
447 struct lapack_traits<std::complex<T> > {
455template<
typename OrdinalType,
typename ScalarType>
459 shouldEquilibrate_(
false),
460 equilibratedA_(
false),
461 equilibratedB_(
false),
464 estimateSolutionErrors_(
false),
465 solutionErrorsEstimated_(
false),
468 reciprocalConditionEstimated_(
false),
469 refineSolution_(
false),
470 solutionRefined_(
false),
491template<
typename OrdinalType,
typename ScalarType>
497template<
typename OrdinalType,
typename ScalarType>
500 LHS_ = Teuchos::null;
501 RHS_ = Teuchos::null;
502 reciprocalConditionEstimated_ =
false;
503 solutionRefined_ =
false;
505 solutionErrorsEstimated_ =
false;
506 equilibratedB_ =
false;
510template<
typename OrdinalType,
typename ScalarType>
511void SerialDenseSolver<OrdinalType,ScalarType>::resetMatrix()
514 equilibratedA_ =
false;
536template<
typename OrdinalType,
typename ScalarType>
544 Min_MN_ = TEUCHOS_MIN(M_,N_);
553template<
typename OrdinalType,
typename ScalarType>
559 "SerialDenseSolver<T>::setVectors: X and B are not the same size!");
561 "SerialDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
563 "SerialDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
565 "SerialDenseSolver<T>::setVectors: B has an invalid stride!");
567 "SerialDenseSolver<T>::setVectors: X has an invalid stride!");
576template<
typename OrdinalType,
typename ScalarType>
579 estimateSolutionErrors_ =
flag;
582 refineSolution_ = refineSolution_ ||
flag;
586template<
typename OrdinalType,
typename ScalarType>
589 if (factored())
return(0);
592 "SerialDenseSolver<T>::factor: Cannot factor an inverted matrix!");
594 ANORM_ = Matrix_->normOne();
601 if (refineSolution_ ) {
603 AF_ = Factor_->values();
604 LDAF_ = Factor_->stride();
608 if (equilibrate_)
ierr = equilibrateMatrix();
612 if ((
int)IPIV_.size()!=Min_MN_) IPIV_.resize( Min_MN_ );
616#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
623 p =
lu_.permutationP();
636 this->GETRF(M_, N_, AF_, LDAF_, &IPIV_[0], &INFO_);
646template<
typename OrdinalType,
typename ScalarType>
660 ierr = equilibrateRHS();
665 "SerialDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
667 "SerialDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
672 "SerialDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
674 std::logic_error,
"SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
678 RHS_->values(), RHS_->stride(), 0.0, LHS_->values(), LHS_->stride());
679 if (INFO_!=0)
return(INFO_);
684 if (!factored()) factor();
687 std::logic_error,
"SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
689 if (RHS_->values()!=LHS_->values()) {
694#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
711 this->GETRS(ETranspChar[TRANS_], N_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
714 if (INFO_!=0)
return(INFO_);
720 if (shouldEquilibrate() && !equilibratedA_)
721 std::cout <<
"WARNING! SerialDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
724 if (refineSolution_ && !inverted())
ierr1 = applyRefinement();
728 if (equilibrate_)
ierr1 = unequilibrateLHS();
733template<
typename OrdinalType,
typename ScalarType>
737 "SerialDenseSolver<T>::applyRefinement: Must have an existing solution!");
739 "SerialDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
741#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
747 FERR_.resize(
NRHS );
748 BERR_.resize(
NRHS );
752 std::vector<typename details::lapack_traits<ScalarType>::iwork_type>
GERFS_WORK( N_ );
753 this->GERFS(ETranspChar[TRANS_], N_,
NRHS, A_, LDA_, AF_, LDAF_, &IPIV_[0],
754 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
755 &FERR_[0], &BERR_[0], &WORK_[0], &
GERFS_WORK[0], &INFO_);
757 solutionErrorsEstimated_ =
true;
758 reciprocalConditionEstimated_ =
true;
759 solutionRefined_ =
true;
767template<
typename OrdinalType,
typename ScalarType>
770 if (R_.size()!=0)
return(0);
776 this->GEEQU (M_, N_, AF_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
781 shouldEquilibrate_ =
true;
788template<
typename OrdinalType,
typename ScalarType>
794 if (equilibratedA_)
return(0);
795 if (R_.size()==0)
ierr = computeEquilibrateScaling();
799 for (
j=0;
j<N_;
j++) {
802 for (
i=0;
i<M_;
i++) {
811 for (
j=0;
j<N_;
j++) {
813 ptr1 = AF_ +
j*LDAF_;
815 for (
i=0;
i<M_;
i++) {
824 equilibratedA_ =
true;
831template<
typename OrdinalType,
typename ScalarType>
837 if (equilibratedB_)
return(0);
838 if (R_.size()==0)
ierr = computeEquilibrateScaling();
841 MagnitudeType *
R_tmp = &R_[0];
842 if (transpose_)
R_tmp = &C_[0];
849 for (
i=0;
i<M_;
i++) {
855 equilibratedB_ =
true;
862template<
typename OrdinalType,
typename ScalarType>
867 if (!equilibratedB_)
return(0);
869 MagnitudeType *
C_tmp = &C_[0];
870 if (transpose_)
C_tmp = &R_[0];
877 for (
i=0;
i<N_;
i++) {
888template<
typename OrdinalType,
typename ScalarType>
892 if (!factored()) factor();
894#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
919 this->GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_[0], LWORK_, &INFO_);
931template<
typename OrdinalType,
typename ScalarType>
934#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
938 if (reciprocalConditionEstimated()) {
946 if (!factored())
ierr = factor();
953 std::vector<typename details::lapack_traits<ScalarType>::iwork_type>
GECON_WORK( 2*N_ );
954 this->GECON(
'1', N_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &
GECON_WORK[0], &INFO_);
956 reciprocalConditionEstimated_ =
true;
964template<
typename OrdinalType,
typename ScalarType>
967 if (Matrix_!=Teuchos::null) os <<
"Solver Matrix" << std::endl << *Matrix_ << std::endl;
968 if (Factor_!=Teuchos::null) os <<
"Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
969 if (LHS_ !=Teuchos::null) os <<
"Solver LHS" << std::endl << *LHS_ << std::endl;
970 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.
Defines basic traits for the scalar field type.
Templated serial dense 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 solving dense linear problems.
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
MagnitudeType ROWCND() const
Ratio of smallest to largest row scale factors for the this matrix (returns -1 if not yet computed).
bool solutionRefined()
Returns true if the current set of vectors has been refined.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
OrdinalType numRows() const
Returns row dimension of system.
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
std::vector< MagnitudeType > C() const
Returns a pointer to the column scale vector used for equilibration.
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i)) of the this matrix.
int equilibrateRHS()
Equilibrates the current RHS.
bool inverted()
Returns true if matrix inverse has been computed (inverse available via getFactoredMatrix()).
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
MagnitudeType COLCND() const
Ratio of smallest to largest column scale factors for the this matrix (returns -1 if not yet computed...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
bool solved()
Returns true if the current set of vectors has been solved.
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
int invert()
Inverts the this matrix.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GETRF.
int equilibrateMatrix()
Equilibrates the this matrix.
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
int applyRefinement()
Apply Iterative Refinement.
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed).
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement.
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the transpose of this matrix,...
virtual ~SerialDenseSolver()
SerialDenseSolver destructor.
SerialDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with 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).
OrdinalType numCols() const
Returns column dimension of system.
bool transpose()
Returns true if transpose of this matrix has and will be used.
std::vector< OrdinalType > IPIV() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
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.
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
#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.