IT++ Logo
qr.cpp
Go to the documentation of this file.
1 
29 #ifndef _MSC_VER
30 # include <itpp/config.h>
31 #else
32 # include <itpp/config_msvc.h>
33 #endif
34 
35 #if defined(HAVE_LAPACK)
36 # include <itpp/base/algebra/lapack.h>
37 #endif
38 
39 #include <itpp/base/algebra/qr.h>
40 #include <itpp/base/specmat.h>
41 
42 
43 namespace itpp
44 {
45 
46 #if defined(HAVE_LAPACK)
47 
48 bool qr(const mat &A, mat &Q, mat &R)
49 {
50  int info;
51  int m = A.rows();
52  int n = A.cols();
53  int lwork = n;
54  int k = std::min(m, n);
55  vec tau(k);
56  vec work(lwork);
57 
58  R = A;
59 
60  // perform workspace query for optimum lwork value
61  int lwork_tmp = -1;
62  dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
63  &info);
64  if (info == 0) {
65  lwork = static_cast<int>(work(0));
66  work.set_size(lwork, false);
67  }
68  dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
69  Q = R;
70  Q.set_size(m, m, true);
71 
72  // construct R
73  for (int i = 0; i < m; i++)
74  for (int j = 0; j < std::min(i, n); j++)
75  R(i, j) = 0;
76 
77  // perform workspace query for optimum lwork value
78  lwork_tmp = -1;
79  dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp,
80  &info);
81  if (info == 0) {
82  lwork = static_cast<int>(work(0));
83  work.set_size(lwork, false);
84  }
85  dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork,
86  &info);
87 
88  return (info == 0);
89 }
90 
91 bool qr(const mat &A, mat &R)
92 {
93  int info;
94  int m = A.rows();
95  int n = A.cols();
96  int lwork = n;
97  int k = std::min(m, n);
98  vec tau(k);
99  vec work(lwork);
100 
101  R = A;
102 
103  // perform workspace query for optimum lwork value
104  int lwork_tmp = -1;
105  dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
106  &info);
107  if (info == 0) {
108  lwork = static_cast<int>(work(0));
109  work.set_size(lwork, false);
110  }
111  dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
112 
113  // construct R
114  for (int i = 0; i < m; i++)
115  for (int j = 0; j < std::min(i, n); j++)
116  R(i, j) = 0;
117 
118  return (info == 0);
119 }
120 
121 bool qr(const mat &A, mat &Q, mat &R, bmat &P)
122 {
123  int info;
124  int m = A.rows();
125  int n = A.cols();
126  int lwork = n;
127  int k = std::min(m, n);
128  vec tau(k);
129  vec work(lwork);
130  ivec jpvt(n);
131  jpvt.zeros();
132 
133  R = A;
134 
135  // perform workspace query for optimum lwork value
136  int lwork_tmp = -1;
137  dgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
138  &lwork_tmp, &info);
139  if (info == 0) {
140  lwork = static_cast<int>(work(0));
141  work.set_size(lwork, false);
142  }
143  dgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
144  &lwork, &info);
145  Q = R;
146  Q.set_size(m, m, true);
147 
148  // construct permutation matrix
149  P = zeros_b(n, n);
150  for (int j = 0; j < n; j++)
151  P(jpvt(j) - 1, j) = 1;
152 
153  // construct R
154  for (int i = 0; i < m; i++)
155  for (int j = 0; j < std::min(i, n); j++)
156  R(i, j) = 0;
157 
158  // perform workspace query for optimum lwork value
159  lwork_tmp = -1;
160  dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp,
161  &info);
162  if (info == 0) {
163  lwork = static_cast<int>(work(0));
164  work.set_size(lwork, false);
165  }
166  dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork,
167  &info);
168 
169  return (info == 0);
170 }
171 
172 
173 
174 bool qr(const cmat &A, cmat &Q, cmat &R)
175 {
176  int info;
177  int m = A.rows();
178  int n = A.cols();
179  int lwork = n;
180  int k = std::min(m, n);
181  cvec tau(k);
182  cvec work(lwork);
183 
184  R = A;
185 
186  // perform workspace query for optimum lwork value
187  int lwork_tmp = -1;
188  zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
189  &info);
190  if (info == 0) {
191  lwork = static_cast<int>(real(work(0)));
192  work.set_size(lwork, false);
193  }
194  zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
195 
196  Q = R;
197  Q.set_size(m, m, true);
198 
199  // construct R
200  for (int i = 0; i < m; i++)
201  for (int j = 0; j < std::min(i, n); j++)
202  R(i, j) = 0;
203 
204  // perform workspace query for optimum lwork value
205  lwork_tmp = -1;
206  zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp,
207  &info);
208  if (info == 0) {
209  lwork = static_cast<int>(real(work(0)));
210  work.set_size(lwork, false);
211  }
212  zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork,
213  &info);
214 
215  return (info == 0);
216 }
217 
218 bool qr(const cmat &A, cmat &R)
219 {
220  int info;
221  int m = A.rows();
222  int n = A.cols();
223  int lwork = n;
224  int k = std::min(m, n);
225  cvec tau(k);
226  cvec work(lwork);
227 
228  R = A;
229 
230  // perform workspace query for optimum lwork value
231  int lwork_tmp = -1;
232  zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
233  &info);
234  if (info == 0) {
235  lwork = static_cast<int>(real(work(0)));
236  work.set_size(lwork, false);
237  }
238  zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
239 
240  // construct R
241  for (int i = 0; i < m; i++)
242  for (int j = 0; j < std::min(i, n); j++)
243  R(i, j) = 0;
244 
245  return (info == 0);
246 }
247 
248 bool qr(const cmat &A, cmat &Q, cmat &R, bmat &P)
249 {
250  int info;
251  int m = A.rows();
252  int n = A.cols();
253  int lwork = n;
254  int k = std::min(m, n);
255  cvec tau(k);
256  cvec work(lwork);
257  vec rwork(std::max(1, 2*n));
258  ivec jpvt(n);
259  jpvt.zeros();
260 
261  R = A;
262 
263  // perform workspace query for optimum lwork value
264  int lwork_tmp = -1;
265  zgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
266  &lwork_tmp, rwork._data(), &info);
267  if (info == 0) {
268  lwork = static_cast<int>(real(work(0)));
269  work.set_size(lwork, false);
270  }
271  zgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
272  &lwork, rwork._data(), &info);
273 
274  Q = R;
275  Q.set_size(m, m, true);
276 
277  // construct permutation matrix
278  P = zeros_b(n, n);
279  for (int j = 0; j < n; j++)
280  P(jpvt(j) - 1, j) = 1;
281 
282  // construct R
283  for (int i = 0; i < m; i++)
284  for (int j = 0; j < std::min(i, n); j++)
285  R(i, j) = 0;
286 
287  // perform workspace query for optimum lwork value
288  lwork_tmp = -1;
289  zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp,
290  &info);
291  if (info == 0) {
292  lwork = static_cast<int>(real(work(0)));
293  work.set_size(lwork, false);
294  }
295  zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork,
296  &info);
297 
298  return (info == 0);
299 }
300 
301 #else
302 
303 bool qr(const mat &A, mat &Q, mat &R)
304 {
305  it_error("LAPACK library is needed to use qr() function");
306  return false;
307 }
308 
309 bool qr(const mat &A, mat &R)
310 {
311  it_error("LAPACK library is needed to use qr() function");
312  return false;
313 }
314 
315 bool qr(const mat &A, mat &Q, mat &R, bmat &P)
316 {
317  it_error("LAPACK library is needed to use qr() function");
318  return false;
319 }
320 
321 bool qr(const cmat &A, cmat &Q, cmat &R)
322 {
323  it_error("LAPACK library is needed to use qr() function");
324  return false;
325 }
326 
327 bool qr(const cmat &A, cmat &R)
328 {
329  it_error("LAPACK library is needed to use qr() function");
330  return false;
331 }
332 
333 bool qr(const cmat &A, cmat &Q, cmat &R, bmat &P)
334 {
335  it_error("LAPACK library is needed to use qr() function");
336  return false;
337 }
338 
339 #endif // HAVE_LAPACK
340 
341 } // namespace itpp
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:19 for IT++ by Doxygen 1.8.2