Mercurial > hg > octave-nkf

comparison liboctave/dbleQRP.cc @ 8597:c86718093c1b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improve & fix QR classes
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 27 Jan 2009 12:40:06 +0100
parents e3c9102431a9
children 20dfb885f877
comparison
equal deleted inserted replaced
8596:8833c0b18eb2 8597:c86718093c1b
28 #include <cassert> 28 #include <cassert>
29 29
30 #include "dbleQRP.h" 30 #include "dbleQRP.h"
31 #include "f77-fcn.h" 31 #include "f77-fcn.h"
32 #include "lo-error.h" 32 #include "lo-error.h"
33 #include "oct-locbuf.h"
33 34
34 extern "C" 35 extern "C"
35 { 36 {
36 F77_RET_T 37 F77_RET_T
37 F77_FUNC (dgeqp3, DGEQP3) (const octave_idx_type&, const octave_idx_type&, double*, 38 F77_FUNC (dgeqp3, DGEQP3) (const octave_idx_type&, const octave_idx_type&, double*,
38 const octave_idx_type&, octave_idx_type*, double*, double*, 39 const octave_idx_type&, octave_idx_type*, double*, double*,
39 const octave_idx_type&, octave_idx_type&); 40 const octave_idx_type&, octave_idx_type&);
40
41 F77_RET_T
42 F77_FUNC (dorgqr, DORGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
43 double*, const octave_idx_type&, double*, double*,
44 const octave_idx_type&, octave_idx_type&);
45 } 41 }
46 42
47 // It would be best to share some of this code with QR class... 43 // It would be best to share some of this code with QR class...
48 44
49 QRP::QRP (const Matrix& a, QR::type qr_type) 45 QRP::QRP (const Matrix& a, QR::type qr_type)
58 assert (qr_type != QR::raw); 54 assert (qr_type != QR::raw);
59 55
60 octave_idx_type m = a.rows (); 56 octave_idx_type m = a.rows ();
61 octave_idx_type n = a.cols (); 57 octave_idx_type n = a.cols ();
62 58
63 if (m == 0 || n == 0)
64 {
65 (*current_liboctave_error_handler) ("QR must have non-empty matrix");
66 return;
67 }
68
69 octave_idx_type min_mn = m < n ? m : n; 59 octave_idx_type min_mn = m < n ? m : n;
70 Array<double> tau (min_mn); 60 OCTAVE_LOCAL_BUFFER (double, tau, min_mn);
71 double *ptau = tau.fortran_vec ();
72 61
73 octave_idx_type info = 0; 62 octave_idx_type info = 0;
74 63
75 Matrix A_fact = a; 64 Matrix afact = a;
76 if (m > n && qr_type != QR::economy) 65 if (m > n && qr_type == QR::std)
77 A_fact.resize (m, m, 0.0); 66 afact.resize (m, m);
78
79 double *tmp_data = A_fact.fortran_vec ();
80 67
81 MArray<octave_idx_type> jpvt (n, 0); 68 MArray<octave_idx_type> jpvt (n, 0);
82 octave_idx_type *pjpvt = jpvt.fortran_vec ();
83 69
84 double rlwork = 0; 70 if (m > 0)
85 // Workspace query... 71 {
86 F77_XFCN (dgeqp3, DGEQP3, (m, n, tmp_data, m, pjpvt, ptau, &rlwork, -1, info)); 72 // workspace query.
73 double rlwork;
74 F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (),
75 tau, &rlwork, -1, info));
87 76
88 octave_idx_type lwork = rlwork; 77 // allocate buffer and do the job.
89 Array<double> work (lwork); 78 octave_idx_type lwork = rlwork; lwork = std::max (lwork, 1);
90 double *pwork = work.fortran_vec (); 79 OCTAVE_LOCAL_BUFFER (double, work, lwork);
91 80 F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (),
92 // Code to enforce a certain permutation could go here... 81 tau, work, lwork, info));
93 82 }
94 F77_XFCN (dgeqp3, DGEQP3, (m, n, tmp_data, m, pjpvt, ptau, pwork, lwork, info)); 83 else
84 for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1;
95 85
96 // Form Permutation matrix (if economy is requested, return the 86 // Form Permutation matrix (if economy is requested, return the
97 // indices only!) 87 // indices only!)
98 88
99 jpvt -= 1; 89 jpvt -= 1;
100 p = PermMatrix (jpvt, true); 90 p = PermMatrix (jpvt, true);
101 91
102 octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m;
103 92
104 if (qr_type == QR::economy && m > n) 93 form (n, afact, tau, qr_type);
105 r.resize (n, n, 0.0);
106 else
107 r.resize (m, n, 0.0);
108
109 for (octave_idx_type j = 0; j < n; j++)
110 {
111 octave_idx_type limit = j < min_mn-1 ? j : min_mn-1;
112 for (octave_idx_type i = 0; i <= limit; i++)
113 r.elem (i, j) = A_fact.elem (i, j);
114 }
115
116 F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau,
117 pwork, lwork, info));
118
119 q = A_fact;
120 q.resize (m, n2);
121 } 94 }
122 95
123 ColumnVector 96 ColumnVector
124 QRP::Pvec (void) const 97 QRP::Pvec (void) const
125 { 98 {