Mercurial > octave-nkf
annotate liboctave/dbleQRP.cc @ 11518:141b3fb5cef7
style fixes
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Thu, 13 Jan 2011 16:52:30 -0500 |
parents | 9ee5a0a1b93d |
children | fd0a3ac60b0e |
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 |
11518 | 39 F77_FUNC (dgeqp3, DGEQP3) (const octave_idx_type&, const octave_idx_type&, |
40 double*, const octave_idx_type&, | |
41 octave_idx_type*, double*, double*, | |
8368
c72c1c9bccdc
call blocked permuted qr factorization routines from LAPACK
Jaroslav Hajek <highegg@gmail.com>
parents:
8367
diff
changeset
|
42 const octave_idx_type&, octave_idx_type&); |
538 | 43 } |
44 | |
45 // It would be best to share some of this code with QR class... | |
46 | |
9713
7918eb15040c
refactor the QR classes onto a templated base
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
47 QRP::QRP (const Matrix& a, qr_type_t qr_type) |
2763 | 48 : QR (), p () |
49 { | |
50 init (a, qr_type); | |
51 } | |
52 | |
53 void | |
9713
7918eb15040c
refactor the QR classes onto a templated base
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
54 QRP::init (const Matrix& a, qr_type_t qr_type) |
538 | 55 { |
9713
7918eb15040c
refactor the QR classes onto a templated base
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
56 assert (qr_type != qr_type_raw); |
538 | 57 |
5275 | 58 octave_idx_type m = a.rows (); |
59 octave_idx_type n = a.cols (); | |
538 | 60 |
5275 | 61 octave_idx_type min_mn = m < n ? m : n; |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
62 OCTAVE_LOCAL_BUFFER (double, tau, min_mn); |
1922 | 63 |
5275 | 64 octave_idx_type info = 0; |
538 | 65 |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
66 Matrix afact = a; |
9713
7918eb15040c
refactor the QR classes onto a templated base
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
67 if (m > n && qr_type == qr_type_std) |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
68 afact.resize (m, m); |
538 | 69 |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
70 MArray<octave_idx_type> jpvt (n, 1, 0); |
538 | 71 |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
72 if (m > 0) |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
73 { |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
74 // workspace query. |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
75 double rlwork; |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
76 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
|
77 tau, &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); |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
83 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
|
84 tau, work, lwork, info)); |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
85 } |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
86 else |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
87 for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; |
1922 | 88 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
89 // 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
|
90 // indices only!) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
91 |
8811 | 92 jpvt -= static_cast<octave_idx_type> (1); |
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
93 p = PermMatrix (jpvt, true); |
1922 | 94 |
95 | |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
96 form (n, afact, tau, qr_type); |
538 | 97 } |
98 | |
10905
9ee5a0a1b93d
Return permutation vector from QR as a row, not column, vector.
Rik <octave@nomad.inbox5.com>
parents:
10521
diff
changeset
|
99 RowVector |
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
100 QRP::Pvec (void) const |
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
101 { |
8375
e3c9102431a9
fix design problems of diag & perm matrix classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8368
diff
changeset
|
102 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
|
103 RowVector pv (MArray<double> (pa) + 1.0); |
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
104 return pv; |
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
105 } |