8 #ifndef MUELU_QR_INTERFACEEX_DEF_HPP_ 9 #define MUELU_QR_INTERFACEEX_DEF_HPP_ 17 template <
class Scalar,
class LocalOrdinal>
18 void LapackQR(Teuchos::LAPACK<LocalOrdinal,Scalar> &lapack, LocalOrdinal myAggSize,
19 int intFineNSDim, ArrayRCP<Scalar> &localQR, ArrayRCP<Scalar> &tau,
20 ArrayRCP<Scalar> &work, LocalOrdinal &workSize, LocalOrdinal &info)
22 lapack.ORGQR(myAggSize, intFineNSDim, intFineNSDim, localQR.getRawPtr(),
23 myAggSize, tau.getRawPtr(), work.getRawPtr(), workSize, &info );
27 template <
class LocalOrdinal>
28 void LapackQR(Teuchos::LAPACK<LocalOrdinal, std::complex<double> > &lapack,
29 LocalOrdinal myAggSize,
int intFineNSDim, ArrayRCP<std::complex<double> > &localQR,
30 ArrayRCP<std::complex<double> > &tau, ArrayRCP<std::complex<double> > &work,
31 LocalOrdinal &workSize, LocalOrdinal &info)
33 lapack.UNGQR(myAggSize, intFineNSDim, intFineNSDim, localQR.getRawPtr(),
34 myAggSize, tau.getRawPtr(), work.getRawPtr(), workSize, &info );
37 template <
class Scalar,
class LocalOrdinal>
39 tau_ = ArrayRCP<Scalar>(NSDim);
40 work_ = ArrayRCP<Scalar>(NSDim);
43 template <
class Scalar,
class LocalOrdinal>
45 bool bIsZeroColumn =
true;
46 for(LocalOrdinal r = 0; r < myAggSize; ++r) {
47 if(localQR[ myAggSize*nspCol + r ]!=0.0) {
48 bIsZeroColumn =
false;
55 template <
class Scalar,
class LocalOrdinal>
58 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
63 for (LocalOrdinal k=0; k<myAggSize; ++k) {dtemp += Teuchos::ScalarTraits<Scalar>::magnitude(localQR[k])*Teuchos::ScalarTraits<Scalar>::magnitude(localQR[k]);}
64 dtemp = Teuchos::ScalarTraits<Magnitude>::squareroot(dtemp);
71 for(LocalOrdinal j = 0; j<NSDim_; ++j) {
72 if(!isZeroNspColumn(myAggSize, localQR, j))
73 nonZeroCols_.push_back(j);
77 internalNSDim_ = Teuchos::as<LocalOrdinal>(nonZeroCols_.size());
80 tau_ = ArrayRCP<Scalar>(internalNSDim_);
81 work_ = ArrayRCP<Scalar>(internalNSDim_);
82 workSize_ = internalNSDim_;
86 internalLocalQR_ = Teuchos::ArrayRCP<Scalar>(myAggSize*internalNSDim_);
87 for(
size_t j = 0; j<internalNSDim_; ++j) {
88 for(
size_t r = 0; r<myAggSize; ++r) {
89 internalLocalQR_[j*myAggSize+r] = localQR[nonZeroCols_[j]*myAggSize+r];
103 lapack_.GEQRF( myAggSize, Teuchos::as<int>(internalNSDim_), internalLocalQR_.getRawPtr(), myAggSize,
104 tau_.getRawPtr(), work_.getRawPtr(), workSize_, &info_ );
106 std::string msg =
"QR_InterfaceEx: dgeqrf (LAPACK QR routine) returned error code " +
Teuchos::toString(info_);
114 if ( Teuchos::ScalarTraits<Scalar>::magnitude(work_[0]) > workSize_) {
115 workSize_ = Teuchos::as<int>(Teuchos::ScalarTraits<Scalar>::magnitude(work_[0]));
116 work_ = ArrayRCP<Scalar>(workSize_);
121 for(
size_t j = 0; j<internalNSDim_; ++j) {
122 for(
size_t r = 0; r<internalNSDim_; ++r) {
123 localQR[nonZeroCols_[j]*myAggSize+nonZeroCols_[r]] = internalLocalQR_[j*myAggSize+r];
138 template <
class Scalar,
class LocalOrdinal>
141 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
144 Magnitude dtemp = Teuchos::ScalarTraits<Scalar>::magnitude(localQR[0]);
145 localQR[0] = tau_[0];
147 for (LocalOrdinal i=0; i<myAggSize; ++i)
154 LapackQR( lapack_, myAggSize, Teuchos::as<int>(internalNSDim_), internalLocalQR_, tau_, work_, workSize_, info_ );
156 std::string msg =
"QR_InterfaceEx: dorgqr (LAPACK auxiliary QR routine) returned error code " +
Teuchos::toString(info_);
162 if ( Teuchos::ScalarTraits<Scalar>::magnitude(work_[0]) > workSize_) {
163 workSize_ = Teuchos::as<int>(Teuchos::ScalarTraits<Scalar>::magnitude(work_[0]));
164 work_ = ArrayRCP<Scalar>(workSize_);
169 for(
size_t j = 0; j<internalNSDim_; ++j) {
170 for(
size_t r = 0; r<myAggSize; ++r) {
171 localQR[nonZeroCols_[j]*myAggSize+r] = internalLocalQR_[j*myAggSize+r];
QR_InterfaceEx(const size_t nullSpaceDimension)
Constructor.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Namespace for MueLu classes and methods.
ArrayRCP< Scalar > tau_
Internal LAPACK variables.
void ExtractQ(LocalOrdinal const &myAggSize, ArrayRCP< Scalar > &localQR)
Calculate the Q factor.
ArrayRCP< Scalar > work_
Temporary work space.
void Compute(LocalOrdinal const &myAggSize, ArrayRCP< Scalar > &localQR)
Compute the QR factorization.
void LapackQR(Teuchos::LAPACK< LocalOrdinal, Scalar > &lapack, LocalOrdinal myAggSize, int intFineNSDim, ArrayRCP< Scalar > &localQR, ArrayRCP< Scalar > &tau, ArrayRCP< Scalar > &work, LocalOrdinal &workSize, LocalOrdinal &info)
Non-member templated function to handle extracting Q from QR factorization.
bool isZeroNspColumn(LocalOrdinal const &myAggSize, ArrayRCP< Scalar > &localQR, LocalOrdinal nspCol)
Exception throws to report errors in the internal logical of the program.