Mercurial > hg > octave-nkf
annotate liboctave/numeric/dbleQRP.cc @ 20761:b7ac1e94266e
maint: Further clean up of functions in ode/private dir.
* AbsRel_Norm.m, fuzzy_compare.m, integrate_adaptive.m, integrate_const.m,
integrate_n_steps.m, ode_struct_value_check.m, odepkg_event_handle.m,
odepkg_structure_check.m, runge_kutta_45_dorpri.m:
Place latest copyright first in file.
Use two spaces before beginning single-line comment.
Use parentheses around variable to be tested in switch stmt.
Use space between function name and opening parenthesis.
author | Rik <rik@octave.org> |
---|---|
date | Mon, 05 Oct 2015 12:03:16 -0700 |
parents | 4197fc428c7d |
children |
rev | line source |
---|---|
538 | 1 /* |
2 | |
19898
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
19039
diff
changeset
|
3 Copyright (C) 1994-2015 John W. Eaton |
10521
4d1fc073fbb7
add some missing copyright stmts
Jaroslav Hajek <highegg@gmail.com>
parents:
10350
diff
changeset
|
4 Copyright (C) 2009 VZLU Prague |
538 | 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 | |
7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
11 option) any later version. | |
538 | 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 | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
538 | 21 |
22 */ | |
23 | |
24 #ifdef HAVE_CONFIG_H | |
1192 | 25 #include <config.h> |
538 | 26 #endif |
27 | |
1367 | 28 #include <cassert> |
538 | 29 |
30 #include "dbleQRP.h" | |
1847 | 31 #include "f77-fcn.h" |
538 | 32 #include "lo-error.h" |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
33 #include "oct-locbuf.h" |
538 | 34 |
35 extern "C" | |
36 { | |
4552 | 37 F77_RET_T |
11518 | 38 F77_FUNC (dgeqp3, DGEQP3) (const octave_idx_type&, const octave_idx_type&, |
39 double*, const octave_idx_type&, | |
40 octave_idx_type*, double*, double*, | |
8368
c72c1c9bccdc
call blocked permuted qr factorization routines from LAPACK
Jaroslav Hajek <highegg@gmail.com>
parents:
8367
diff
changeset
|
41 const octave_idx_type&, octave_idx_type&); |
538 | 42 } |
43 | |
44 // It would be best to share some of this code with QR class... | |
45 | |
9713
7918eb15040c
refactor the QR classes onto a templated base
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
46 QRP::QRP (const Matrix& a, qr_type_t qr_type) |
2763 | 47 : QR (), p () |
48 { | |
49 init (a, qr_type); | |
50 } | |
51 | |
52 void | |
9713
7918eb15040c
refactor the QR classes onto a templated base
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
53 QRP::init (const Matrix& a, qr_type_t qr_type) |
538 | 54 { |
9713
7918eb15040c
refactor the QR classes onto a templated base
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
55 assert (qr_type != qr_type_raw); |
538 | 56 |
5275 | 57 octave_idx_type m = a.rows (); |
58 octave_idx_type n = a.cols (); | |
538 | 59 |
5275 | 60 octave_idx_type min_mn = m < n ? m : n; |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
61 OCTAVE_LOCAL_BUFFER (double, tau, min_mn); |
1922 | 62 |
5275 | 63 octave_idx_type info = 0; |
538 | 64 |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
65 Matrix afact = a; |
9713
7918eb15040c
refactor the QR classes onto a templated base
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
66 if (m > n && qr_type == qr_type_std) |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
67 afact.resize (m, m); |
538 | 68 |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
69 MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); |
538 | 70 |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
71 if (m > 0) |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
72 { |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
73 // workspace query. |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
74 double rlwork; |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
75 F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
76 m, jpvt.fortran_vec (), tau, |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
77 &rlwork, -1, info)); |
8368
c72c1c9bccdc
call blocked permuted qr factorization routines from LAPACK
Jaroslav Hajek <highegg@gmail.com>
parents:
8367
diff
changeset
|
78 |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
79 // allocate buffer and do the job. |
8811 | 80 octave_idx_type lwork = rlwork; |
81 lwork = std::max (lwork, static_cast<octave_idx_type> (1)); | |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
82 OCTAVE_LOCAL_BUFFER (double, work, lwork); |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
83 F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
84 m, jpvt.fortran_vec (), tau, |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
85 work, lwork, info)); |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
86 } |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
87 else |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
88 for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; |
1922 | 89 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
90 // Form Permutation matrix (if economy is requested, return the |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
91 // indices only!) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
92 |
8811 | 93 jpvt -= static_cast<octave_idx_type> (1); |
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
94 p = PermMatrix (jpvt, true); |
1922 | 95 |
96 | |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
97 form (n, afact, tau, qr_type); |
538 | 98 } |
99 | |
10905
9ee5a0a1b93d
Return permutation vector from QR as a row, not column, vector.
Rik <octave@nomad.inbox5.com>
parents:
10521
diff
changeset
|
100 RowVector |
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
101 QRP::Pvec (void) const |
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
102 { |
19039
aa9ca67f09fb
make all permutation matrices column permutations (bug #42418)
David Spies <dnspies@gmail.com>
parents:
17769
diff
changeset
|
103 Array<double> pa (p.col_perm_vec ()); |
10905
9ee5a0a1b93d
Return permutation vector from QR as a row, not column, vector.
Rik <octave@nomad.inbox5.com>
parents:
10521
diff
changeset
|
104 RowVector pv (MArray<double> (pa) + 1.0); |
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
105 return pv; |
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
106 } |