MueLu  Version of the Day
MueLu_QR_InterfaceEx_def.hpp
Go to the documentation of this file.
1 /*
2  * MueLu_QR_InterfaceEx_def.hpp
3  *
4  * Created on: Jan 31, 2013
5  * Author: tobias
6  */
7 
8 #ifndef MUELU_QR_INTERFACEEX_DEF_HPP_
9 #define MUELU_QR_INTERFACEEX_DEF_HPP_
10 
12 #include "MueLu_Exceptions.hpp"
13 
14 namespace MueLu {
15 
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)
21  {
22  lapack.ORGQR(myAggSize, intFineNSDim, intFineNSDim, localQR.getRawPtr(),
23  myAggSize, tau.getRawPtr(), work.getRawPtr(), workSize, &info );
24  }
25 
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)
32  {
33  lapack.UNGQR(myAggSize, intFineNSDim, intFineNSDim, localQR.getRawPtr(),
34  myAggSize, tau.getRawPtr(), work.getRawPtr(), workSize, &info );
35  }
36 
37  template <class Scalar, class LocalOrdinal>
38  QR_InterfaceEx<Scalar,LocalOrdinal>::QR_InterfaceEx(const size_t NSDim) : workSize_(NSDim), NSDim_(NSDim), info_(0) {
39  tau_ = ArrayRCP<Scalar>(NSDim);
40  work_ = ArrayRCP<Scalar>(NSDim);
41  }
42 
43  template <class Scalar, class LocalOrdinal>
44  bool QR_InterfaceEx<Scalar,LocalOrdinal>::isZeroNspColumn(LocalOrdinal const &myAggSize, ArrayRCP<Scalar> &localQR, LocalOrdinal nspCol) {
45  bool bIsZeroColumn = true;
46  for(LocalOrdinal r = 0; r < myAggSize; ++r) {
47  if(localQR[ myAggSize*nspCol + r ]!=0.0) {
48  bIsZeroColumn = false;
49  break;
50  }
51  }
52  return bIsZeroColumn;
53  }
54 
55  template <class Scalar, class LocalOrdinal>
56  void QR_InterfaceEx<Scalar,LocalOrdinal>::Compute(LocalOrdinal const &myAggSize, ArrayRCP<Scalar> &localQR)
57  {
58  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
59  if (NSDim_ == 1) {
60  //only one nullspace vector, so normalize by hand
61  Magnitude dtemp=0;
62  //scalar type might be complex, so take absolute value.
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);
65  tau_[0] = localQR[0];
66  localQR[0] = dtemp;
67  } else {
68 
69  // check whether there are zero columns in localQR
70  nonZeroCols_.clear();
71  for(LocalOrdinal j = 0; j<NSDim_; ++j) {
72  if(!isZeroNspColumn(myAggSize, localQR, j))
73  nonZeroCols_.push_back(j);
74  }
75 
76  // set internal null space dimension used for QR decomposition
77  internalNSDim_ = Teuchos::as<LocalOrdinal>(nonZeroCols_.size());
78 
79  // reset internal variables
80  tau_ = ArrayRCP<Scalar>(internalNSDim_);
81  work_ = ArrayRCP<Scalar>(internalNSDim_);
83  info_ = 0;
84 
85  // build new input array
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];
90  }
91  }
92 
93  // print new input matrix
94  /*std::cout << " internal input matrix for GEQRF: " << std::endl;
95  for (int k=0; k<myAggSize; ++k) { // loop over rows
96  for (size_t j=0; j<internalNSDim_; ++j) { // loop over columns
97  std::cout << internalLocalQR_[j*myAggSize+k] << "\t";
98  }
99  std::cout << std::endl;
100  }*/
101 
102  // perform lapack call
103  lapack_.GEQRF( myAggSize, Teuchos::as<int>(internalNSDim_), internalLocalQR_.getRawPtr(), myAggSize,
104  tau_.getRawPtr(), work_.getRawPtr(), workSize_, &info_ );
105  if (info_ != 0) {
106  std::string msg = "QR_InterfaceEx: dgeqrf (LAPACK QR routine) returned error code " + Teuchos::toString(info_);
107  throw(Exceptions::RuntimeError(msg));
108  }
109 
110  // LAPACK may have determined a better length for the work array. Returns it in work[0],
111  // so we cast to avoid compiler warnings. Taking a look at the NETLIB reference implementation
112  // CGEQRF (complex), work[0] is assigned an integer, so it's safe to take the magnitude.
113  // Scalar type might be complex, so take absolute value.
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_);
117  }
118 
119  // write back results to localQR
120  // build new input array
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];
124  }
125  }
126 
127  // print new input matrix
128  /*std::cout << " localQR after internal GEQRF call: " << std::endl;
129  for (int k=0; k<myAggSize; ++k) { // loop over rows
130  for (size_t j=0; j<NSDim_; ++j) { // loop over columns
131  std::cout << localQR[j*myAggSize+k] << "\t";
132  }
133  std::cout << std::endl;
134  }*/
135  }
136  } //Compute()
137 
138  template <class Scalar, class LocalOrdinal>
139  void QR_InterfaceEx<Scalar,LocalOrdinal>::ExtractQ(LocalOrdinal const &myAggSize, ArrayRCP<Scalar> &localQR)
140  {
141  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
142  if (NSDim_ == 1) {
143  //again, only one nullspace vector, so calculate Q by hand
144  Magnitude dtemp = Teuchos::ScalarTraits<Scalar>::magnitude(localQR[0]);
145  localQR[0] = tau_[0];
146  dtemp = 1 / dtemp;
147  for (LocalOrdinal i=0; i<myAggSize; ++i)
148  localQR[i] *= dtemp;
149  } else {
150  // do LAPACK call on internalLocalQR_ array
151  //lapack_.ORGQR(myAggSize, Teuchos::as<int>(workSize_), Teuchos::as<int>(workSize_), localQR.getRawPtr(),
152  // myAggSize, tau_.getRawPtr(), work_.getRawPtr(), workSize_, &info_ );
153  //call nonmember function (perhaps specialized)
154  LapackQR( lapack_, myAggSize, Teuchos::as<int>(internalNSDim_), internalLocalQR_, tau_, work_, workSize_, info_ );
155  if (info_ != 0) {
156  std::string msg = "QR_InterfaceEx: dorgqr (LAPACK auxiliary QR routine) returned error code " + Teuchos::toString(info_);
157  throw(Exceptions::RuntimeError(msg));
158  }
159  // LAPACK may have determined a better length for the work array. Returns it in work[0],
160  // so we cast to avoid compiler warnings.
161  // Scalar type might be complex, so take absolute value.
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_);
165  }
166 
167  // write back results to localQR
168  // build new input array
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];
172  }
173  }
174 
175  // print new input matrix
176  /*std::cout << " localQR after internal ORGQR call: " << std::endl;
177  for (int k=0; k<myAggSize; ++k) { // loop over rows
178  for (size_t j=0; j<NSDim_; ++j) { // loop over columns
179  std::cout << localQR[j*myAggSize+k] << "\t";
180  }
181  std::cout << std::endl;
182  }*/
183  }
184  } //ExtractQ2()
185 
186 } //namespace MueLu
187 
188 #endif /* MUELU_QR_INTERFACEEX_DEF_HPP_ */
QR_InterfaceEx(const size_t nullSpaceDimension)
Constructor.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
LocalOrdinal info_
(out) =0: success; =i, i<0: i-th argument has illegal value
Teuchos::LAPACK< LocalOrdinal, Scalar > lapack_
Teuchos LAPACK wrapper.
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.
std::vector< LocalOrdinal > nonZeroCols_
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.
LocalOrdinal NSDim_
Dimension of nullspace.
bool isZeroNspColumn(LocalOrdinal const &myAggSize, ArrayRCP< Scalar > &localQR, LocalOrdinal nspCol)
Exception throws to report errors in the internal logical of the program.
LocalOrdinal workSize_
Length of work vectors. Must be at least dimension of nullspace.
ArrayRCP< Scalar > internalLocalQR_
internal local QR (without zero columns)