Mercurial > octave
annotate liboctave/dbleQRP.cc @ 10905:9ee5a0a1b93d
Return permutation vector from QR as a row, not column, vector.
Fixes Bug #30832.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Sat, 21 Aug 2010 13:02:34 -0700 |
parents | 4d1fc073fbb7 |
children | 141b3fb5cef7 |
rev | line source |
---|---|
538 | 1 /* |
2 | |
8920 | 3 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007, |
4 2008, 2009 John W. Eaton | |
10521
4d1fc073fbb7
add some missing copyright stmts
Jaroslav Hajek <highegg@gmail.com>
parents:
10350
diff
changeset
|
5 Copyright (C) 2009 VZLU Prague |
538 | 6 |
7 This file is part of Octave. | |
8 | |
9 Octave is free software; you can redistribute it and/or modify it | |
10 under the terms of the GNU General Public License as published by the | |
7016 | 11 Free Software Foundation; either version 3 of the License, or (at your |
12 option) any later version. | |
538 | 13 |
14 Octave is distributed in the hope that it will be useful, but WITHOUT | |
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
17 for more details. | |
18 | |
19 You should have received a copy of the GNU General Public License | |
7016 | 20 along with Octave; see the file COPYING. If not, see |
21 <http://www.gnu.org/licenses/>. | |
538 | 22 |
23 */ | |
24 | |
25 #ifdef HAVE_CONFIG_H | |
1192 | 26 #include <config.h> |
538 | 27 #endif |
28 | |
1367 | 29 #include <cassert> |
538 | 30 |
31 #include "dbleQRP.h" | |
1847 | 32 #include "f77-fcn.h" |
538 | 33 #include "lo-error.h" |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
34 #include "oct-locbuf.h" |
538 | 35 |
36 extern "C" | |
37 { | |
4552 | 38 F77_RET_T |
8368
c72c1c9bccdc
call blocked permuted qr factorization routines from LAPACK
Jaroslav Hajek <highegg@gmail.com>
parents:
8367
diff
changeset
|
39 F77_FUNC (dgeqp3, DGEQP3) (const octave_idx_type&, const octave_idx_type&, double*, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
40 const octave_idx_type&, 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 |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
69 MArray<octave_idx_type> jpvt (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; |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
75 F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (), |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
76 tau, &rlwork, -1, info)); |
8368
c72c1c9bccdc
call blocked permuted qr factorization routines from LAPACK
Jaroslav Hajek <highegg@gmail.com>
parents:
8367
diff
changeset
|
77 |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
78 // allocate buffer and do the job. |
8811 | 79 octave_idx_type lwork = rlwork; |
80 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
|
81 OCTAVE_LOCAL_BUFFER (double, work, lwork); |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
82 F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (), |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
83 tau, work, lwork, info)); |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
84 } |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
85 else |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
86 for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; |
1922 | 87 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
88 // 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
|
89 // indices only!) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
90 |
8811 | 91 jpvt -= static_cast<octave_idx_type> (1); |
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
92 p = PermMatrix (jpvt, true); |
1922 | 93 |
94 | |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
95 form (n, afact, tau, qr_type); |
538 | 96 } |
97 | |
10905
9ee5a0a1b93d
Return permutation vector from QR as a row, not column, vector.
Rik <octave@nomad.inbox5.com>
parents:
10521
diff
changeset
|
98 RowVector |
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
99 QRP::Pvec (void) const |
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
100 { |
8375
e3c9102431a9
fix design problems of diag & perm matrix classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8368
diff
changeset
|
101 Array<double> pa (p.pvec ()); |
10905
9ee5a0a1b93d
Return permutation vector from QR as a row, not column, vector.
Rik <octave@nomad.inbox5.com>
parents:
10521
diff
changeset
|
102 RowVector pv (MArray<double> (pa) + 1.0); |
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
103 return pv; |
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
104 } |