Mercurial > octave

comparison liboctave/dbleQR.cc @ 9713:7918eb15040c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
refactor the QR classes onto a templated base
author Jaroslav Hajek <highegg@gmail.com>
date Mon, 12 Oct 2009 10:42:49 +0200
parents b03953732530
children 1a9508872af0
comparison
equal deleted inserted replaced
9712:25e1e368618c 9713:7918eb15040c
32 #include "lo-error.h" 32 #include "lo-error.h"
33 #include "Range.h" 33 #include "Range.h"
34 #include "idx-vector.h" 34 #include "idx-vector.h"
35 #include "oct-locbuf.h" 35 #include "oct-locbuf.h"
36 36
37 #include "base-qr.cc"
38
39 template class base_qr<Matrix>;
40
37 extern "C" 41 extern "C"
38 { 42 {
39 F77_RET_T 43 F77_RET_T
40 F77_FUNC (dgeqrf, DGEQRF) (const octave_idx_type&, const octave_idx_type&, double*, const octave_idx_type&, 44 F77_FUNC (dgeqrf, DGEQRF) (const octave_idx_type&, const octave_idx_type&, double*, const octave_idx_type&,
41 double*, double*, const octave_idx_type&, octave_idx_type&); 45 double*, double*, const octave_idx_type&, octave_idx_type&);
78 double*); 82 double*);
79 83
80 #endif 84 #endif
81 } 85 }
82 86
83 QR::QR (const Matrix& a, QR::type qr_type) 87 QR::QR (const Matrix& a, qr_type_t qr_type)
84 : q (), r ()
85 { 88 {
86 init (a, qr_type); 89 init (a, qr_type);
87 } 90 }
88 91
89 void 92 void
90 QR::init (const Matrix& a, QR::type qr_type) 93 QR::init (const Matrix& a, qr_type_t qr_type)
91 { 94 {
92 octave_idx_type m = a.rows (); 95 octave_idx_type m = a.rows ();
93 octave_idx_type n = a.cols (); 96 octave_idx_type n = a.cols ();
94 97
95 octave_idx_type min_mn = m < n ? m : n; 98 octave_idx_type min_mn = m < n ? m : n;
96 OCTAVE_LOCAL_BUFFER (double, tau, min_mn); 99 OCTAVE_LOCAL_BUFFER (double, tau, min_mn);
97 100
98 octave_idx_type info = 0; 101 octave_idx_type info = 0;
99 102
100 Matrix afact = a; 103 Matrix afact = a;
101 if (m > n && qr_type == QR::std) 104 if (m > n && qr_type == qr_type_std)
102 afact.resize (m, m); 105 afact.resize (m, m);
103 106
104 if (m > 0) 107 if (m > 0)
105 { 108 {
106 // workspace query. 109 // workspace query.
116 119
117 form (n, afact, tau, qr_type); 120 form (n, afact, tau, qr_type);
118 } 121 }
119 122
120 void QR::form (octave_idx_type n, Matrix& afact, 123 void QR::form (octave_idx_type n, Matrix& afact,
121 double *tau, QR::type qr_type) 124 double *tau, qr_type_t qr_type)
122 { 125 {
123 octave_idx_type m = afact.rows (), min_mn = std::min (m, n); 126 octave_idx_type m = afact.rows (), min_mn = std::min (m, n);
124 octave_idx_type info; 127 octave_idx_type info;
125 128
126 if (qr_type == QR::raw) 129 if (qr_type == qr_type_raw)
127 { 130 {
128 for (octave_idx_type j = 0; j < min_mn; j++) 131 for (octave_idx_type j = 0; j < min_mn; j++)
129 { 132 {
130 octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; 133 octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
131 for (octave_idx_type i = limit + 1; i < m; i++) 134 for (octave_idx_type i = limit + 1; i < m; i++)
139 // Attempt to minimize copying. 142 // Attempt to minimize copying.
140 if (m >= n) 143 if (m >= n)
141 { 144 {
142 // afact will become q. 145 // afact will become q.
143 q = afact; 146 q = afact;
144 octave_idx_type k = qr_type == QR::economy ? n : m; 147 octave_idx_type k = qr_type == qr_type_economy ? n : m;
145 r = Matrix (k, n); 148 r = Matrix (k, n);
146 for (octave_idx_type j = 0; j < n; j++) 149 for (octave_idx_type j = 0; j < n; j++)
147 { 150 {
148 octave_idx_type i = 0; 151 octave_idx_type i = 0;
149 for (; i <= j; i++) 152 for (; i <= j; i++)
183 work, lwork, info)); 186 work, lwork, info));
184 } 187 }
185 } 188 }
186 } 189 }
187 190
188 QR::QR (const Matrix& q_arg, const Matrix& r_arg)
189 {
190 octave_idx_type qr = q_arg.rows (), qc = q_arg.columns ();
191 octave_idx_type rr = r_arg.rows (), rc = r_arg.columns ();
192 if (qc == rr && (qr == qc || (qr > qc && rr == rc)))
193 {
194 q = q_arg;
195 r = r_arg;
196 }
197 else
198 (*current_liboctave_error_handler) ("QR dimensions mismatch");
199 }
200
201 QR::type
202 QR::get_type (void) const
203 {
204 QR::type retval;
205 if (!q.is_empty () && q.is_square ())
206 retval = QR::std;
207 else if (q.rows () > q.columns () && r.is_square ())
208 retval = QR::economy;
209 else
210 retval = QR::raw;
211 return retval;
212 }
213
214 #ifdef HAVE_QRUPDATE 191 #ifdef HAVE_QRUPDATE
215 192
216 void 193 void
217 QR::update (const ColumnVector& u, const ColumnVector& v) 194 QR::update (const ColumnVector& u, const ColumnVector& v)
218 { 195 {