42#ifndef _TEUCHOS_SERIALQRDENSESOLVER_UQ_PCE_HPP_
43#define _TEUCHOS_SERIALQRDENSESOLVER_UQ_PCE_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>
318 cijk_type cijk = Kokkos::getGlobalCijkTensor<cijk_type>();
320 new (v_pce+i)
pce_type(cijk, vector_size, v+i*vector_size,
false);
329 using namespace details;
333template<
typename OrdinalType,
typename Storage>
337 Object(
"Teuchos::SerialQRDenseSolver"),
346template<
typename OrdinalType,
typename Storage>
350 LHS_ = Teuchos::null;
351 RHS_ = Teuchos::null;
355template<
typename OrdinalType,
typename Storage>
365template<
typename OrdinalType,
typename Storage>
366RCP<SerialDenseMatrix<OrdinalType, typename Sacado::UQ::PCE<Storage>::value_type> >
369 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& mat)
const
373 mat->values(), mat->stride()*mat->numCols());
374 RCP<BaseMatrixType> base_mat =
377 mat->stride()*SacadoSize_,
378 mat->numRows()*SacadoSize_,
384template<
typename OrdinalType,
typename Storage>
385RCP<SerialDenseMatrix<OrdinalType, Sacado::UQ::PCE<Storage> > >
388 const RCP<SerialDenseMatrix<OrdinalType,BaseScalarType> >& base_mat)
const
392 base_mat->values(), base_mat->stride()*base_mat->numCols(), SacadoSize_);
393 RCP<MatrixType> mat =
396 base_mat->stride()/SacadoSize_,
397 base_mat->numRows()/SacadoSize_,
398 base_mat->numCols()));
403template<
typename OrdinalType,
typename Storage>
405setMatrix(
const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)
407 TEUCHOS_TEST_FOR_EXCEPTION(
408 A->numRows() < A->numCols(), std::invalid_argument,
409 "SerialQRDenseSolver<T>::setMatrix: the matrix A must have A.numRows() >= A.numCols()!");
418 if (Storage::is_static)
419 SacadoSize_ = Storage::static_size;
420 else if (M_ > 0 && N_ > 0)
421 SacadoSize_ = (*A)(0,0).size();
425 return base_QR_.setMatrix( createBaseMatrix(A) );
429template<
typename OrdinalType,
typename Storage>
432 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
433 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
436 TEUCHOS_TEST_FOR_EXCEPTION(
437 B->numCols() != X->numCols(), std::invalid_argument,
438 "SerialQRDenseSolver<T>::setVectors: X and B have different numbers of columns!");
439 TEUCHOS_TEST_FOR_EXCEPTION(
440 B->values()==0, std::invalid_argument,
441 "SerialQRDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
442 TEUCHOS_TEST_FOR_EXCEPTION(
443 X->values()==0, std::invalid_argument,
444 "SerialQRDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
445 TEUCHOS_TEST_FOR_EXCEPTION(
446 B->stride()<1, std::invalid_argument,
447 "SerialQRDenseSolver<T>::setVectors: B has an invalid stride!");
448 TEUCHOS_TEST_FOR_EXCEPTION(
449 X->stride()<1, std::invalid_argument,
450 "SerialQRDenseSolver<T>::setVectors: X has an invalid stride!");
456 return base_QR_.setVectors( createBaseMatrix(X), createBaseMatrix(B) );
461template<
typename OrdinalType,
typename Storage>
464 int ret = base_QR_.formQ();
465 FactorQ_ = createMatrix( base_QR_.getQ() );
471template<
typename OrdinalType,
typename Storage>
474 int ret = base_QR_.formR();
475 FactorR_ = createMatrix( base_QR_.getR() );
476 Factor_ = createMatrix( base_QR_.getFactoredMatrix() );
482template<
typename OrdinalType,
typename Storage>
484multiplyQ(ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C)
486 return base_QR_.multiplyQ(transq, createBaseMatrix(C));
491template<
typename OrdinalType,
typename Storage>
493solveR(ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C)
495 return base_QR_.solveR(transr, createBaseMatrix(C));
500template<
typename OrdinalType,
typename Storage>
502Print(std::ostream& os)
const {
504 if (Matrix_!=Teuchos::null) os <<
"Solver Matrix" << std::endl << *Matrix_ << std::endl;
505 if (Factor_!=Teuchos::null) os <<
"Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
506 if (LHS_ !=Teuchos::null) os <<
"Solver LHS" << std::endl << *LHS_ << std::endl;
507 if (RHS_ !=Teuchos::null) os <<
"Solver RHS" << std::endl << *RHS_ << std::endl;
SerialQRDenseSolver< OrdinalType, BaseScalarType > BaseQRType
SerialQRDenseSolver & operator=(const SerialQRDenseSolver &Source)
RCP< BaseMatrixType > Base_RHS_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
int equilibrateRHS()
Equilibrates the current RHS.
Sacado::UQ::PCE< Storage > ScalarType
bool formedR()
Returns true if R has been formed explicitly.
bool formedQ()
Returns true if Q has been formed explicitly.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
ScalarTraits< ScalarType >::magnitudeType MagnitudeType
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
int computeEquilibrateScaling()
Determines if this matrix should be scaled.
RCP< BaseMatrixType > Base_Factor_
RCP< BaseMatrixType > Base_Matrix_
OrdinalType numCols() const
Returns column dimension of system.
bool solved()
Returns true if the current set of vectors has been solved.
int equilibrateMatrix()
Equilibrates the this matrix.
virtual ~SerialQRDenseSolver()
SerialQRDenseSolver destructor.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
RCP< BaseMatrixType > Base_LHS_
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
SerialDenseMatrix< OrdinalType, BaseScalarType > BaseMatrixType
SerialDenseMatrix< OrdinalType, ScalarType > MatrixType
RCP< MatrixType > Matrix_
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the adjoint of this matrix,...
std::vector< ScalarType > tau() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
ScalarType::value_type BaseScalarType
RCP< MatrixType > FactorR_
bool transpose()
Returns true if adjoint of this matrix has and will be used.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
SerialQRDenseSolver(const SerialQRDenseSolver &Source)
MagnitudeType ANORM() const
Returns the absolute value of the largest element of this matrix (returns -1 if not yet computed).
RCP< BaseMatrixType > Base_FactorR_
OrdinalType numRows() const
Returns row dimension of system.
RCP< MatrixType > FactorQ_
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
RCP< BaseMatrixType > Base_FactorQ_
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
RCP< MatrixType > Factor_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getR() const
Returns pointer to R (assuming factorization has been performed).
Specialization for Sacado::UQ::PCE< Storage<...> >
Sacado::ETPCE::OrthogPoly< double, Stokhos::StandardStorage< int, double > > pce_type
Top-level namespace for Stokhos classes and functions.
static pce_type * getPCEArray(Value *v, Ordinal len, Ordinal vector_size)
Stokhos::DynamicStorage< Ordinal, Value, Device > Storage
Sacado::UQ::PCE< Storage > pce_type
static Value * getValueArray(pce_type *v, Ordinal len)
pce_type::cijk_type cijk_type