Mercurial > octave-nkf

comparison liboctave/dbleQR.cc @ 8547:d66c9b6e506a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
imported patch qrupdate.diff
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 20 Jan 2009 21:16:42 +0100
parents 445d27d79f4e
children a6edd5c23cb5
comparison
equal deleted inserted replaced
8546:3d8a914c580e 8547:d66c9b6e506a
1 /* 1 /*
2 2
3 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 3 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
4 John W. Eaton 4 John W. Eaton
5 Copyright (C) 2008 Jaroslav Hajek 5 Copyright (C) 2008, 2009 Jaroslav Hajek
6 6
7 This file is part of Octave. 7 This file is part of Octave.
8 8
9 Octave is free software; you can redistribute it and/or modify it 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 10 under the terms of the GNU General Public License as published by the
29 #include "dbleQR.h" 29 #include "dbleQR.h"
30 #include "f77-fcn.h" 30 #include "f77-fcn.h"
31 #include "lo-error.h" 31 #include "lo-error.h"
32 #include "Range.h" 32 #include "Range.h"
33 #include "idx-vector.h" 33 #include "idx-vector.h"
34 #include "oct-locbuf.h"
34 35
35 extern "C" 36 extern "C"
36 { 37 {
37 F77_RET_T 38 F77_RET_T
38 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&,
40 41
41 F77_RET_T 42 F77_RET_T
42 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*,
43 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&);
44 45
45 // these come from qrupdate 46 #ifdef HAVE_QRUPDATE
46 47
47 F77_RET_T 48 F77_RET_T
48 F77_FUNC (dqr1up, DQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 49 F77_FUNC (dqr1up, DQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
49 double*, double*, const double*, const double*); 50 double*, const octave_idx_type&, double*, const octave_idx_type&,
51 double*, double*, double*);
50 52
51 F77_RET_T 53 F77_RET_T
52 F77_FUNC (dqrinc, DQRINC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 54 F77_FUNC (dqrinc, DQRINC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
53 double*, const double*, double*, const octave_idx_type&, const double*); 55 double*, const octave_idx_type&, double*, const octave_idx_type&,
56 const octave_idx_type&, const double*, double*);
54 57
55 F77_RET_T 58 F77_RET_T
56 F77_FUNC (dqrdec, DQRDEC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 59 F77_FUNC (dqrdec, DQRDEC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
57 double*, const double*, double*, const octave_idx_type&); 60 double*, const octave_idx_type&, double*, const octave_idx_type&,
61 const octave_idx_type&, double*);
58 62
59 F77_RET_T 63 F77_RET_T
60 F77_FUNC (dqrinr, DQRINR) (const octave_idx_type&, const octave_idx_type&, 64 F77_FUNC (dqrinr, DQRINR) (const octave_idx_type&, const octave_idx_type&,
61 const double*, double*, const double*, double*, 65 double*, const octave_idx_type&, double*, const octave_idx_type&,
62 const octave_idx_type&, const double*); 66 const octave_idx_type&, const double*, double*);
63 67
64 F77_RET_T 68 F77_RET_T
65 F77_FUNC (dqrder, DQRDER) (const octave_idx_type&, const octave_idx_type&, 69 F77_FUNC (dqrder, DQRDER) (const octave_idx_type&, const octave_idx_type&,
66 const double*, double*, const double*, double *, 70 double*, const octave_idx_type&, double*, const octave_idx_type&,
67 const octave_idx_type&); 71 const octave_idx_type&, double*);
68 72
69 F77_RET_T 73 F77_RET_T
70 F77_FUNC (dqrshc, DQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 74 F77_FUNC (dqrshc, DQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
71 double*, double*, const octave_idx_type&, const octave_idx_type&); 75 double*, const octave_idx_type&, double*, const octave_idx_type&,
76 const octave_idx_type&, const octave_idx_type&,
77 double*);
78
79 #endif
72 } 80 }
73 81
74 QR::QR (const Matrix& a, QR::type qr_type) 82 QR::QR (const Matrix& a, QR::type qr_type)
75 : q (), r () 83 : q (), r ()
76 { 84 {
159 167
160 this->q = q_arg; 168 this->q = q_arg;
161 this->r = r_arg; 169 this->r = r_arg;
162 } 170 }
163 171
172 #ifdef HAVE_QRUPDATE
173
174 void
175 QR::update (const ColumnVector& u, const ColumnVector& v)
176 {
177 octave_idx_type m = q.rows ();
178 octave_idx_type n = r.columns ();
179 octave_idx_type k = q.columns ();
180
181 if (u.length () == m && v.length () == n)
182 {
183 ColumnVector utmp = u, vtmp = v;
184 OCTAVE_LOCAL_BUFFER (double, w, 2*k);
185 F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k,
186 utmp.fortran_vec (), vtmp.fortran_vec (), w));
187 }
188 else
189 (*current_liboctave_error_handler) ("QR update dimensions mismatch");
190 }
191
164 void 192 void
165 QR::update (const Matrix& u, const Matrix& v) 193 QR::update (const Matrix& u, const Matrix& v)
166 { 194 {
167 octave_idx_type m = q.rows (); 195 octave_idx_type m = q.rows ();
168 octave_idx_type n = r.columns (); 196 octave_idx_type n = r.columns ();
169 octave_idx_type k = q.columns (); 197 octave_idx_type k = q.columns ();
170 198
171 if (u.length () == m && v.length () == n) 199 if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
172 F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), 200 {
173 u.data (), v.data ())); 201 OCTAVE_LOCAL_BUFFER (double, w, 2*k);
174 else 202 for (octave_idx_type i = 0; i < u.cols (); i++)
175 (*current_liboctave_error_handler) ("QR update dimensions mismatch"); 203 {
176 } 204 ColumnVector utmp = u.column (i), vtmp = v.column (i);
177 205 F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k,
178 void 206 utmp.fortran_vec (), vtmp.fortran_vec (), w));
179 QR::insert_col (const Matrix& u, octave_idx_type j) 207 }
208 }
209 else
210 (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
211 }
212
213 void
214 QR::insert_col (const ColumnVector& u, octave_idx_type j)
180 { 215 {
181 octave_idx_type m = q.rows (); 216 octave_idx_type m = q.rows ();
182 octave_idx_type n = r.columns (); 217 octave_idx_type n = r.columns ();
183 octave_idx_type k = q.columns (); 218 octave_idx_type k = q.columns ();
184 219
185 if (u.length () != m) 220 if (u.length () != m)
186 (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); 221 (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
187 else if (j < 0 || j > n) 222 else if (j < 0 || j > n)
188 (*current_liboctave_error_handler) ("QR insert index out of range"); 223 (*current_liboctave_error_handler) ("qrinsert: index out of range");
189 else 224 else
190 { 225 {
191 Matrix r1 (m, n+1); 226 if (k < m)
192 227 {
193 F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), r.data (), 228 q.resize (m, k+1);
194 r1.fortran_vec (), j+1, u.data ())); 229 r.resize (k+1, n+1);
195 230 }
196 r = r1; 231 else
232 {
233 r.resize (k, n+1);
234 }
235
236 ColumnVector utmp = u;
237 OCTAVE_LOCAL_BUFFER (double, w, k);
238 F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), q.rows (),
239 r.fortran_vec (), r.rows (), j + 1,
240 utmp.data (), w));
241 }
242 }
243
244 void
245 QR::insert_col (const Matrix& u, const Array<octave_idx_type>& j)
246 {
247 octave_idx_type m = q.rows ();
248 octave_idx_type n = r.columns ();
249 octave_idx_type k = q.columns ();
250
251 Array<octave_idx_type> jsi;
252 Array<octave_idx_type> js = j.sort (jsi, ASCENDING);
253 octave_idx_type nj = js.length ();
254 bool dups = false;
255 for (octave_idx_type i = 0; i < nj - 1; i++)
256 dups = dups && js(i) == js(i+1);
257
258 if (dups)
259 (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
260 else if (u.length () != m || u.columns () != nj)
261 (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
262 else if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
263 (*current_liboctave_error_handler) ("qrinsert: index out of range");
264 else if (nj > 0)
265 {
266 octave_idx_type kmax = std::min (k + nj, m);
267 if (k < m)
268 {
269 q.resize (m, kmax);
270 r.resize (kmax, n + nj);
271 }
272 else
273 {
274 r.resize (k, n + nj);
275 }
276
277 OCTAVE_LOCAL_BUFFER (double, w, kmax);
278 for (octave_idx_type i = 0; i < js.length (); i++)
279 {
280 ColumnVector utmp = u.column (jsi(i));
281 F77_XFCN (dqrinc, DQRINC, (m, n + i, std::min (kmax, k + i),
282 q.fortran_vec (), q.rows (),
283 r.fortran_vec (), r.rows (), js(i) + 1,
284 utmp.data (), w));
285 }
197 } 286 }
198 } 287 }
199 288
200 void 289 void
201 QR::delete_col (octave_idx_type j) 290 QR::delete_col (octave_idx_type j)
202 { 291 {
203 octave_idx_type m = q.rows (); 292 octave_idx_type m = q.rows ();
204 octave_idx_type k = r.rows (); 293 octave_idx_type k = r.rows ();
205 octave_idx_type n = r.columns (); 294 octave_idx_type n = r.columns ();
206 295
207 if (k < m && k < n) 296 if (j < 0 || j > n-1)
208 (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); 297 (*current_liboctave_error_handler) ("qrdelete: index out of range");
209 else if (j < 0 || j > n-1) 298 else
210 (*current_liboctave_error_handler) ("QR delete index out of range"); 299 {
211 else 300 OCTAVE_LOCAL_BUFFER (double, w, k);
212 { 301 F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
213 Matrix r1 (k, n-1); 302 r.fortran_vec (), r.rows (), j + 1, w));
214 303
215 F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), r.data (), 304 if (k < m)
216 r1.fortran_vec (), j+1)); 305 {
217 306 q.resize (m, k-1);
218 r = r1; 307 r.resize (k-1, n-1);
219 } 308 }
220 } 309 else
221 310 {
222 void 311 r.resize (k, n-1);
223 QR::insert_row (const Matrix& u, octave_idx_type j) 312 }
313 }
314 }
315
316 void
317 QR::delete_col (const Array<octave_idx_type>& j)
318 {
319 octave_idx_type m = q.rows ();
320 octave_idx_type n = r.columns ();
321 octave_idx_type k = q.columns ();
322
323 Array<octave_idx_type> jsi;
324 Array<octave_idx_type> js = j.sort (jsi, DESCENDING);
325 octave_idx_type nj = js.length ();
326 bool dups = false;
327 for (octave_idx_type i = 0; i < nj - 1; i++)
328 dups = dups && js(i) == js(i+1);
329
330 if (dups)
331 (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
332 else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
333 (*current_liboctave_error_handler) ("qrinsert: index out of range");
334 else if (nj > 0)
335 {
336 OCTAVE_LOCAL_BUFFER (double, w, k);
337 for (octave_idx_type i = 0; i < js.length (); i++)
338 {
339 F77_XFCN (dqrdec, DQRDEC, (m, n - i, k == m ? k : k - i,
340 q.fortran_vec (), q.rows (),
341 r.fortran_vec (), r.rows (), js(i) + 1, w));
342 }
343 if (k < m)
344 {
345 q.resize (m, k - nj);
346 r.resize (k - nj, n - nj);
347 }
348 else
349 {
350 r.resize (k, n - nj);
351 }
352
353 }
354 }
355
356 void
357 QR::insert_row (const RowVector& u, octave_idx_type j)
224 { 358 {
225 octave_idx_type m = r.rows (); 359 octave_idx_type m = r.rows ();
226 octave_idx_type n = r.columns (); 360 octave_idx_type n = r.columns ();
361 octave_idx_type k = std::min (m, n);
227 362
228 if (! q.is_square () || u.length () != n) 363 if (! q.is_square () || u.length () != n)
229 (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); 364 (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
230 else if (j < 0 || j > m) 365 else if (j < 0 || j > m)
231 (*current_liboctave_error_handler) ("QR insert index out of range"); 366 (*current_liboctave_error_handler) ("qrinsert: index out of range");
232 else 367 else
233 { 368 {
234 Matrix q1 (m+1, m+1); 369 q.resize (m + 1, m + 1);
235 Matrix r1 (m+1, n); 370 r.resize (m + 1, n);
236 371 RowVector utmp = u;
237 F77_XFCN (dqrinr, DQRINR, (m, n, q.data (), q1.fortran_vec (), 372 OCTAVE_LOCAL_BUFFER (double, w, k);
238 r.data (), r1.fortran_vec (), j+1, u.data ())); 373 F77_XFCN (dqrinr, DQRINR, (m, n, q.fortran_vec (), q.rows (),
239 374 r.fortran_vec (), r.rows (),
240 q = q1; 375 j + 1, utmp.fortran_vec (), w));
241 r = r1; 376
242 } 377 }
243 } 378 }
244 379
245 void 380 void
246 QR::delete_row (octave_idx_type j) 381 QR::delete_row (octave_idx_type j)
247 { 382 {
248 octave_idx_type m = r.rows (); 383 octave_idx_type m = r.rows ();
249 octave_idx_type n = r.columns (); 384 octave_idx_type n = r.columns ();
250 385
251 if (! q.is_square ()) 386 if (! q.is_square ())
252 (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); 387 (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
253 else if (j < 0 || j > m-1) 388 else if (j < 0 || j > m-1)
254 (*current_liboctave_error_handler) ("QR delete index out of range"); 389 (*current_liboctave_error_handler) ("qrdelete: index out of range");
255 else 390 else
256 { 391 {
257 Matrix q1 (m-1, m-1); 392 OCTAVE_LOCAL_BUFFER (double, w, 2*m);
258 Matrix r1 (m-1, n); 393 F77_XFCN (dqrder, DQRDER, (m, n, q.fortran_vec (), q.rows (),
259 394 r.fortran_vec (), r.rows (), j + 1,
260 F77_XFCN (dqrder, DQRDER, (m, n, q.data (), q1.fortran_vec (), 395 w));
261 r.data (), r1.fortran_vec (), j+1 )); 396
262 397 q.resize (m - 1, m - 1);
263 q = q1; 398 r.resize (m - 1, n);
264 r = r1;
265 } 399 }
266 } 400 }
267 401
268 void 402 void
269 QR::shift_cols (octave_idx_type i, octave_idx_type j) 403 QR::shift_cols (octave_idx_type i, octave_idx_type j)
271 octave_idx_type m = q.rows (); 405 octave_idx_type m = q.rows ();
272 octave_idx_type k = r.rows (); 406 octave_idx_type k = r.rows ();
273 octave_idx_type n = r.columns (); 407 octave_idx_type n = r.columns ();
274 408
275 if (i < 0 || i > n-1 || j < 0 || j > n-1) 409 if (i < 0 || i > n-1 || j < 0 || j > n-1)
276 (*current_liboctave_error_handler) ("QR shift index out of range"); 410 (*current_liboctave_error_handler) ("qrshift: index out of range");
277 else 411 else
278 F77_XFCN (dqrshc, DQRSHC, (m, n, k, q.fortran_vec (), r.fortran_vec (), i+1, j+1)); 412 {
279 } 413 OCTAVE_LOCAL_BUFFER (double, w, 2*k);
280 414 F77_XFCN (dqrshc, DQRSHC, (m, n, k,
281 void 415 q.fortran_vec (), q.rows (),
282 QR::economize (void) 416 r.fortran_vec (), r.rows (),
283 { 417 i + 1, j + 1, w));
284 octave_idx_type r_nc = r.columns (); 418 }
285 419 }
286 if (r.rows () > r_nc) 420
287 { 421 #endif
288 q.resize (q.rows (), r_nc);
289 r.resize (r_nc, r_nc);
290 }
291 }
292
293 422
294 /* 423 /*
295 ;;; Local Variables: *** 424 ;;; Local Variables: ***
296 ;;; mode: C++ *** 425 ;;; mode: C++ ***
297 ;;; End: *** 426 ;;; End: ***