42#ifndef _TEUCHOS_SERIALQRDENSESOLVER_MP_VECTOR_HPP_
43#define _TEUCHOS_SERIALQRDENSESOLVER_MP_VECTOR_HPP_
49#include "Teuchos_SerialQRDenseSolver.hpp"
59 template<
typename OrdinalType,
typename Storage>
85 int setMatrix(
const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);
91 int setVectors(
const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
92 const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
102 base_QR_.factorWithEquilibration(flag);
107 base_QR_.solveWithTranspose(flag);
112 base_QR_.solveWithTranspose(trans);
124 int factor() {
return base_QR_.factor(); }
130 int solve() {
return base_QR_.solve(); }
137 return base_QR_.computeEquilibrateScaling();
177 int multiplyQ (ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C);
183 int solveR (ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C);
205 bool solved() {
return base_QR_.solved(); }
219 RCP<SerialDenseMatrix<OrdinalType, ScalarType> >
getMatrix()
const {
return(Matrix_);}
225 RCP<SerialDenseMatrix<OrdinalType, ScalarType> >
getQ()
const {
return(FactorQ_);}
228 RCP<SerialDenseMatrix<OrdinalType, ScalarType> >
getR()
const {
return(FactorR_);}
231 RCP<SerialDenseMatrix<OrdinalType, ScalarType> >
getLHS()
const {
return(LHS_);}
234 RCP<SerialDenseMatrix<OrdinalType, ScalarType> >
getRHS()
const {
return(RHS_);}
243 std::vector<ScalarType>
tau()
const {
return base_QR_.tau(); }
253 void Print(std::ostream& os)
const;
260 typedef SerialDenseMatrix<OrdinalType, ScalarType>
MatrixType;
264 RCP<BaseMatrixType> createBaseMatrix(
const RCP<MatrixType>& mat)
const;
265 RCP<MatrixType> createMatrix(
const RCP<BaseMatrixType>& base_mat)
const;
300 template <
typename Ordinal,
typename Value,
typename Device>
317 new (v_vec+i)
Vector(vector_size, v+i*vector_size,
false);
323 template <
typename Ordinal,
typename Value,
typename Device,
int Num>
330 return reinterpret_cast<Value*
>(v);
335 return reinterpret_cast<Vector*
>(v);
342 using namespace details;
346template<
typename OrdinalType,
typename Storage>
350 Object(
"Teuchos::SerialQRDenseSolver"),
359template<
typename OrdinalType,
typename Storage>
363 LHS_ = Teuchos::null;
364 RHS_ = Teuchos::null;
368template<
typename OrdinalType,
typename Storage>
378template<
typename OrdinalType,
typename Storage>
379RCP<SerialDenseMatrix<OrdinalType, typename Sacado::MP::Vector<Storage>::value_type> >
382 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& mat)
const
386 mat->values(), mat->stride()*mat->numCols());
387 RCP<BaseMatrixType> base_mat =
390 mat->stride()*SacadoSize_,
391 mat->numRows()*SacadoSize_,
397template<
typename OrdinalType,
typename Storage>
398RCP<SerialDenseMatrix<OrdinalType, Sacado::MP::Vector<Storage> > >
401 const RCP<SerialDenseMatrix<OrdinalType,BaseScalarType> >& base_mat)
const
405 base_mat->values(), base_mat->stride()*base_mat->numCols(), SacadoSize_);
406 RCP<MatrixType> mat =
409 base_mat->stride()/SacadoSize_,
410 base_mat->numRows()/SacadoSize_,
411 base_mat->numCols()));
416template<
typename OrdinalType,
typename Storage>
418setMatrix(
const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)
420 TEUCHOS_TEST_FOR_EXCEPTION(
421 A->numRows() < A->numCols(), std::invalid_argument,
422 "SerialQRDenseSolver<T>::setMatrix: the matrix A must have A.numRows() >= A.numCols()!");
431 if (Storage::is_static)
432 SacadoSize_ = Storage::static_size;
433 else if (M_ > 0 && N_ > 0)
434 SacadoSize_ = (*A)(0,0).size();
438 return base_QR_.setMatrix( createBaseMatrix(A) );
442template<
typename OrdinalType,
typename Storage>
445 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
446 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
449 TEUCHOS_TEST_FOR_EXCEPTION(
450 B->numCols() != X->numCols(), std::invalid_argument,
451 "SerialQRDenseSolver<T>::setVectors: X and B have different numbers of columns!");
452 TEUCHOS_TEST_FOR_EXCEPTION(
453 B->values()==0, std::invalid_argument,
454 "SerialQRDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
455 TEUCHOS_TEST_FOR_EXCEPTION(
456 X->values()==0, std::invalid_argument,
457 "SerialQRDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
458 TEUCHOS_TEST_FOR_EXCEPTION(
459 B->stride()<1, std::invalid_argument,
460 "SerialQRDenseSolver<T>::setVectors: B has an invalid stride!");
461 TEUCHOS_TEST_FOR_EXCEPTION(
462 X->stride()<1, std::invalid_argument,
463 "SerialQRDenseSolver<T>::setVectors: X has an invalid stride!");
469 return base_QR_.setVectors( createBaseMatrix(X), createBaseMatrix(B) );
474template<
typename OrdinalType,
typename Storage>
477 int ret = base_QR_.formQ();
478 FactorQ_ = createMatrix( base_QR_.getQ() );
484template<
typename OrdinalType,
typename Storage>
487 int ret = base_QR_.formR();
488 FactorR_ = createMatrix( base_QR_.getR() );
489 Factor_ = createMatrix( base_QR_.getFactoredMatrix() );
495template<
typename OrdinalType,
typename Storage>
497multiplyQ(ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C)
499 return base_QR_.multiplyQ(transq, createBaseMatrix(C));
504template<
typename OrdinalType,
typename Storage>
506solveR(ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C)
508 return base_QR_.solveR(transr, createBaseMatrix(C));
513template<
typename OrdinalType,
typename Storage>
515Print(std::ostream& os)
const {
517 if (Matrix_!=Teuchos::null) os <<
"Solver Matrix" << std::endl << *Matrix_ << std::endl;
518 if (Factor_!=Teuchos::null) os <<
"Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
519 if (LHS_ !=Teuchos::null) os <<
"Solver LHS" << std::endl << *LHS_ << std::endl;
520 if (RHS_ !=Teuchos::null) os <<
"Solver RHS" << std::endl << *RHS_ << std::endl;
Statically allocated storage class.
SerialQRDenseSolver & operator=(const SerialQRDenseSolver &Source)
SerialDenseMatrix< OrdinalType, BaseScalarType > BaseMatrixType
SerialDenseMatrix< OrdinalType, ScalarType > MatrixType
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the adjoint of this matrix,...
OrdinalType numRows() const
Returns row dimension of system.
RCP< MatrixType > Matrix_
RCP< BaseMatrixType > Base_RHS_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
ScalarType::value_type BaseScalarType
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
RCP< MatrixType > FactorR_
virtual ~SerialQRDenseSolver()
SerialQRDenseSolver destructor.
RCP< MatrixType > Factor_
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
int equilibrateRHS()
Equilibrates the current RHS.
ScalarTraits< ScalarType >::magnitudeType MagnitudeType
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
RCP< BaseMatrixType > Base_Factor_
bool formedR()
Returns true if R has been formed explicitly.
OrdinalType numCols() const
Returns column dimension of system.
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
Sacado::MP::Vector< Storage > ScalarType
bool solved()
Returns true if the current set of vectors has been solved.
RCP< BaseMatrixType > Base_LHS_
RCP< BaseMatrixType > Base_FactorQ_
MagnitudeType ANORM() const
Returns the absolute value of the largest element of this matrix (returns -1 if not yet computed).
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
RCP< BaseMatrixType > Base_FactorR_
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
int computeEquilibrateScaling()
Determines if this matrix should be scaled.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
bool transpose()
Returns true if adjoint of this matrix has and will be used.
RCP< BaseMatrixType > Base_Matrix_
RCP< MatrixType > FactorQ_
SerialQRDenseSolver(const SerialQRDenseSolver &Source)
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
int equilibrateMatrix()
Equilibrates the this matrix.
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getR() const
Returns pointer to R (assuming factorization has been performed).
SerialQRDenseSolver< OrdinalType, BaseScalarType > BaseQRType
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
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.
Specialization for Sacado::UQ::PCE< Storage<...> >
Top-level namespace for Stokhos classes and functions.
static Vector * getVectorArray(Value *v, Ordinal len, Ordinal vector_size)
Stokhos::DynamicStorage< Ordinal, Value, Device > Storage
Sacado::MP::Vector< Storage > Vector
static Value * getValueArray(Vector *v, Ordinal len)
static Value * getValueArray(Vector *v, Ordinal len)
Sacado::MP::Vector< Storage > Vector
static Vector * getVectorArray(Value *v, Ordinal len, Ordinal vector_size)
Stokhos::StaticFixedStorage< Ordinal, Value, Num, Device > Storage