Mercurial > hg > octave-jordi

comparison liboctave/dbleQR.cc @ 7553:56be6f31dd4e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implementation of QR factorization updating
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 04 Mar 2008 21:47:11 -0500
parents 29980c6b8604
children 40574114c514
comparison
equal deleted inserted replaced
7552:6070c3bd69c4 7553:56be6f31dd4e
19 along with Octave; see the file COPYING. If not, see 19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>. 20 <http://www.gnu.org/licenses/>.
21 21
22 */ 22 */
23 23
24 // updating/downdating by Jaroslav Hajek 2008
25
24 #ifdef HAVE_CONFIG_H 26 #ifdef HAVE_CONFIG_H
25 #include <config.h> 27 #include <config.h>
26 #endif 28 #endif
27 29
28 #include "dbleQR.h" 30 #include "dbleQR.h"
29 #include "f77-fcn.h" 31 #include "f77-fcn.h"
30 #include "lo-error.h" 32 #include "lo-error.h"
33 #include "Range.h"
34 #include "idx-vector.h"
31 35
32 extern "C" 36 extern "C"
33 { 37 {
34 F77_RET_T 38 F77_RET_T
35 F77_FUNC (dgeqrf, DGEQRF) (const octave_idx_type&, const octave_idx_type&, double*, const octave_idx_type&, 39 F77_FUNC (dgeqrf, DGEQRF) (const octave_idx_type&, const octave_idx_type&, double*, const octave_idx_type&,
36 double*, double*, const octave_idx_type&, octave_idx_type&); 40 double*, double*, const octave_idx_type&, octave_idx_type&);
37 41
38 F77_RET_T 42 F77_RET_T
39 F77_FUNC (dorgqr, DORGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, double*, 43 F77_FUNC (dorgqr, DORGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, double*,
40 const octave_idx_type&, double*, double*, const octave_idx_type&, octave_idx_type&); 44 const octave_idx_type&, double*, double*, const octave_idx_type&, octave_idx_type&);
45
46 // these come from qrupdate
47
48 F77_RET_T
49 F77_FUNC (dqr1up, DQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
50 double*, double*, const double*, const double*);
51
52 F77_RET_T
53 F77_FUNC (dqrinc, DQRINC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
54 double*, const double*, double*, const octave_idx_type&, const double*);
55
56 F77_RET_T
57 F77_FUNC (dqrdec, DQRDEC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
58 double*, const double*, double*, const octave_idx_type&);
59
60 F77_RET_T
61 F77_FUNC (dqrinr, DQRINR) (const octave_idx_type&, const octave_idx_type&,
62 const double*, double*, const double*, double*,
63 const octave_idx_type&, const double*);
64
65 F77_RET_T
66 F77_FUNC (dqrder, DQRDER) (const octave_idx_type&, const octave_idx_type&,
67 const double*, double*, const double*, double *,
68 const octave_idx_type&);
41 } 69 }
42 70
43 QR::QR (const Matrix& a, QR::type qr_type) 71 QR::QR (const Matrix& a, QR::type qr_type)
44 : q (), r () 72 : q (), r ()
45 { 73 {
116 q = A_fact; 144 q = A_fact;
117 q.resize (m, n2); 145 q.resize (m, n2);
118 } 146 }
119 } 147 }
120 148
149 QR::QR (const Matrix& q, const Matrix& r)
150 {
151 if (q.columns () != r.rows ())
152 {
153 (*current_liboctave_error_handler) ("QR dimensions mismatch");
154 return;
155 }
156
157 this->q = q;
158 this->r = r;
159 }
160
161 void
162 QR::update (const Matrix& u, const Matrix& v)
163 {
164 octave_idx_type m = q.rows ();
165 octave_idx_type n = r.columns ();
166 octave_idx_type k = q.columns ();
167
168 if (u.length () == m && v.length () == n)
169 F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (),
170 u.data (), v.data ()));
171 else
172 (*current_liboctave_error_handler) ("QR update dimensions mismatch");
173 }
174
175 void
176 QR::insert_col (const Matrix& u, octave_idx_type j)
177 {
178 octave_idx_type m = q.rows ();
179 octave_idx_type n = r.columns ();
180 octave_idx_type k = q.columns ();
181
182 if (u.length () != m)
183 (*current_liboctave_error_handler) ("QR insert dimensions mismatch");
184 else if (j < 1 || j > n+1)
185 (*current_liboctave_error_handler) ("QR insert index out of range");
186 else
187 {
188 Matrix r1 (m, n+1);
189
190 F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), r.data (),
191 r1.fortran_vec (), j, u.data ()));
192
193 r = r1;
194 }
195 }
196
197 void
198 QR::delete_col (octave_idx_type j)
199 {
200 octave_idx_type m = q.rows ();
201 octave_idx_type k = r.rows ();
202 octave_idx_type n = r.columns ();
203
204 if (k < m && k < n)
205 (*current_liboctave_error_handler) ("QR delete dimensions mismatch");
206 else if (j < 1 || j > n)
207 (*current_liboctave_error_handler) ("QR delete index out of range");
208 else
209 {
210 Matrix r1 (k, n-1);
211
212 F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), r.data (),
213 r1.fortran_vec (), j));
214
215 r = r1;
216 }
217 }
218
219 void
220 QR::insert_row (const Matrix& u, octave_idx_type j)
221 {
222 octave_idx_type m = r.rows ();
223 octave_idx_type n = r.columns ();
224
225 if (! q.is_square () || u.length () != n)
226 (*current_liboctave_error_handler) ("QR insert dimensions mismatch");
227 else if (j < 1 || j > m+1)
228 (*current_liboctave_error_handler) ("QR insert index out of range");
229 else
230 {
231 Matrix q1 (m+1, m+1);
232 Matrix r1 (m+1, n);
233
234 F77_XFCN (dqrinr, DQRINR, (m, n, q.data (), q1.fortran_vec (),
235 r.data (), r1.fortran_vec (), j, u.data ()));
236
237 q = q1;
238 r = r1;
239 }
240 }
241
242 void
243 QR::delete_row (octave_idx_type j)
244 {
245 octave_idx_type m = r.rows ();
246 octave_idx_type n = r.columns ();
247
248 if (! q.is_square ())
249 (*current_liboctave_error_handler) ("QR insert dimensions mismatch");
250 else if (j < 1 || j > n)
251 (*current_liboctave_error_handler) ("QR delete index out of range");
252 else
253 {
254 Matrix q1 (m-1, m-1);
255 Matrix r1 (m-1, n);
256
257 F77_XFCN (dqrder, DQRDER, (m, n, q.data (), q1.fortran_vec (),
258 r.data (), r1.fortran_vec (), j ));
259
260 q = q1;
261 r = r1;
262 }
263 }
264
265 void
266 QR::economize (void)
267 {
268 idx_vector c (':'), i (Range (1, r.columns ()));
269 q = Matrix (q.index (c, i));
270 r = Matrix (r.index (i, c));
271 #if 0
272 octave_idx_type r_nc = r.columns ();
273
274 if (r.rows () > r_nc)
275 {
276 q.resize (q.rows (), r_nc);
277 r.resize (r_nc, r_nc);
278 }
279 #endif
280 }
281
282
121 /* 283 /*
122 ;;; Local Variables: *** 284 ;;; Local Variables: ***
123 ;;; mode: C++ *** 285 ;;; mode: C++ ***
124 ;;; End: *** 286 ;;; End: ***
125 */ 287 */