GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
qrp.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1994-2024 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include <cassert>
31 
32 #include <algorithm>
33 
34 #include "Array.h"
35 #include "CMatrix.h"
36 #include "MArray.h"
37 #include "dMatrix.h"
38 #include "dRowVector.h"
39 #include "fCMatrix.h"
40 #include "fMatrix.h"
41 #include "fRowVector.h"
42 #include "lo-lapack-proto.h"
43 #include "oct-locbuf.h"
44 #include "qrp.h"
45 
47 
49 
50 // Specialization.
51 
52 template <>
54 void
55 qrp<Matrix>::init (const Matrix& a, type qr_type)
56 {
57  assert (qr_type != qr<Matrix>::raw);
58 
59  F77_INT m = to_f77_int (a.rows ());
60  F77_INT n = to_f77_int (a.cols ());
61 
62  F77_INT min_mn = (m < n ? m : n);
63  OCTAVE_LOCAL_BUFFER (double, tau, min_mn);
64 
65  F77_INT info = 0;
66 
67  Matrix afact = a;
68  if (m > n && qr_type == qr<Matrix>::std)
69  afact.resize (m, m);
70 
71  MArray<F77_INT> jpvt (dim_vector (n, 1), 0);
72 
73  if (m > 0)
74  {
75  // workspace query.
76  double rlwork;
77  F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (),
78  m, jpvt.fortran_vec (), tau,
79  &rlwork, -1, info));
80 
81  // allocate buffer and do the job.
82  F77_INT lwork = static_cast<F77_INT> (rlwork);
83  lwork = std::max (lwork, static_cast<F77_INT> (1));
84  OCTAVE_LOCAL_BUFFER (double, work, lwork);
85 
86  F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (),
87  m, jpvt.fortran_vec (), tau,
88  work, lwork, info));
89  }
90  else
91  {
92  for (F77_INT i = 0; i < n; i++)
93  jpvt(i) = i+1;
94  }
95 
96  // Form Permutation matrix (if economy is requested, return the
97  // indices only!)
98 
99  jpvt -= static_cast<F77_INT> (1);
100  m_p = PermMatrix (jpvt, true);
101 
102  form (n, afact, tau, qr_type);
103 }
104 
105 template <>
107 qrp<Matrix>::qrp (const Matrix& a, type qr_type)
108  : qr<Matrix> (), m_p ()
109 {
110  init (a, qr_type);
111 }
112 
113 template <>
115 RowVector
117 {
118  Array<double> pa (m_p.col_perm_vec ());
119  RowVector pv (MArray<double> (pa) + 1.0);
120  return pv;
121 }
122 
123 template <>
125 void
127 {
128  assert (qr_type != qr<FloatMatrix>::raw);
129 
130  F77_INT m = to_f77_int (a.rows ());
131  F77_INT n = to_f77_int (a.cols ());
132 
133  F77_INT min_mn = (m < n ? m : n);
134  OCTAVE_LOCAL_BUFFER (float, tau, min_mn);
135 
136  F77_INT info = 0;
137 
138  FloatMatrix afact = a;
139  if (m > n && qr_type == qr<FloatMatrix>::std)
140  afact.resize (m, m);
141 
142  MArray<F77_INT> jpvt (dim_vector (n, 1), 0);
143 
144  if (m > 0)
145  {
146  // workspace query.
147  float rlwork;
148  F77_XFCN (sgeqp3, SGEQP3, (m, n, afact.fortran_vec (),
149  m, jpvt.fortran_vec (), tau,
150  &rlwork, -1, info));
151 
152  // allocate buffer and do the job.
153  F77_INT lwork = static_cast<F77_INT> (rlwork);
154  lwork = std::max (lwork, static_cast<F77_INT> (1));
155  OCTAVE_LOCAL_BUFFER (float, work, lwork);
156 
157  F77_XFCN (sgeqp3, SGEQP3, (m, n, afact.fortran_vec (),
158  m, jpvt.fortran_vec (), tau,
159  work, lwork, info));
160  }
161  else
162  {
163  for (F77_INT i = 0; i < n; i++)
164  jpvt(i) = i+1;
165  }
166 
167  // Form Permutation matrix (if economy is requested, return the
168  // indices only!)
169 
170  jpvt -= static_cast<F77_INT> (1);
171  m_p = PermMatrix (jpvt, true);
172 
173  form (n, afact, tau, qr_type);
174 }
175 
176 template <>
179  : qr<FloatMatrix> (), m_p ()
180 {
181  init (a, qr_type);
182 }
183 
184 template <>
188 {
189  Array<float> pa (m_p.col_perm_vec ());
190  FloatRowVector pv (MArray<float> (pa) + 1.0f);
191  return pv;
192 }
193 
194 template <>
196 void
198 {
199  assert (qr_type != qr<ComplexMatrix>::raw);
200 
201  F77_INT m = to_f77_int (a.rows ());
202  F77_INT n = to_f77_int (a.cols ());
203 
204  F77_INT min_mn = (m < n ? m : n);
205  OCTAVE_LOCAL_BUFFER (Complex, tau, min_mn);
206 
207  F77_INT info = 0;
208 
209  ComplexMatrix afact = a;
210  if (m > n && qr_type == qr<ComplexMatrix>::std)
211  afact.resize (m, m);
212 
213  MArray<F77_INT> jpvt (dim_vector (n, 1), 0);
214 
215  if (m > 0)
216  {
217  OCTAVE_LOCAL_BUFFER (double, rwork, 2*n);
218 
219  // workspace query.
220  Complex clwork;
221  F77_XFCN (zgeqp3, ZGEQP3, (m, n,
222  F77_DBLE_CMPLX_ARG (afact.fortran_vec ()),
223  m, jpvt.fortran_vec (),
224  F77_DBLE_CMPLX_ARG (tau),
225  F77_DBLE_CMPLX_ARG (&clwork),
226  -1, rwork, info));
227 
228  // allocate buffer and do the job.
229  F77_INT lwork = static_cast<F77_INT> (clwork.real ());
230  lwork = std::max (lwork, static_cast<F77_INT> (1));
231  OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
232 
233  F77_XFCN (zgeqp3, ZGEQP3, (m, n,
234  F77_DBLE_CMPLX_ARG (afact.fortran_vec ()),
235  m, jpvt.fortran_vec (),
236  F77_DBLE_CMPLX_ARG (tau),
237  F77_DBLE_CMPLX_ARG (work),
238  lwork, rwork, info));
239  }
240  else
241  {
242  for (F77_INT i = 0; i < n; i++)
243  jpvt(i) = i+1;
244  }
245 
246  // Form Permutation matrix (if economy is requested, return the
247  // indices only!)
248 
249  jpvt -= static_cast<F77_INT> (1);
250  m_p = PermMatrix (jpvt, true);
251 
252  form (n, afact, tau, qr_type);
253 }
254 
255 template <>
258  : qr<ComplexMatrix> (), m_p ()
259 {
260  init (a, qr_type);
261 }
262 
263 template <>
265 RowVector
267 {
268  Array<double> pa (m_p.col_perm_vec ());
269  RowVector pv (MArray<double> (pa) + 1.0);
270  return pv;
271 }
272 
273 template <>
275 void
277 {
278  assert (qr_type != qr<FloatComplexMatrix>::raw);
279 
280  F77_INT m = to_f77_int (a.rows ());
281  F77_INT n = to_f77_int (a.cols ());
282 
283  F77_INT min_mn = (m < n ? m : n);
284  OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn);
285 
286  F77_INT info = 0;
287 
288  FloatComplexMatrix afact = a;
289  if (m > n && qr_type == qr<FloatComplexMatrix>::std)
290  afact.resize (m, m);
291 
292  MArray<F77_INT> jpvt (dim_vector (n, 1), 0);
293 
294  if (m > 0)
295  {
296  OCTAVE_LOCAL_BUFFER (float, rwork, 2*n);
297 
298  // workspace query.
299  FloatComplex clwork;
300  F77_XFCN (cgeqp3, CGEQP3, (m, n,
301  F77_CMPLX_ARG (afact.fortran_vec ()),
302  m, jpvt.fortran_vec (),
303  F77_CMPLX_ARG (tau),
304  F77_CMPLX_ARG (&clwork),
305  -1, rwork, info));
306 
307  // allocate buffer and do the job.
308  F77_INT lwork = static_cast<F77_INT> (clwork.real ());
309  lwork = std::max (lwork, static_cast<F77_INT> (1));
310  OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork);
311 
312  F77_XFCN (cgeqp3, CGEQP3, (m, n,
313  F77_CMPLX_ARG (afact.fortran_vec ()),
314  m, jpvt.fortran_vec (),
315  F77_CMPLX_ARG (tau),
316  F77_CMPLX_ARG (work),
317  lwork, rwork, info));
318  }
319  else
320  {
321  for (F77_INT i = 0; i < n; i++)
322  jpvt(i) = i+1;
323  }
324 
325  // Form Permutation matrix (if economy is requested, return the
326  // indices only!)
327 
328  jpvt -= static_cast<F77_INT> (1);
329  m_p = PermMatrix (jpvt, true);
330 
331  form (n, afact, tau, qr_type);
332 }
333 
334 template <>
337  : qr<FloatComplexMatrix> (), m_p ()
338 {
339  init (a, qr_type);
340 }
341 
342 template <>
346 {
347  Array<float> pa (m_p.col_perm_vec ());
348  FloatRowVector pv (MArray<float> (pa) + 1.0f);
349  return pv;
350 }
351 
352 OCTAVE_END_NAMESPACE(math)
353 OCTAVE_END_NAMESPACE(octave)
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:230
T * fortran_vec()
Size of the specified dimension.
Definition: Array-base.cc:1764
octave_idx_type rows() const
Definition: Array.h:459
octave_idx_type cols() const
Definition: Array.h:469
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Definition: CMatrix.h:193
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
Definition: fCMatrix.h:201
void resize(octave_idx_type nr, octave_idx_type nc, float rfv=0)
Definition: fMatrix.h:158
Template for N-dimensional array classes with like-type math operators.
Definition: MArray.h:63
Definition: dMatrix.h:42
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:158
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
Definition: qr.h:40
type
Definition: qr.h:48
void init(const T &, type=qr< T >::std)
RV_T Pvec() const
qrp()
Definition: qrp.h:48
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:316
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:310
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:45
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
#define OCTAVE_API
Definition: main.cc:55
T octave_idx_type m
Definition: mx-inlines.cc:781
octave_idx_type n
Definition: mx-inlines.cc:761
std::complex< double > Complex
Definition: oct-cmplx.h:33
std::complex< float > FloatComplex
Definition: oct-cmplx.h:34
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:44