41 #ifndef STOKHOS_MUELU_QR_INTERFACE_DEF_HPP 42 #define STOKHOS_MUELU_QR_INTERFACE_DEF_HPP 44 #include "MueLu_QR_Interface_decl.hpp" 45 #include "MueLu_Exceptions.hpp" 49 #if defined(HAVE_STOKHOS_MUELU) 51 template <
class Scalar,
class Storage,
class LocalOrdinal>
52 QR_Interface< Sacado::PCE::OrthogPoly<Scalar, Storage>,
LocalOrdinal>::QR_Interface(
const size_t NSDim) : workSize_(NSDim), NSDim_(NSDim), info_(0) {
53 tau_ = ArrayRCP<Scalar>(NSDim);
54 work_ = ArrayRCP<Scalar>(NSDim);
57 template <
class Scalar,
class Storage,
class LocalOrdinal>
59 int pce_sz = localQR[0].size();
60 if (localQR.size()*pce_sz > localQR_.size())
61 localQR_.resize(localQR.size()*pce_sz);
63 for (
int i=0; i<localQR.size(); ++i) {
64 for (
int j=0;
j<pce_sz; ++
j) {
65 localQR_[i*pce_sz+
j] = (localQR[i]).coeff(
j);
68 lapack_.GEQRF( myAggSize*pce_sz, Teuchos::as<int>(NSDim_), localQR_.getRawPtr(), myAggSize*pce_sz,
69 tau_.getRawPtr(), work_.getRawPtr(), workSize_, &info_ );
71 std::string msg =
"QR_Interface: dgeqrf (LAPACK QR routine) returned error code " + Teuchos::toString(info_);
72 throw(Exceptions::RuntimeError(msg));
75 for (
int i=0; i<localQR.size(); ++i) {
76 localQR[i].copyForWrite();
77 if (localQR[i].size() < pce_sz)
78 localQR[i].reset(localQR[0].expansion());
79 for (
int j=0;
j<pce_sz; ++
j)
80 localQR[i].fastAccessCoeff(
j) = localQR_[i*pce_sz+
j];
86 if (
std::abs(work_[0]) > workSize_) {
88 work_ = ArrayRCP<Scalar>(workSize_);
92 template <
class Scalar,
class Storage,
class LocalOrdinal>
97 int pce_sz = localQR[0].size();
98 LapackQR( lapack_, myAggSize*pce_sz, Teuchos::as<int>(NSDim_), localQR_, tau_, work_, workSize_, info_ );
100 std::string msg =
"QR_Interface: dorgqr (LAPACK auxiliary QR routine) returned error code " + Teuchos::toString(info_);
101 throw(Exceptions::RuntimeError(msg));
104 for (
int i=0; i<localQR.size(); ++i) {
105 localQR[i].copyForWrite();
106 if (localQR[i].size() < pce_sz)
107 localQR[i].reset(localQR[0].expansion());
108 for (
int j=0;
j<pce_sz; ++
j)
109 localQR[i].fastAccessCoeff(
j) = localQR_[i*pce_sz+
j];
115 if (
std::abs(work_[0]) > workSize_) {
117 work_ = ArrayRCP<Scalar>(workSize_);
120 #endif //if defined(HAVE_STOKHOS_MUELU) 125 #endif // STOKHOS_MUELU_QR_INTERFACE_DEF_HPP
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)