42 #ifndef _TEUCHOS_SERIALBANDDENSESOLVER_HPP_ 43 #define _TEUCHOS_SERIALBANDDENSESOLVER_HPP_ 59 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 60 #include "Eigen/Dense" 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;
177 typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
178 typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
179 typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
180 typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
181 typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
182 typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
183 typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
184 typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
185 typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
186 typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
187 typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
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_);}
390 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_;
463 template<
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),
500 template<
typename OrdinalType,
typename ScalarType>
506 template<
typename OrdinalType,
typename ScalarType>
511 reciprocalConditionEstimated_ =
false;
512 solutionRefined_ =
false;
514 solutionErrorsEstimated_ =
false;
515 equilibratedB_ =
false;
519 template<
typename OrdinalType,
typename ScalarType>
523 equilibratedA_ =
false;
545 template<
typename OrdinalType,
typename ScalarType>
549 OrdinalType m = AB->numRows();
550 OrdinalType
n = AB->numCols();
551 OrdinalType kl = AB->lowerBandwidth();
552 OrdinalType ku = AB->upperBandwidth();
556 "SerialBandDenseSolver<T>::setMatrix: A is an empty SerialBandDenseMatrix<T>!");
575 template<
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!");
598 template<
typename OrdinalType,
typename ScalarType>
601 estimateSolutionErrors_ = flag;
604 refineSolution_ = refineSolution_ || flag;
608 template<
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();
628 if (ierr!=0)
return(ierr);
630 if (IPIV_.size()==0) IPIV_.resize( N_ );
634 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 635 EigenMatrixMap aMap( AF_+KL_, KL_+KU_+1, N_, EigenOuterStride(LDAF_) );
636 EigenMatrix fullMatrix(M_, N_);
637 for (OrdinalType j=0; j<N_; j++) {
639 fullMatrix(i,j) = aMap(KU_+i-j, j);
643 EigenPermutationMatrix p;
644 EigenOrdinalArray ind;
645 EigenOrdinalArrayMap ipivMap( &IPIV_[0], Min_MN_ );
647 lu_.compute(fullMatrix);
648 p = lu_.permutationP();
651 for (OrdinalType i=0; i<ipivMap.innerSize(); i++) {
656 this->GBTRF(M_, N_, KL_, KU_, AF_, LDAF_, &IPIV_[0], &INFO_);
666 template<
typename OrdinalType,
typename ScalarType>
676 ierr = equilibrateRHS();
677 equilibratedB_ =
true;
679 if (ierr != 0)
return(ierr);
682 std::logic_error,
"SerialBandDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
684 "SerialBandDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
686 "SerialBandDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
688 if (shouldEquilibrate() && !equilibratedA_)
689 std::cout <<
"WARNING! SerialBandDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
691 if (!factored()) factor();
693 if (RHS_->values()!=LHS_->values()) {
698 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 699 EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
700 EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
702 lhsMap = lu_.solve(rhsMap);
704 lu_.matrixLU().template triangularView<Eigen::Upper>().transpose().solveInPlace(lhsMap);
705 lu_.matrixLU().template triangularView<Eigen::UnitLower>().transpose().solveInPlace(lhsMap);
706 EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
709 lu_.matrixLU().template triangularView<Eigen::Upper>().adjoint().solveInPlace(lhsMap);
710 lu_.matrixLU().template triangularView<Eigen::UnitLower>().adjoint().solveInPlace(lhsMap);
711 EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
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();
732 template<
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 745 OrdinalType NRHS = RHS_->numCols();
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;
766 template<
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;
787 template<
typename OrdinalType,
typename ScalarType>
793 if (equilibratedA_)
return(0);
794 if (R_.size()==0) ierr = computeEquilibrateScaling();
795 if (ierr!=0)
return(ierr);
800 for (j=0; j<N_; j++) {
802 ScalarType s1 = C_[j];
804 *ptr = *ptr*s1*R_[i];
813 for (j=0; j<N_; j++) {
815 ScalarType s1 = C_[j];
817 *ptr = *ptr*s1*R_[i];
821 for (j=0; j<N_; j++) {
823 ptrL = AF_ + KL_ + KU_ + 1 + j*LDAF_;
824 ScalarType s1 = C_[j];
826 *ptrU = *ptrU*s1*R_[i];
830 *ptrL = *ptrL*s1*R_[i];
836 equilibratedA_ =
true;
843 template<
typename OrdinalType,
typename ScalarType>
849 if (equilibratedB_)
return(0);
850 if (R_.size()==0) ierr = computeEquilibrateScaling();
851 if (ierr!=0)
return(ierr);
854 if (transpose_) R_tmp = &C_[0];
856 OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
857 ScalarType *
B = RHS_->values();
859 for (j=0; j<NRHS; j++) {
861 for (i=0; i<M_; i++) {
862 *ptr = *ptr*R_tmp[i];
867 equilibratedB_ =
true;
875 template<
typename OrdinalType,
typename ScalarType>
880 if (!equilibratedB_)
return(0);
883 if (transpose_) C_tmp = &R_[0];
885 OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
886 ScalarType * X = LHS_->values();
888 for (j=0; j<NLHS; j++) {
890 for (i=0; i<N_; i++) {
891 *ptr = *ptr*C_tmp[i];
901 template<
typename OrdinalType,
typename ScalarType>
904 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 909 if (reciprocalConditionEstimated()) {
917 if (!factored()) ierr = factor();
918 if (ierr!=0)
return(ierr);
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;
934 template<
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;
bool equilibratedA()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
OrdinalType numUpper() const
Returns upper bandwidth of system.
std::vector< int > IWORK_
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
T magnitudeType
Mandatory typedef for result of magnitude.
Templated serial dense matrix class.
int applyRefinement()
Apply Iterative Refinement.
SerialBandDenseSolver & operator=(const SerialBandDenseSolver< OrdinalType, ScalarType > &Source)
void solveToRefinedSolution(bool flag)
Set whether or not to use iterative refinement to improve solutions to linear systems.
int setMatrix(const RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > &AB)
Sets the pointer for coefficient matrix.
Templated interface class to BLAS routines.
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
MagnitudeType ROWCND() const
Ratio of smallest to largest row scale factors for the this matrix (returns -1 if not yet computed)...
std::vector< MagnitudeType > C_
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
OrdinalType numCols() const
Returns column dimension of system.
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > RHS_
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
std::vector< OrdinalType > IPIV() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
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 unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
OrdinalType numRows() const
Returns row dimension of system.
int equilibrateMatrix()
Equilibrates the this matrix.
Object for storing data and providing functionality that is common to all computational classes...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
This structure defines some basic traits for a scalar field type.
bool reciprocalConditionEstimated_
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i)) of the this matrix.
bool solutionRefined()
Returns true if the current set of vectors has been refined.
std::vector< MagnitudeType > FERR_
Templated class for solving dense linear problems.
Functionality and data that is common to all computational classes.
Templated interface class to LAPACK routines.
Templated serial dense matrix class.
bool estimateSolutionErrors_
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
Teuchos implementation details.
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
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.
void solveWithTranspose(bool flag)
This class creates and provides basic support for banded dense matrices of templated type...
MagnitudeType COLCND() const
Ratio of smallest to largest column scale factors for the this matrix (returns -1 if not yet computed...
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
A class for representing and solving banded dense linear systems.
The Templated LAPACK Wrapper Class.
#define TEUCHOS_MAX(x, y)
bool transpose()
Returns true if transpose of this matrix has and will be used.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > Factor_
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > Matrix_
bool factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed)...
bool solved()
Returns true if the current set of vectors has been solved.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GBTRF.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system...
virtual ~SerialBandDenseSolver()
SerialBandDenseSolver destructor.
std::vector< OrdinalType > IPIV_
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
SerialBandDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors()...
OrdinalType numLower() const
Returns lower bandwidth of system.
Defines basic traits for the scalar field type.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETranspChar[]
void solveWithTransposeFlag(Teuchos::ETransp trans)
std::vector< MagnitudeType > R_
Smart reference counting pointer class for automatic garbage collection.
#define TEUCHOS_MIN(x, y)
int equilibrateRHS()
Equilibrates the current RHS.
std::vector< ScalarType > WORK_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > LHS_
Reference-counted pointer class and non-member templated function implementations.
static T one()
Returns representation of one for this scalar type.
bool solutionErrorsEstimated_
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed)...
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
void factorWithEquilibration(bool flag)
This class creates and provides basic support for dense rectangular matrix of templated type...
std::vector< MagnitudeType > BERR_