42#ifndef _TEUCHOS_SERIALBANDDENSESOLVER_HPP_
43#define _TEUCHOS_SERIALBANDDENSESOLVER_HPP_
59#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
165 template<
typename OrdinalType,
typename ScalarType>
167 public LAPACK<OrdinalType, ScalarType>
173#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
174 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor>
EigenMatrix;
175 typedef typename Eigen::internal::BandMatrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor>
EigenBandMatrix;
176 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1>
EigenVector;
179 typedef typename Eigen::Map<EigenVector,0,EigenInnerStride >
EigenVectorMap;
181 typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride >
EigenMatrixMap;
183 typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType>
EigenPermutationMatrix;
186 typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1>
EigenScalarArray;
376 std::vector<MagnitudeType>
FERR()
const {
return(
FERR_);}
379 std::vector<MagnitudeType>
BERR()
const {
return(
BERR_);}
382 std::vector<MagnitudeType>
R()
const {
return(
R_);}
385 std::vector<MagnitudeType>
C()
const {
return(
C_);}
391 void Print(std::ostream&
os)
const;
444 std::vector<MagnitudeType>
R_;
445 std::vector<MagnitudeType>
C_;
446#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
447 Eigen::PartialPivLU<EigenMatrix>
lu_;
463template<
typename OrdinalType,
typename ScalarType>
467 shouldEquilibrate_(
false),
468 equilibratedA_(
false),
469 equilibratedB_(
false),
472 estimateSolutionErrors_(
false),
473 solutionErrorsEstimated_(
false),
475 reciprocalConditionEstimated_(
false),
476 refineSolution_(
false),
477 solutionRefined_(
false),
500template<
typename OrdinalType,
typename ScalarType>
506template<
typename OrdinalType,
typename ScalarType>
511 reciprocalConditionEstimated_ =
false;
512 solutionRefined_ =
false;
514 solutionErrorsEstimated_ =
false;
515 equilibratedB_ =
false;
519template<
typename OrdinalType,
typename ScalarType>
523 equilibratedA_ =
false;
545template<
typename OrdinalType,
typename ScalarType>
556 "SerialBandDenseSolver<T>::setMatrix: A is an empty SerialBandDenseMatrix<T>!");
575template<
typename OrdinalType,
typename ScalarType>
581 "SerialBandDenseSolver<T>::setVectors: X and B are not the same size!");
583 "SerialBandDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
585 "SerialBandDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
587 "SerialBandDenseSolver<T>::setVectors: B has an invalid stride!");
589 "SerialBandDenseSolver<T>::setVectors: X has an invalid stride!");
598template<
typename OrdinalType,
typename ScalarType>
601 estimateSolutionErrors_ =
flag;
604 refineSolution_ = refineSolution_ ||
flag;
608template<
typename OrdinalType,
typename ScalarType>
611 if (factored())
return(0);
613 ANORM_ = Matrix_->normOne();
619 if (refineSolution_ ) {
621 AF_ = Factor_->values();
622 LDAF_ = Factor_->stride();
626 if (equilibrate_)
ierr = equilibrateMatrix();
630 if (IPIV_.size()==0) IPIV_.resize( N_ );
634#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
648 p =
lu_.permutationP();
656 this->GBTRF(M_, N_, KL_, KU_, AF_, LDAF_, &IPIV_[0], &INFO_);
666template<
typename OrdinalType,
typename ScalarType>
676 ierr = equilibrateRHS();
681 "SerialBandDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
683 "SerialBandDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
685 if (!factored()) factor();
688 std::logic_error,
"SerialBandDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
690 if (shouldEquilibrate() && !equilibratedA_)
691 std::cout <<
"WARNING! SerialBandDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
693 if (RHS_->values()!=LHS_->values()) {
698#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
715 this->GBTRS(
ETranspChar[TRANS_], N_, KL_, KU_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
718 if (INFO_!=0)
return(INFO_);
722 if (refineSolution_)
ierr1 = applyRefinement();
726 if (equilibrate_)
ierr1 = unequilibrateLHS();
732template<
typename OrdinalType,
typename ScalarType>
736 "SerialBandDenseSolver<T>::applyRefinement: Must have an existing solution!");
738 "SerialBandDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
740#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
746 FERR_.resize(
NRHS );
747 BERR_.resize(
NRHS );
751 std::vector<typename details::lapack_traits<ScalarType>::iwork_type>
GBRFS_WORK( N_ );
752 this->GBRFS(
ETranspChar[TRANS_], N_, KL_, KU_,
NRHS, A_+KL_, LDA_, AF_, LDAF_, &IPIV_[0],
753 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
754 &FERR_[0], &BERR_[0], &WORK_[0], &
GBRFS_WORK[0], &INFO_);
756 solutionErrorsEstimated_ =
true;
757 reciprocalConditionEstimated_ =
true;
758 solutionRefined_ =
true;
766template<
typename OrdinalType,
typename ScalarType>
769 if (R_.size()!=0)
return(0);
775 this->GBEQU (M_, N_, KL_, KU_, AF_+KL_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
780 shouldEquilibrate_ =
true;
787template<
typename OrdinalType,
typename ScalarType>
793 if (equilibratedA_)
return(0);
794 if (R_.size()==0)
ierr = computeEquilibrateScaling();
800 for (
j=0;
j<N_;
j++) {
804 *ptr = *ptr*
s1*R_[
i];
813 for (
j=0;
j<N_;
j++) {
817 *ptr = *ptr*
s1*R_[
i];
821 for (
j=0;
j<N_;
j++) {
823 ptrL = AF_ + KL_ + KU_ + 1 +
j*LDAF_;
836 equilibratedA_ =
true;
843template<
typename OrdinalType,
typename ScalarType>
849 if (equilibratedB_)
return(0);
850 if (R_.size()==0)
ierr = computeEquilibrateScaling();
854 if (transpose_)
R_tmp = &C_[0];
861 for (
i=0;
i<M_;
i++) {
867 equilibratedB_ =
true;
875template<
typename OrdinalType,
typename ScalarType>
880 if (!equilibratedB_)
return(0);
883 if (transpose_)
C_tmp = &R_[0];
890 for (
i=0;
i<N_;
i++) {
901template<
typename OrdinalType,
typename ScalarType>
904#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
909 if (reciprocalConditionEstimated()) {
917 if (!factored())
ierr = factor();
924 std::vector<typename details::lapack_traits<ScalarType>::iwork_type>
GBCON_WORK( N_ );
925 this->GBCON(
'1', N_, KL_, KU_, AF_, LDAF_, &IPIV_[0], ANORM_, &RCOND_, &WORK_[0], &
GBCON_WORK[0], &INFO_);
926 reciprocalConditionEstimated_ =
true;
934template<
typename OrdinalType,
typename ScalarType>
937 if (Matrix_!=
Teuchos::null)
os <<
"Solver Matrix" << std::endl << *Matrix_ << std::endl;
938 if (Factor_!=
Teuchos::null)
os <<
"Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
939 if (LHS_ !=
Teuchos::null)
os <<
"Solver LHS" << std::endl << *LHS_ << std::endl;
940 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.
#define TEUCHOS_MIN(x, y)
#define TEUCHOS_MAX(x, y)
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.
Templated serial dense matrix class.
Templated class for solving dense linear problems.
Functionality and data that is common to all computational classes.
The Templated LAPACK Wrapper Class.
Smart reference counting pointer class for automatic garbage collection.
A class for representing and solving banded dense linear systems.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
void solveWithTranspose(bool flag)
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
SerialBandDenseSolver(const SerialBandDenseSolver< OrdinalType, ScalarType > &Source)
MagnitudeType ROWCND() const
Ratio of smallest to largest row scale factors for the this matrix (returns -1 if not yet computed).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
std::vector< MagnitudeType > R_
OrdinalType numRows() const
Returns row dimension of system.
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).
std::vector< MagnitudeType > FERR_
bool factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
bool solutionErrorsEstimated_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > RHS_
std::vector< OrdinalType > IPIV() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
int applyRefinement()
Apply Iterative Refinement.
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
std::vector< int > IWORK_
OrdinalType numLower() const
Returns lower bandwidth of system.
bool solved()
Returns true if the current set of vectors has been solved.
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed).
OrdinalType numUpper() const
Returns upper bandwidth of system.
virtual ~SerialBandDenseSolver()
SerialBandDenseSolver destructor.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > LHS_
int setMatrix(const RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > &AB)
Sets the pointer for coefficient matrix.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
void factorWithEquilibration(bool flag)
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GBTRF.
bool estimateSolutionErrors_
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
SerialBandDenseSolver & operator=(const SerialBandDenseSolver< OrdinalType, ScalarType > &Source)
std::vector< OrdinalType > IPIV_
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > Matrix_
bool transpose()
Returns true if transpose of this matrix has and will be used.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed).
std::vector< MagnitudeType > C() const
Returns a pointer to the column scale vector used for equilibration.
MagnitudeType COLCND() const
Ratio of smallest to largest column scale factors for the this matrix (returns -1 if not yet computed...
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > Factor_
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i)) of the this matrix.
std::vector< ScalarType > WORK_
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
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.
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
std::vector< MagnitudeType > C_
int equilibrateMatrix()
Equilibrates the this matrix.
bool solutionRefined()
Returns true if the current set of vectors has been refined.
std::vector< MagnitudeType > BERR_
OrdinalType numCols() const
Returns column dimension of system.
bool reciprocalConditionEstimated_
int equilibrateRHS()
Equilibrates the current RHS.
void solveWithTransposeFlag(Teuchos::ETransp trans)
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
void solveToRefinedSolution(bool flag)
Set whether or not to use iterative refinement to improve solutions to linear systems.
SerialBandDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
Concrete serial communicator subclass.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETranspChar[]
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
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.