42#ifndef _TEUCHOS_SERIALQRDENSESOLVER_HPP_
43#define _TEUCHOS_SERIALQRDENSESOLVER_HPP_
57#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
131 template<
typename OrdinalType,
typename ScalarType>
133 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;
144 typedef typename Eigen::Map<EigenVector,0,EigenInnerStride >
EigenVectorMap;
146 typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride >
EigenMatrixMap;
148 typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType>
EigenPermutationMatrix;
151 typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1>
EigenScalarArray;
320 std::vector<ScalarType>
tau()
const {
return(
TAU_);}
330 void Print(std::ostream&
os)
const;
382#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
383 Eigen::HouseholderQR<EigenMatrix>
qr_;
399template<
typename OrdinalType,
typename ScalarType>
403 shouldEquilibrate_(
false),
404 equilibratedA_(
false),
405 equilibratedB_(
false),
406 equilibrationModeA_(0),
407 equilibrationModeB_(0),
411 matrixIsZero_(
false),
435template<
typename OrdinalType,
typename ScalarType>
441template<
typename OrdinalType,
typename ScalarType>
448 equilibratedB_ =
false;
452template<
typename OrdinalType,
typename ScalarType>
456 equilibratedA_ =
false;
457 equilibrationModeA_ = 0;
458 equilibrationModeB_ = 0;
460 matrixIsZero_ =
false;
481template<
typename OrdinalType,
typename ScalarType>
485 "SerialQRDenseSolver<T>::setMatrix: the matrix A must have A.numRows() >= A.numCols()!");
508template<
typename OrdinalType,
typename ScalarType>
514 "SerialQRDenseSolver<T>::setVectors: X and B have different numbers of columns!");
516 "SerialQRDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
518 "SerialQRDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
520 "SerialQRDenseSolver<T>::setVectors: B has an invalid stride!");
522 "SerialQRDenseSolver<T>::setVectors: X has an invalid stride!");
532template<
typename OrdinalType,
typename ScalarType>
535 if (factored())
return(0);
539 if (equilibrate_)
ierr = equilibrateMatrix();
543 if ((
int)TAU_.size()!=Min_MN_) TAU_.resize( Min_MN_ );
548#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
564 this->GEQRF(M_, N_, AF_, LDAF_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);
574template<
typename OrdinalType,
typename ScalarType>
586 ierr = equilibrateRHS();
591 "SerialQRDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
593 "SerialQRDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
596 "SerialQRDenseSolver<T>::solve: No transpose: must have LHS_->numRows() = N_");
599 "SerialQRDenseSolver<T>::solve: Transpose: must have LHS_->numRows() = M_");
602 if (shouldEquilibrate() && !equilibratedA_)
603 std::cout <<
"WARNING! SerialQRDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
606 if (!factored()) factor();
609 std::logic_error,
"SerialQRDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
614 (*TMP_)(
i,
j) = (*RHS_)(
i,
j);
641 (*LHS_)(
i,
j) = (*TMP_)(
i,
j);
649 ierr = unequilibrateLHS();
660template<
typename OrdinalType,
typename ScalarType>
675 for (
j = 0;
j < N_;
j++) {
676 for (
i = 0;
i < M_;
i++) {
685 shouldEquilibrate_ =
true;
686 }
else if (ANORM_ >
bignum) {
688 shouldEquilibrate_ =
true;
689 }
else if (ANORM_ ==
mZero) {
691 matrixIsZero_ =
true;
699template<
typename OrdinalType,
typename ScalarType>
702 if (equilibratedA_)
return(0);
717 for (
j = 0;
j < N_;
j++) {
718 for (
i = 0;
i < M_;
i++) {
728 equilibrationModeA_ = 1;
729 }
else if (ANORM_ >
bignum) {
732 equilibrationModeA_ = 2;
733 }
else if (ANORM_ ==
mZero) {
735 matrixIsZero_ =
true;
738 equilibratedA_ =
true;
745template<
typename OrdinalType,
typename ScalarType>
748 if (equilibratedB_)
return(0);
763 for (
j = 0;
j <RHS_->numCols();
j++) {
764 for (
i = 0;
i < RHS_->numRows();
i++) {
775 RHS_->values(), RHS_->stride(), &INFO_);
776 equilibrationModeB_ = 1;
777 }
else if (BNORM_ >
bignum) {
780 RHS_->values(), RHS_->stride(), &INFO_);
781 equilibrationModeB_ = 2;
784 equilibratedB_ =
true;
791template<
typename OrdinalType,
typename ScalarType>
794 if (!equilibratedB_)
return(0);
808 if (equilibrationModeA_ == 1) {
810 LHS_->values(), LHS_->stride(), &INFO_);
811 }
else if (equilibrationModeA_ == 2) {
813 LHS_->values(), LHS_->stride(), &INFO_);
815 if (equilibrationModeB_ == 1) {
817 LHS_->values(), LHS_->stride(), &INFO_);
818 }
else if (equilibrationModeB_ == 2) {
820 LHS_->values(), LHS_->stride(), &INFO_);
828template<
typename OrdinalType,
typename ScalarType>
832 if (!factored()) factor();
837 Q_ = FactorQ_->values();
838 LDQ_ = FactorQ_->stride();
844#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
853 this->UNGQR(M_, N_, N_, Q_, LDQ_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);
856 if (INFO_!=0)
return(INFO_);
865template<
typename OrdinalType,
typename ScalarType>
869 if (!factored()) factor();
874 R_ = FactorR_->values();
875 LDR_ = FactorR_->stride();
881 (*FactorR_)(
i,
j) = (*Factor_)(
i,
j);
892template<
typename OrdinalType,
typename ScalarType>
898 "SerialQRDenseSolver<T>::multiplyQ: C.numRows() != M_!");
900 "SerialQRDenseSolver<T>::multiplyQ: C is an empty SerialDenseMatrix<T>!");
903 if (!factored()) factor();
908#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
931 &TAU_[0],
C.values(),
C.stride(), &WORK_[0], LWORK_, &INFO_);
937 &TAU_[0],
C.values(),
C.stride(), &WORK_[0], LWORK_, &INFO_);
953template<
typename OrdinalType,
typename ScalarType>
959 "SerialQRDenseSolver<T>::solveR: must have N_ < C.numRows() < M_!");
961 "SerialQRDenseSolver<T>::solveR: C is an empty SerialDenseMatrix<T>!");
964 if (!factored()) factor();
969#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
1015template<
typename OrdinalType,
typename ScalarType>
1018 if (Matrix_!=
Teuchos::null)
os <<
"Solver Matrix" << std::endl << *Matrix_ << std::endl;
1019 if (Factor_!=
Teuchos::null)
os <<
"Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
1020 if (Q_!=
Teuchos::null)
os <<
"Solver Factor Q" << std::endl << *Q_ << std::endl;
1021 if (LHS_ !=
Teuchos::null)
os <<
"Solver LHS" << std::endl << *LHS_ << std::endl;
1022 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 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.
Concrete serial communicator subclass.
A class for solving dense linear problems.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > FactorR_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > LHS_
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
SerialQRDenseSolver & operator=(const SerialQRDenseSolver< OrdinalType, ScalarType > &Source)
int solveR(ETransp transr, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Solve input matrix on the left with the upper triangular matrix R or its adjoint.
bool solved()
Returns true if the current set of vectors has been solved.
std::vector< ScalarType > WORK_
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
bool formedR()
Returns true if R has been formed explicitly.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > FactorQ_
int formR()
Explicitly forms the upper triangular matrix R.
ScalarTraits< ScalarType >::magnitudeType MagnitudeType
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the 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).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > Matrix_
SerialQRDenseSolver(const SerialQRDenseSolver< OrdinalType, ScalarType > &Source)
OrdinalType equilibrationModeA_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the adjoint of this matrix,...
bool transpose()
Returns true if adjoint of this matrix has and will be used.
int multiplyQ(ETransp transq, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Left multiply the input matrix by the unitary matrix Q or its adjoint.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > TMP_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > RHS_
MagnitudeType ANORM() const
Returns the absolute value of the largest element of this matrix (returns -1 if not yet computed).
OrdinalType numCols() const
Returns column dimension of system.
bool formedQ()
Returns true if Q has been formed explicitly.
std::vector< ScalarType > tau() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
OrdinalType equilibrationModeB_
OrdinalType numRows() const
Returns row dimension of system.
virtual ~SerialQRDenseSolver()
SerialQRDenseSolver destructor.
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > Factor_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
SerialQRDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
std::vector< ScalarType > TAU_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getR() const
Returns pointer to R (assuming factorization has been performed).
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
int formQ()
Explicitly forms the unitary matrix Q.
int equilibrateMatrix()
Equilibrates the this matrix.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
int equilibrateRHS()
Equilibrates the current RHS.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
int computeEquilibrateScaling()
Determines if this matrix should be scaled.
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
#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.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char EUploChar[]
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETypeChar[]
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char EDiagChar[]
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ESideChar[]
Teuchos implementation details.
This structure defines some basic traits for a scalar field type.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
static T one()
Returns representation of one for this scalar type.
static T zero()
Returns representation of zero for this scalar type.
static magnitudeType sfmin()
Returns safe minimum (sfmin), such that 1/sfmin does not overflow.
static magnitudeType prec()
Returns eps*base.