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
40#include <itpp/base/specmat.h>
41
42
43namespace itpp
44{
45
46#if defined(HAVE_LAPACK)
47
48bool 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
91bool 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
121bool 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
174bool 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
218bool 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
248bool 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
303bool 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
309bool qr(const mat &A, mat &R)
310{
311 it_error("LAPACK library is needed to use qr() function");
312 return false;
313}
314
315bool 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
321bool 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
327bool qr(const cmat &A, cmat &R)
328{
329 it_error("LAPACK library is needed to use qr() function");
330 return false;
331}
332
333bool 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
#define it_error(s)
Abort unconditionally.
Definition: itassert.h:126
bool qr(const mat &A, mat &Q, mat &R)
QR factorisation of real matrix.
Definition: qr.cpp:303
vec real(const cvec &data)
Real part of complex values.
Definition: elem_math.cpp:157
ITPP_EXPORT bvec zeros_b(int size)
A Binary vector of zeros.
Mat< bin > bmat
bin matrix
Definition: mat.h:508
itpp namespace
Definition: itmex.h:37
Definitions of QR factorisation functions.
Definitions of special vectors and matrices.
SourceForge Logo

Generated on Tue Jan 24 2023 00:00:00 for IT++ by Doxygen 1.9.6