Mercurial > hg > octave-lyh
annotate liboctave/CmplxQRP.cc @ 10521:4d1fc073fbb7
add some missing copyright stmts
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Wed, 14 Apr 2010 12:23:13 +0200 |
parents | 12884915a8e4 |
children | 9ee5a0a1b93d |
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 "CmplxQRP.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 (zgeqp3, ZGEQP3) (const octave_idx_type&, const octave_idx_type&, Complex*, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
40 const octave_idx_type&, octave_idx_type*, Complex*, Complex*, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
41 const octave_idx_type&, double*, octave_idx_type&); |
538 | 42 } |
43 | |
44 // It would be best to share some of this code with ComplexQR class... | |
45 | |
9713
7918eb15040c
refactor the QR classes onto a templated base
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
46 ComplexQRP::ComplexQRP (const ComplexMatrix& a, qr_type_t qr_type) |
2763 | 47 : ComplexQR (), 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 ComplexQRP::init (const ComplexMatrix& 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 (Complex, 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 ComplexMatrix 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 OCTAVE_LOCAL_BUFFER (double, rwork, 2*n); |
8368
c72c1c9bccdc
call blocked permuted qr factorization routines from LAPACK
Jaroslav Hajek <highegg@gmail.com>
parents:
8367
diff
changeset
|
74 |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
75 // workspace query. |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
76 Complex clwork; |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
77 F77_XFCN (zgeqp3, ZGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (), |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
78 tau, &clwork, -1, rwork, info)); |
8368
c72c1c9bccdc
call blocked permuted qr factorization routines from LAPACK
Jaroslav Hajek <highegg@gmail.com>
parents:
8367
diff
changeset
|
79 |
8597
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
80 // allocate buffer and do the job. |
8811 | 81 octave_idx_type lwork = clwork.real (); |
82 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
|
83 OCTAVE_LOCAL_BUFFER (Complex, work, lwork); |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
84 F77_XFCN (zgeqp3, ZGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (), |
c86718093c1b
improve & fix QR classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8375
diff
changeset
|
85 tau, work, lwork, rwork, info)); |
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 | |
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
100 ColumnVector |
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
101 ComplexQRP::Pvec (void) const |
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
102 { |
8375
e3c9102431a9
fix design problems of diag & perm matrix classes
Jaroslav Hajek <highegg@gmail.com>
parents:
8368
diff
changeset
|
103 Array<double> pa (p.pvec ()); |
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
104 ColumnVector pv (MArray<double> (pa) + 1.0); |
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 } |