43 #ifndef TEUCHOS_SERIALSPDDENSESOLVER_H 44 #define TEUCHOS_SERIALSPDDENSESOLVER_H 59 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 60 #include "Eigen/Dense" 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;
142 typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
143 typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
144 typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
145 typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
146 typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
147 typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
148 typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
149 typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
150 typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
151 typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
152 typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
311 OrdinalType
numRows()
const {
return(numRowCols_);}
314 OrdinalType
numCols()
const {
return(numRowCols_);}
317 MagnitudeType
ANORM()
const {
return(ANORM_);}
320 MagnitudeType
RCOND()
const {
return(RCOND_);}
325 MagnitudeType
SCOND() {
return(SCOND_);};
328 MagnitudeType
AMAX()
const {
return(AMAX_);}
331 std::vector<MagnitudeType>
FERR()
const {
return(FERR_);}
334 std::vector<MagnitudeType>
BERR()
const {
return(BERR_);}
337 std::vector<MagnitudeType>
R()
const {
return(R_);}
343 void Print(std::ostream& os)
const;
349 void allocateWORK() { LWORK_ = 4*numRowCols_; WORK_.resize( LWORK_ );
return;}
350 void allocateIWORK() { IWORK_.resize( numRowCols_ );
return;}
355 bool shouldEquilibrate_;
360 bool estimateSolutionErrors_;
361 bool solutionErrorsEstimated_;
364 bool reciprocalConditionEstimated_;
365 bool refineSolution_;
366 bool solutionRefined_;
368 OrdinalType numRowCols_;
374 std::vector<int> IWORK_;
376 MagnitudeType ANORM_;
377 MagnitudeType RCOND_;
378 MagnitudeType SCOND_;
388 std::vector<MagnitudeType> FERR_;
389 std::vector<MagnitudeType> BERR_;
390 std::vector<ScalarType> WORK_;
391 std::vector<MagnitudeType> R_;
392 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 393 Eigen::LLT<EigenMatrix,Eigen::Upper> lltu_;
394 Eigen::LLT<EigenMatrix,Eigen::Lower> lltl_;
410 template<
typename OrdinalType,
typename ScalarType>
414 shouldEquilibrate_(false),
415 equilibratedA_(false),
416 equilibratedB_(false),
419 estimateSolutionErrors_(false),
420 solutionErrorsEstimated_(false),
423 reciprocalConditionEstimated_(false),
424 refineSolution_(false),
425 solutionRefined_(false),
442 template<
typename OrdinalType,
typename ScalarType>
448 template<
typename OrdinalType,
typename ScalarType>
451 LHS_ = Teuchos::null;
452 RHS_ = Teuchos::null;
453 reciprocalConditionEstimated_ =
false;
454 solutionRefined_ =
false;
456 solutionErrorsEstimated_ =
false;
457 equilibratedB_ =
false;
461 template<
typename OrdinalType,
typename ScalarType>
465 equilibratedA_ =
false;
483 template<
typename OrdinalType,
typename ScalarType>
489 numRowCols_ = A->numRows();
498 template<
typename OrdinalType,
typename ScalarType>
504 "SerialSpdDenseSolver<T>::setVectors: X and B are not the same size!");
506 "SerialSpdDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
508 "SerialSpdDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
510 "SerialSpdDenseSolver<T>::setVectors: B has an invalid stride!");
512 "SerialSpdDenseSolver<T>::setVectors: X has an invalid stride!");
521 template<
typename OrdinalType,
typename ScalarType>
524 estimateSolutionErrors_ = flag;
527 refineSolution_ = refineSolution_ || flag;
531 template<
typename OrdinalType,
typename ScalarType>
537 "SerialSpdDenseSolver<T>::factor: Cannot factor an inverted matrix!");
539 ANORM_ = Matrix_->normOne();
546 if (refineSolution_ ) {
548 AF_ = Factor_->values();
549 LDAF_ = Factor_->stride();
555 if (ierr!=0)
return(ierr);
559 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 560 EigenMatrixMap aMap( AF_, numRowCols_, numRowCols_, EigenOuterStride(LDAF_) );
562 if (Matrix_->upper())
568 this->
POTRF(Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
578 template<
typename OrdinalType,
typename ScalarType>
593 equilibratedB_ =
true;
595 if (ierr != 0)
return(ierr);
598 std::logic_error,
"SerialSpdDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
600 "SerialSpdDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
602 "SerialSpdDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
605 std::cout <<
"WARNING! SerialSpdDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
610 "SerialSpdDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
614 numRowCols_, 1.0, AF_, LDAF_, RHS_->values(), RHS_->stride(), 0.0,
615 LHS_->values(), LHS_->stride());
616 if (INFO_!=0)
return(INFO_);
623 if (RHS_->values()!=LHS_->values()) {
628 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 629 EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
630 EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
632 if (Matrix_->upper())
633 lhsMap = lltu_.solve(rhsMap);
635 lhsMap = lltl_.solve(rhsMap);
638 this->
POTRS(Matrix_->UPLO(), numRowCols_, RHS_->numCols(), AF_, LDAF_, LHS_->values(), LHS_->stride(), &INFO_);
641 if (INFO_!=0)
return(INFO_);
655 template<
typename OrdinalType,
typename ScalarType>
659 "SerialSpdDenseSolver<T>::applyRefinement: Must have an existing solution!");
661 "SerialSpdDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
664 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 668 OrdinalType NRHS = RHS_->numCols();
669 FERR_.resize( NRHS );
670 BERR_.resize( NRHS );
675 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> PORFS_WORK( numRowCols_ );
676 this->
PORFS(Matrix_->UPLO(), numRowCols_, NRHS, A_, LDA_, AF_, LDAF_,
677 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
678 &FERR_[0], &BERR_[0], &WORK_[0], &PORFS_WORK[0], &INFO_);
680 solutionErrorsEstimated_ =
true;
681 reciprocalConditionEstimated_ =
true;
682 solutionRefined_ =
true;
690 template<
typename OrdinalType,
typename ScalarType>
693 if (R_.size()!=0)
return(0);
695 R_.resize( numRowCols_ );
698 this->
POEQU (numRowCols_, AF_, LDAF_, &R_[0], &SCOND_, &AMAX_, &INFO_);
702 shouldEquilibrate_ =
true;
709 template<
typename OrdinalType,
typename ScalarType>
715 if (equilibratedA_)
return(0);
717 if (ierr!=0)
return(ierr);
718 if (Matrix_->upper()) {
721 for (j=0; j<numRowCols_; j++) {
723 ScalarType s1 = R_[j];
724 for (i=0; i<=j; i++) {
725 *ptr = *ptr*s1*R_[i];
733 for (j=0; j<numRowCols_; j++) {
735 ptr1 = AF_ + j*LDAF_;
736 ScalarType s1 = R_[j];
737 for (i=0; i<=j; i++) {
738 *ptr = *ptr*s1*R_[i];
740 *ptr1 = *ptr1*s1*R_[i];
749 for (j=0; j<numRowCols_; j++) {
750 ptr = A_ + j + j*LDA_;
751 ScalarType s1 = R_[j];
752 for (i=j; i<numRowCols_; i++) {
753 *ptr = *ptr*s1*R_[i];
761 for (j=0; j<numRowCols_; j++) {
762 ptr = A_ + j + j*LDA_;
763 ptr1 = AF_ + j + j*LDAF_;
764 ScalarType s1 = R_[j];
765 for (i=j; i<numRowCols_; i++) {
766 *ptr = *ptr*s1*R_[i];
768 *ptr1 = *ptr1*s1*R_[i];
775 equilibratedA_ =
true;
782 template<
typename OrdinalType,
typename ScalarType>
788 if (equilibratedB_)
return(0);
790 if (ierr!=0)
return(ierr);
792 OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
793 ScalarType * B = RHS_->values();
795 for (j=0; j<NRHS; j++) {
797 for (i=0; i<numRowCols_; i++) {
803 equilibratedB_ =
true;
810 template<
typename OrdinalType,
typename ScalarType>
815 if (!equilibratedB_)
return(0);
817 OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
818 ScalarType * X = LHS_->values();
820 for (j=0; j<NLHS; j++) {
822 for (i=0; i<numRowCols_; i++) {
833 template<
typename OrdinalType,
typename ScalarType>
836 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 843 this->
POTRI( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
846 if (Matrix_->upper())
860 template<
typename OrdinalType,
typename ScalarType>
863 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 876 if (ierr!=0)
return(ierr);
883 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> POCON_WORK( numRowCols_ );
884 this->
POCON( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &POCON_WORK[0], &INFO_);
885 reciprocalConditionEstimated_ =
true;
893 template<
typename OrdinalType,
typename ScalarType>
896 if (Matrix_!=Teuchos::null) os <<
"Solver Matrix" << std::endl << *Matrix_ << std::endl;
897 if (Factor_!=Teuchos::null) os <<
"Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
898 if (LHS_ !=Teuchos::null) os <<
"Solver LHS" << std::endl << *LHS_ << std::endl;
899 if (RHS_ !=Teuchos::null) os <<
"Solver RHS" << std::endl << *RHS_ << std::endl;
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed)...
void POCON(const char UPLO, const OrdinalType n, const ScalarType *A, const OrdinalType lda, const ScalarType anorm, ScalarType *rcond, ScalarType *WORK, OrdinalType *IWORK, OrdinalType *info) const
Estimates the reciprocal of the condition number (1-norm) of a real symmetric positive definite matri...
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement. ...
A class for constructing and using Hermitian positive definite dense matrices.
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
T magnitudeType
Mandatory typedef for result of magnitude.
MagnitudeType SCOND()
Ratio of smallest to largest equilibration scale factors for the this matrix (returns -1 if not yet c...
OrdinalType numCols() const
Returns column dimension of system.
Templated serial dense matrix class.
bool transpose()
Returns true if transpose of this matrix has and will be used.
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
Templated interface class to BLAS routines.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
int factor()
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF...
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
bool factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
int invert()
Inverts the this matrix.
Object for storing data and providing functionality that is common to all computational classes...
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed)...
This structure defines some basic traits for a scalar field type.
int equilibrateMatrix()
Equilibrates the this matrix.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
Templated class for solving dense linear problems.
Functionality and data that is common to all computational classes.
Templated interface class to LAPACK routines.
void POTRF(const char UPLO, const OrdinalType n, ScalarType *A, const OrdinalType lda, OrdinalType *info) const
Computes Cholesky factorization of a real symmetric positive definite matrix A.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system...
void POTRS(const char UPLO, const OrdinalType n, const OrdinalType nrhs, const ScalarType *A, const OrdinalType lda, ScalarType *B, const OrdinalType ldb, OrdinalType *info) const
Solves a system of linear equations A*X=B, where A is a symmetric positive definite matrix factored b...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
int setMatrix(const RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > &A_in)
Sets the pointers for coefficient matrix.
virtual ~SerialSpdDenseSolver()
SerialSpdDenseSolver destructor.
Teuchos implementation details.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to 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).
void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const alpha_type alpha, const A_type *A, const OrdinalType lda, const B_type *B, const OrdinalType ldb, const beta_type beta, ScalarType *C, const OrdinalType ldc) const
General matrix-matrix multiply.
void PORFS(const char UPLO, const OrdinalType n, const OrdinalType nrhs, const ScalarType *A, const OrdinalType lda, const ScalarType *AF, const OrdinalType ldaf, const ScalarType *B, const OrdinalType ldb, ScalarType *X, const OrdinalType ldx, ScalarType *FERR, ScalarType *BERR, ScalarType *WORK, OrdinalType *IWORK, OrdinalType *info) const
Improves the computed solution to a system of linear equations when the coefficient matrix is symmetr...
int applyRefinement()
Apply Iterative Refinement.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
The Templated LAPACK Wrapper Class.
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
bool solved()
Returns true if the current set of vectors has been solved.
Templated serial, dense, symmetric matrix class.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
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...
bool inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
OrdinalType numRows() const
Returns row dimension of system.
Smart reference counting pointer class for automatic garbage collection.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the this matrix.
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
bool solutionRefined()
Returns true if the current set of vectors has been refined.
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
Reference-counted pointer class and non-member templated function implementations.
void POTRI(const char UPLO, const OrdinalType n, ScalarType *A, const OrdinalType lda, OrdinalType *info) const
Computes the inverse of a real symmetric positive definite matrix A using the Cholesky factorization ...
void POEQU(const OrdinalType n, const ScalarType *A, const OrdinalType lda, MagnitudeType *S, MagnitudeType *scond, MagnitudeType *amax, OrdinalType *info) const
Computes row and column scalings intended to equilibrate a symmetric positive definite matrix A and r...
static T one()
Returns representation of one for this scalar type.
SerialSpdDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors()...
int equilibrateRHS()
Equilibrates the current RHS.
This class creates and provides basic support for dense rectangular matrix of templated type...