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(); }
252 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>
320 new (v_pce+i)
pce_type(
cijk, vector_size, v+i*vector_size,
false);
329 using namespace details;
333 template<
typename OrdinalType,
typename Storage>
334 SerialQRDenseSolver<OrdinalType,Sacado::UQ::PCE<Storage> >::
335 SerialQRDenseSolver()
337 Object(
"Teuchos::SerialQRDenseSolver"),
346 template<
typename OrdinalType,
typename Storage>
350 LHS_ = Teuchos::null;
351 RHS_ = Teuchos::null;
355 template<
typename OrdinalType,
typename Storage>
365 template<
typename OrdinalType,
typename Storage>
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_,
384 template<
typename OrdinalType,
typename Storage>
385 RCP<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()));
403 template<
typename OrdinalType,
typename Storage>
405 setMatrix(
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) );
429 template<
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) );
461 template<
typename OrdinalType,
typename Storage>
464 int ret = base_QR_.formQ();
465 FactorQ_ = createMatrix( base_QR_.getQ() );
471 template<
typename OrdinalType,
typename Storage>
474 int ret = base_QR_.formR();
475 FactorR_ = createMatrix( base_QR_.getR() );
476 Factor_ = createMatrix( base_QR_.getFactoredMatrix() );
482 template<
typename OrdinalType,
typename Storage>
484 multiplyQ(ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C)
486 return base_QR_.multiplyQ(transq, createBaseMatrix(C));
491 template<
typename OrdinalType,
typename Storage>
493 solveR(ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C)
495 return base_QR_.solveR(transr, createBaseMatrix(C));
500 template<
typename OrdinalType,
typename Storage>
502 Print(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;
static Value * getValueArray(pce_type *v, Ordinal len)
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getR() const
Returns pointer to R (assuming factorization has been performed).
RCP< BaseMatrixType > Base_Matrix_
Stokhos::DynamicStorage< Ordinal, Value, Device > Storage
pce_type::cijk_type cijk_type
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
Sacado::ETPCE::OrthogPoly< double, Stokhos::StandardStorage< int, double > > pce_type
ScalarType::value_type BaseScalarType
static pce_type * getPCEArray(Value *v, Ordinal len, Ordinal vector_size)
OrdinalType numRows() const
Returns row dimension of system.
RCP< MatrixType > FactorQ_
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
RCP< MatrixType > Matrix_
virtual ~SerialQRDenseSolver()
SerialQRDenseSolver destructor.
RCP< MatrixType > Factor_
bool transpose()
Returns true if adjoint of this matrix has and will be used.
int equilibrateRHS()
Equilibrates the current RHS.
ScalarTraits< ScalarType >::magnitudeType MagnitudeType
SerialDenseMatrix< OrdinalType, BaseScalarType > BaseMatrixType
RCP< MatrixType > FactorR_
bool solved()
Returns true if the current set of vectors has been solved.
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
bool formedR()
Returns true if R has been formed explicitly.
RCP< BaseMatrixType > Base_LHS_
Specialization for Sacado::UQ::PCE< Storage<...> >
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
RCP< BaseMatrixType > Base_Factor_
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
MagnitudeType ANORM() const
Returns the absolute value of the largest element of this matrix (returns -1 if not yet computed)...
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor...
Sacado::UQ::PCE< Storage > pce_type
bool formedQ()
Returns true if Q has been formed explicitly.
Top-level namespace for Stokhos classes and functions.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
int computeEquilibrateScaling()
Determines if this matrix should be scaled.
RCP< BaseMatrixType > Base_FactorQ_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Sacado::UQ::PCE< Storage > ScalarType
std::vector< ScalarType > tau() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
int equilibrateMatrix()
Equilibrates the this matrix.
SerialQRDenseSolver< OrdinalType, BaseScalarType > BaseQRType
OrdinalType numCols() const
Returns column dimension of system.
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the adjoint of this matrix...
SerialDenseMatrix< OrdinalType, ScalarType > MatrixType
RCP< BaseMatrixType > Base_FactorR_
RCP< BaseMatrixType > Base_RHS_
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system...
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...