42 #ifndef _TEUCHOS_SERIALDENSESOLVER_HPP_ 43 #define _TEUCHOS_SERIALDENSESOLVER_HPP_ 56 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 57 #include "Eigen/Dense" 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;
147 typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
148 typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
149 typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
150 typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
151 typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
152 typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
153 typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
154 typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
155 typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
156 typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
157 typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
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_);}
371 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> > {
455 template<
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),
491 template<
typename OrdinalType,
typename ScalarType>
497 template<
typename OrdinalType,
typename ScalarType>
500 LHS_ = Teuchos::null;
501 RHS_ = Teuchos::null;
502 reciprocalConditionEstimated_ =
false;
503 solutionRefined_ =
false;
505 solutionErrorsEstimated_ =
false;
506 equilibratedB_ =
false;
510 template<
typename OrdinalType,
typename ScalarType>
511 void SerialDenseSolver<OrdinalType,ScalarType>::resetMatrix()
514 equilibratedA_ =
false;
536 template<
typename OrdinalType,
typename ScalarType>
544 Min_MN_ = TEUCHOS_MIN(M_,N_);
553 template<
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!");
576 template<
typename OrdinalType,
typename ScalarType>
579 estimateSolutionErrors_ = flag;
582 refineSolution_ = refineSolution_ || flag;
586 template<
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();
610 if (ierr!=0)
return(ierr);
612 if ((
int)IPIV_.size()!=Min_MN_) IPIV_.resize( Min_MN_ );
616 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 617 EigenMatrixMap aMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
618 EigenPermutationMatrix p;
619 EigenOrdinalArray ind;
620 EigenOrdinalArrayMap ipivMap( &IPIV_[0], Min_MN_ );
623 p = lu_.permutationP();
626 EigenMatrix luMat = lu_.matrixLU();
627 for (OrdinalType j=0; j<aMap.outerSize(); j++) {
628 for (OrdinalType i=0; i<aMap.innerSize(); i++) {
629 aMap(i,j) = luMat(i,j);
632 for (OrdinalType i=0; i<ipivMap.innerSize(); i++) {
636 this->GETRF(M_, N_, AF_, LDAF_, &IPIV_[0], &INFO_);
646 template<
typename OrdinalType,
typename ScalarType>
660 ierr = equilibrateRHS();
661 equilibratedB_ =
true;
663 if (ierr != 0)
return(ierr);
666 std::logic_error,
"SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
668 "SerialDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
670 "SerialDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
672 if (shouldEquilibrate() && !equilibratedA_)
673 std::cout <<
"WARNING! SerialDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
678 "SerialDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
682 RHS_->values(), RHS_->stride(), 0.0, LHS_->values(), LHS_->stride());
683 if (INFO_!=0)
return(INFO_);
688 if (!factored()) factor();
690 if (RHS_->values()!=LHS_->values()) {
695 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 696 EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
697 EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
699 lhsMap = lu_.solve(rhsMap);
701 lu_.matrixLU().template triangularView<Eigen::Upper>().transpose().solveInPlace(lhsMap);
702 lu_.matrixLU().template triangularView<Eigen::UnitLower>().transpose().solveInPlace(lhsMap);
703 EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
706 lu_.matrixLU().template triangularView<Eigen::Upper>().adjoint().solveInPlace(lhsMap);
707 lu_.matrixLU().template triangularView<Eigen::UnitLower>().adjoint().solveInPlace(lhsMap);
708 EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
712 this->GETRS(ETranspChar[TRANS_], N_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
715 if (INFO_!=0)
return(INFO_);
720 if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
724 if (equilibrate_) ierr1 = unequilibrateLHS();
729 template<
typename OrdinalType,
typename ScalarType>
733 "SerialDenseSolver<T>::applyRefinement: Must have an existing solution!");
735 "SerialDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
737 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 742 OrdinalType NRHS = RHS_->numCols();
743 FERR_.resize( NRHS );
744 BERR_.resize( NRHS );
748 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GERFS_WORK( N_ );
749 this->GERFS(ETranspChar[TRANS_], N_, NRHS, A_, LDA_, AF_, LDAF_, &IPIV_[0],
750 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
751 &FERR_[0], &BERR_[0], &WORK_[0], &GERFS_WORK[0], &INFO_);
753 solutionErrorsEstimated_ =
true;
754 reciprocalConditionEstimated_ =
true;
755 solutionRefined_ =
true;
763 template<
typename OrdinalType,
typename ScalarType>
766 if (R_.size()!=0)
return(0);
772 this->GEEQU (M_, N_, AF_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
777 shouldEquilibrate_ =
true;
784 template<
typename OrdinalType,
typename ScalarType>
790 if (equilibratedA_)
return(0);
791 if (R_.size()==0) ierr = computeEquilibrateScaling();
792 if (ierr!=0)
return(ierr);
795 for (j=0; j<N_; j++) {
797 ScalarType s1 = C_[j];
798 for (i=0; i<M_; i++) {
799 *ptr = *ptr*s1*R_[i];
807 for (j=0; j<N_; j++) {
809 ptr1 = AF_ + j*LDAF_;
810 ScalarType s1 = C_[j];
811 for (i=0; i<M_; i++) {
812 *ptr = *ptr*s1*R_[i];
814 *ptr1 = *ptr1*s1*R_[i];
820 equilibratedA_ =
true;
827 template<
typename OrdinalType,
typename ScalarType>
833 if (equilibratedB_)
return(0);
834 if (R_.size()==0) ierr = computeEquilibrateScaling();
835 if (ierr!=0)
return(ierr);
837 MagnitudeType * R_tmp = &R_[0];
838 if (transpose_) R_tmp = &C_[0];
840 OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
841 ScalarType * B = RHS_->values();
843 for (j=0; j<NRHS; j++) {
845 for (i=0; i<M_; i++) {
846 *ptr = *ptr*R_tmp[i];
851 equilibratedB_ =
true;
858 template<
typename OrdinalType,
typename ScalarType>
863 if (!equilibratedB_)
return(0);
865 MagnitudeType * C_tmp = &C_[0];
866 if (transpose_) C_tmp = &R_[0];
868 OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
869 ScalarType * X = LHS_->values();
871 for (j=0; j<NLHS; j++) {
873 for (i=0; i<N_; i++) {
874 *ptr = *ptr*C_tmp[i];
884 template<
typename OrdinalType,
typename ScalarType>
888 if (!factored()) factor();
890 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 891 EigenMatrixMap invMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
892 EigenMatrix invMat = lu_.inverse();
893 for (OrdinalType j=0; j<invMap.outerSize(); j++) {
894 for (OrdinalType i=0; i<invMap.innerSize(); i++) {
895 invMap(i,j) = invMat(i,j);
915 this->GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_[0], LWORK_, &INFO_);
927 template<
typename OrdinalType,
typename ScalarType>
930 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 934 if (reciprocalConditionEstimated()) {
942 if (!factored()) ierr = factor();
943 if (ierr!=0)
return(ierr);
949 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GECON_WORK( 2*N_ );
950 this->GECON(
'1', N_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &GECON_WORK[0], &INFO_);
952 reciprocalConditionEstimated_ =
true;
960 template<
typename OrdinalType,
typename ScalarType>
963 if (Matrix_!=Teuchos::null) os <<
"Solver Matrix" << std::endl << *Matrix_ << std::endl;
964 if (Factor_!=Teuchos::null) os <<
"Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
965 if (LHS_ !=Teuchos::null) os <<
"Solver LHS" << std::endl << *LHS_ << std::endl;
966 if (RHS_ !=Teuchos::null) os <<
"Solver RHS" << std::endl << *RHS_ << std::endl;
std::vector< MagnitudeType > C() const
Returns a pointer to the column scale vector used for equilibration.
int equilibrateMatrix()
Equilibrates the this matrix.
T magnitudeType
Mandatory typedef for result of magnitude.
Templated serial dense matrix class.
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the transpose of this matrix...
Templated interface class to BLAS routines.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed)...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
bool solved()
Returns true if the current set of vectors has been solved.
bool solutionRefined()
Returns true if the current set of vectors has been refined.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
Object for storing data and providing functionality that is common to all computational classes...
OrdinalType numRows() const
Returns row dimension of system.
This structure defines some basic traits for a scalar field type.
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
MagnitudeType COLCND() const
Ratio of smallest to largest column scale factors for 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. ...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
Functionality and data that is common to all computational classes.
Templated interface class to LAPACK routines.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
Teuchos implementation details.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system...
bool inverted()
Returns true if matrix inverse has been computed (inverse available via getFactoredMatrix()).
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor...
std::vector< OrdinalType > IPIV() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
int invert()
Inverts the this matrix.
bool transpose()
Returns true if transpose of this matrix has and will be used.
int applyRefinement()
Apply Iterative Refinement.
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i)) of the this matrix.
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
SerialDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors()...
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
The Templated LAPACK Wrapper Class.
MagnitudeType ROWCND() const
Ratio of smallest to largest row scale factors for the this matrix (returns -1 if not yet computed)...
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
virtual ~SerialDenseSolver()
SerialDenseSolver destructor.
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
int equilibrateRHS()
Equilibrates the current RHS.
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.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed)...
Defines basic traits for the scalar field type.
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
Smart reference counting pointer class for automatic garbage collection.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GETRF.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
Reference-counted pointer class and non-member templated function implementations.
static T one()
Returns representation of one for this scalar type.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
A class for solving dense linear problems.
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
This class creates and provides basic support for dense rectangular matrix of templated type...