Mercurial > octave
comparison liboctave/numeric/qrp.cc @ 21280:ebdf74c15722
better use of templates for qrp classes
* liboctave/numeric/qrp.h, liboctave/numeric/qrp.cc: New files for qrp
classes generated from CmplxQRP.cc, CmplxQRP.h, dbleQRP.cc, dbleQRP.h,
fCmplxQRP.cc, fCmplxQRP.h, floatQRP.cc, and floatQRP.h with classes
converted to templates.
* liboctave/numeric/module.mk: Update.
* qr.cc, mx-defs.h: Use new classes.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 17 Feb 2016 02:54:00 -0500 |
parents | liboctave/numeric/dbleQRP.cc@eb1524b07fe3 |
children | 40de9f8f23a6 |
comparison
equal
deleted
inserted
replaced
21279:eb1524b07fe3 | 21280:ebdf74c15722 |
---|---|
1 /* | |
2 | |
3 Copyright (C) 1994-2015 John W. Eatonn | |
4 Copyright (C) 2009 VZLU Prague | |
5 | |
6 This file is part of Octave. | |
7 | |
8 Octave is free software; you can redistribute it and/or modify it | |
9 under the terms of the GNU General Public License as published by the | |
10 Free Software Foundation; either version 3 of the License, or (at your | |
11 option) any later version. | |
12 | |
13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
16 for more details. | |
17 | |
18 You should have received a copy of the GNU General Public License | |
19 along with Octave; see the file COPYING. If not, see | |
20 <http://www.gnu.org/licenses/>. | |
21 | |
22 */ | |
23 | |
24 #ifdef HAVE_CONFIG_H | |
25 # include <config.h> | |
26 #endif | |
27 | |
28 #include <cassert> | |
29 | |
30 #include "CMatrix.h" | |
31 #include "dMatrix.h" | |
32 #include "dRowVector.h" | |
33 #include "f77-fcn.h" | |
34 #include "fCMatrix.h" | |
35 #include "fMatrix.h" | |
36 #include "fRowVector.h" | |
37 #include "lo-error.h" | |
38 #include "oct-locbuf.h" | |
39 #include "qrp.h" | |
40 | |
41 extern "C" | |
42 { | |
43 F77_RET_T | |
44 F77_FUNC (dgeqp3, DGEQP3) (const octave_idx_type&, const octave_idx_type&, | |
45 double*, const octave_idx_type&, | |
46 octave_idx_type*, double*, double*, | |
47 const octave_idx_type&, octave_idx_type&); | |
48 | |
49 F77_RET_T | |
50 F77_FUNC (sgeqp3, SGEQP3) (const octave_idx_type&, const octave_idx_type&, | |
51 float*, const octave_idx_type&, octave_idx_type*, | |
52 float*, float*, const octave_idx_type&, | |
53 octave_idx_type&); | |
54 F77_RET_T | |
55 F77_FUNC (zgeqp3, ZGEQP3) (const octave_idx_type&, const octave_idx_type&, | |
56 Complex*, const octave_idx_type&, | |
57 octave_idx_type*, Complex*, Complex*, | |
58 const octave_idx_type&, double*, | |
59 octave_idx_type&); | |
60 F77_RET_T | |
61 F77_FUNC (cgeqp3, CGEQP3) (const octave_idx_type&, const octave_idx_type&, | |
62 FloatComplex*, const octave_idx_type&, | |
63 octave_idx_type*, FloatComplex*, FloatComplex*, | |
64 const octave_idx_type&, float*, octave_idx_type&); | |
65 } | |
66 | |
67 // Specialization. | |
68 | |
69 template <> | |
70 void | |
71 qrp<Matrix>::init (const Matrix& a, type qr_type) | |
72 { | |
73 assert (qr_type != qr<Matrix>::raw); | |
74 | |
75 octave_idx_type m = a.rows (); | |
76 octave_idx_type n = a.cols (); | |
77 | |
78 octave_idx_type min_mn = m < n ? m : n; | |
79 OCTAVE_LOCAL_BUFFER (double, tau, min_mn); | |
80 | |
81 octave_idx_type info = 0; | |
82 | |
83 Matrix afact = a; | |
84 if (m > n && qr_type == qr<Matrix>::std) | |
85 afact.resize (m, m); | |
86 | |
87 MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); | |
88 | |
89 if (m > 0) | |
90 { | |
91 // workspace query. | |
92 double rlwork; | |
93 F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), | |
94 m, jpvt.fortran_vec (), tau, | |
95 &rlwork, -1, info)); | |
96 | |
97 // allocate buffer and do the job. | |
98 octave_idx_type lwork = rlwork; | |
99 lwork = std::max (lwork, static_cast<octave_idx_type> (1)); | |
100 OCTAVE_LOCAL_BUFFER (double, work, lwork); | |
101 F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), | |
102 m, jpvt.fortran_vec (), tau, | |
103 work, lwork, info)); | |
104 } | |
105 else | |
106 for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; | |
107 | |
108 // Form Permutation matrix (if economy is requested, return the | |
109 // indices only!) | |
110 | |
111 jpvt -= static_cast<octave_idx_type> (1); | |
112 p = PermMatrix (jpvt, true); | |
113 | |
114 | |
115 form (n, afact, tau, qr_type); | |
116 } | |
117 | |
118 template <> | |
119 qrp<Matrix>::qrp (const Matrix& a, type qr_type) | |
120 : qr<Matrix> (), p () | |
121 { | |
122 init (a, qr_type); | |
123 } | |
124 | |
125 template <> | |
126 RowVector | |
127 qrp<Matrix>::Pvec (void) const | |
128 { | |
129 Array<double> pa (p.col_perm_vec ()); | |
130 RowVector pv (MArray<double> (pa) + 1.0); | |
131 return pv; | |
132 } | |
133 | |
134 template <> | |
135 void | |
136 qrp<FloatMatrix>::init (const FloatMatrix& a, type qr_type) | |
137 { | |
138 assert (qr_type != qr<FloatMatrix>::raw); | |
139 | |
140 octave_idx_type m = a.rows (); | |
141 octave_idx_type n = a.cols (); | |
142 | |
143 octave_idx_type min_mn = m < n ? m : n; | |
144 OCTAVE_LOCAL_BUFFER (float, tau, min_mn); | |
145 | |
146 octave_idx_type info = 0; | |
147 | |
148 FloatMatrix afact = a; | |
149 if (m > n && qr_type == qr<FloatMatrix>::std) | |
150 afact.resize (m, m); | |
151 | |
152 MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); | |
153 | |
154 if (m > 0) | |
155 { | |
156 // workspace query. | |
157 float rlwork; | |
158 F77_XFCN (sgeqp3, SGEQP3, (m, n, afact.fortran_vec (), | |
159 m, jpvt.fortran_vec (), tau, | |
160 &rlwork, -1, info)); | |
161 | |
162 // allocate buffer and do the job. | |
163 octave_idx_type lwork = rlwork; | |
164 lwork = std::max (lwork, static_cast<octave_idx_type> (1)); | |
165 OCTAVE_LOCAL_BUFFER (float, work, lwork); | |
166 F77_XFCN (sgeqp3, SGEQP3, (m, n, afact.fortran_vec (), | |
167 m, jpvt.fortran_vec (), tau, | |
168 work, lwork, info)); | |
169 } | |
170 else | |
171 for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; | |
172 | |
173 // Form Permutation matrix (if economy is requested, return the | |
174 // indices only!) | |
175 | |
176 jpvt -= static_cast<octave_idx_type> (1); | |
177 p = PermMatrix (jpvt, true); | |
178 | |
179 | |
180 form (n, afact, tau, qr_type); | |
181 } | |
182 | |
183 template <> | |
184 qrp<FloatMatrix>::qrp (const FloatMatrix& a, type qr_type) | |
185 : qr<FloatMatrix> (), p () | |
186 { | |
187 init (a, qr_type); | |
188 } | |
189 | |
190 template <> | |
191 FloatRowVector | |
192 qrp<FloatMatrix>::Pvec (void) const | |
193 { | |
194 Array<float> pa (p.col_perm_vec ()); | |
195 FloatRowVector pv (MArray<float> (pa) + 1.0f); | |
196 return pv; | |
197 } | |
198 | |
199 template <> | |
200 void | |
201 qrp<ComplexMatrix>::init (const ComplexMatrix& a, type qr_type) | |
202 { | |
203 assert (qr_type != qr<ComplexMatrix>::raw); | |
204 | |
205 octave_idx_type m = a.rows (); | |
206 octave_idx_type n = a.cols (); | |
207 | |
208 octave_idx_type min_mn = m < n ? m : n; | |
209 OCTAVE_LOCAL_BUFFER (Complex, tau, min_mn); | |
210 | |
211 octave_idx_type info = 0; | |
212 | |
213 ComplexMatrix afact = a; | |
214 if (m > n && qr_type == qr<ComplexMatrix>::std) | |
215 afact.resize (m, m); | |
216 | |
217 MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); | |
218 | |
219 if (m > 0) | |
220 { | |
221 OCTAVE_LOCAL_BUFFER (double, rwork, 2*n); | |
222 | |
223 // workspace query. | |
224 Complex clwork; | |
225 F77_XFCN (zgeqp3, ZGEQP3, (m, n, afact.fortran_vec (), | |
226 m, jpvt.fortran_vec (), tau, | |
227 &clwork, -1, rwork, info)); | |
228 | |
229 // allocate buffer and do the job. | |
230 octave_idx_type lwork = clwork.real (); | |
231 lwork = std::max (lwork, static_cast<octave_idx_type> (1)); | |
232 OCTAVE_LOCAL_BUFFER (Complex, work, lwork); | |
233 F77_XFCN (zgeqp3, ZGEQP3, (m, n, afact.fortran_vec (), | |
234 m, jpvt.fortran_vec (), tau, | |
235 work, lwork, rwork, info)); | |
236 } | |
237 else | |
238 for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; | |
239 | |
240 // Form Permutation matrix (if economy is requested, return the | |
241 // indices only!) | |
242 | |
243 jpvt -= static_cast<octave_idx_type> (1); | |
244 p = PermMatrix (jpvt, true); | |
245 | |
246 | |
247 form (n, afact, tau, qr_type); | |
248 } | |
249 | |
250 template <> | |
251 qrp<ComplexMatrix>::qrp (const ComplexMatrix& a, type qr_type) | |
252 : qr<ComplexMatrix> (), p () | |
253 { | |
254 init (a, qr_type); | |
255 } | |
256 | |
257 template <> | |
258 RowVector | |
259 qrp<ComplexMatrix>::Pvec (void) const | |
260 { | |
261 Array<double> pa (p.col_perm_vec ()); | |
262 RowVector pv (MArray<double> (pa) + 1.0); | |
263 return pv; | |
264 } | |
265 | |
266 template <> | |
267 void | |
268 qrp<FloatComplexMatrix>::init (const FloatComplexMatrix& a, type qr_type) | |
269 { | |
270 assert (qr_type != qr<FloatComplexMatrix>::raw); | |
271 | |
272 octave_idx_type m = a.rows (); | |
273 octave_idx_type n = a.cols (); | |
274 | |
275 octave_idx_type min_mn = m < n ? m : n; | |
276 OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn); | |
277 | |
278 octave_idx_type info = 0; | |
279 | |
280 FloatComplexMatrix afact = a; | |
281 if (m > n && qr_type == qr<FloatComplexMatrix>::std) | |
282 afact.resize (m, m); | |
283 | |
284 MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); | |
285 | |
286 if (m > 0) | |
287 { | |
288 OCTAVE_LOCAL_BUFFER (float, rwork, 2*n); | |
289 | |
290 // workspace query. | |
291 FloatComplex clwork; | |
292 F77_XFCN (cgeqp3, CGEQP3, (m, n, afact.fortran_vec (), | |
293 m, jpvt.fortran_vec (), tau, | |
294 &clwork, -1, rwork, info)); | |
295 | |
296 // allocate buffer and do the job. | |
297 octave_idx_type lwork = clwork.real (); | |
298 lwork = std::max (lwork, static_cast<octave_idx_type> (1)); | |
299 OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork); | |
300 F77_XFCN (cgeqp3, CGEQP3, (m, n, afact.fortran_vec (), | |
301 m, jpvt.fortran_vec (), tau, | |
302 work, lwork, rwork, info)); | |
303 } | |
304 else | |
305 for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; | |
306 | |
307 // Form Permutation matrix (if economy is requested, return the | |
308 // indices only!) | |
309 | |
310 jpvt -= static_cast<octave_idx_type> (1); | |
311 p = PermMatrix (jpvt, true); | |
312 | |
313 | |
314 form (n, afact, tau, qr_type); | |
315 } | |
316 | |
317 template <> | |
318 qrp<FloatComplexMatrix>::qrp (const FloatComplexMatrix& a, type qr_type) | |
319 : qr<FloatComplexMatrix> (), p () | |
320 { | |
321 init (a, qr_type); | |
322 } | |
323 | |
324 template <> | |
325 FloatRowVector | |
326 qrp<FloatComplexMatrix>::Pvec (void) const | |
327 { | |
328 Array<float> pa (p.col_perm_vec ()); | |
329 FloatRowVector pv (MArray<float> (pa) + 1.0f); | |
330 return pv; | |
331 } |