GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
dbleQR.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2013 John W. Eaton
4 Copyright (C) 2008-2009 Jaroslav Hajek
5 Copyright (C) 2009 VZLU Prague
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
11 Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
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
20 along with Octave; see the file COPYING. If not, see
21 <http://www.gnu.org/licenses/>.
22 
23 */
24 
25 #ifdef HAVE_CONFIG_H
26 #include <config.h>
27 #endif
28 
29 #include "dbleQR.h"
30 #include "f77-fcn.h"
31 #include "lo-error.h"
32 #include "Range.h"
33 #include "idx-vector.h"
34 #include "oct-locbuf.h"
35 
36 #include "base-qr.cc"
37 
38 template class base_qr<Matrix>;
39 
40 extern "C"
41 {
42  F77_RET_T
43  F77_FUNC (dgeqrf, DGEQRF) (const octave_idx_type&, const octave_idx_type&,
44  double*, const octave_idx_type&, double*,
45  double*, const octave_idx_type&,
46  octave_idx_type&);
47 
48  F77_RET_T
49  F77_FUNC (dorgqr, DORGQR) (const octave_idx_type&, const octave_idx_type&,
50  const octave_idx_type&, double*,
51  const octave_idx_type&, double*, double*,
52  const octave_idx_type&, octave_idx_type&);
53 
54 #ifdef HAVE_QRUPDATE
55 
56  F77_RET_T
57  F77_FUNC (dqr1up, DQR1UP) (const octave_idx_type&, const octave_idx_type&,
58  const octave_idx_type&, double*,
59  const octave_idx_type&, double*,
60  const octave_idx_type&, double*, double*, double*);
61 
62  F77_RET_T
63  F77_FUNC (dqrinc, DQRINC) (const octave_idx_type&, const octave_idx_type&,
64  const octave_idx_type&, double*,
65  const octave_idx_type&, double*,
66  const octave_idx_type&, const octave_idx_type&,
67  const double*, double*);
68 
69  F77_RET_T
70  F77_FUNC (dqrdec, DQRDEC) (const octave_idx_type&, const octave_idx_type&,
71  const octave_idx_type&, double*,
72  const octave_idx_type&, double*,
73  const octave_idx_type&, const octave_idx_type&,
74  double*);
75 
76  F77_RET_T
77  F77_FUNC (dqrinr, DQRINR) (const octave_idx_type&, const octave_idx_type&,
78  double*, const octave_idx_type&, double*,
79  const octave_idx_type&, const octave_idx_type&,
80  const double*, double*);
81 
82  F77_RET_T
83  F77_FUNC (dqrder, DQRDER) (const octave_idx_type&, const octave_idx_type&,
84  double*, const octave_idx_type&, double*,
85  const octave_idx_type&, const octave_idx_type&,
86  double*);
87 
88  F77_RET_T
89  F77_FUNC (dqrshc, DQRSHC) (const octave_idx_type&, const octave_idx_type&,
90  const octave_idx_type&, double*,
91  const octave_idx_type&, double*,
92  const octave_idx_type&, const octave_idx_type&,
93  const octave_idx_type&, double*);
94 
95 #endif
96 }
97 
99 
100 QR::QR (const Matrix& a, qr_type_t qr_type)
101 {
102  init (a, qr_type);
103 }
104 
105 void
106 QR::init (const Matrix& a, qr_type_t qr_type)
107 {
108  octave_idx_type m = a.rows ();
109  octave_idx_type n = a.cols ();
110 
111  octave_idx_type min_mn = m < n ? m : n;
112  OCTAVE_LOCAL_BUFFER (double, tau, min_mn);
113 
114  octave_idx_type info = 0;
115 
116  Matrix afact = a;
117  if (m > n && qr_type == qr_type_std)
118  afact.resize (m, m);
119 
120  if (m > 0)
121  {
122  // workspace query.
123  double rlwork;
124  F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,
125  &rlwork, -1, info));
126 
127  // allocate buffer and do the job.
128  octave_idx_type lwork = rlwork;
129  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
130  OCTAVE_LOCAL_BUFFER (double, work, lwork);
131  F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,
132  work, lwork, info));
133  }
134 
135  form (n, afact, tau, qr_type);
136 }
137 
139  double *tau, qr_type_t qr_type)
140 {
141  octave_idx_type m = afact.rows (), min_mn = std::min (m, n);
142  octave_idx_type info;
143 
144  if (qr_type == qr_type_raw)
145  {
146  for (octave_idx_type j = 0; j < min_mn; j++)
147  {
148  octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
149  for (octave_idx_type i = limit + 1; i < m; i++)
150  afact.elem (i, j) *= tau[j];
151  }
152 
153  r = afact;
154  }
155  else
156  {
157  // Attempt to minimize copying.
158  if (m >= n)
159  {
160  // afact will become q.
161  q = afact;
162  octave_idx_type k = qr_type == qr_type_economy ? n : m;
163  r = Matrix (k, n);
164  for (octave_idx_type j = 0; j < n; j++)
165  {
166  octave_idx_type i = 0;
167  for (; i <= j; i++)
168  r.xelem (i, j) = afact.xelem (i, j);
169  for (; i < k; i++)
170  r.xelem (i, j) = 0;
171  }
172  afact = Matrix (); // optimize memory
173  }
174  else
175  {
176  // afact will become r.
177  q = Matrix (m, m);
178  for (octave_idx_type j = 0; j < m; j++)
179  for (octave_idx_type i = j + 1; i < m; i++)
180  {
181  q.xelem (i, j) = afact.xelem (i, j);
182  afact.xelem (i, j) = 0;
183  }
184  r = afact;
185  }
186 
187 
188  if (m > 0)
189  {
190  octave_idx_type k = q.columns ();
191  // workspace query.
192  double rlwork;
193  F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
194  &rlwork, -1, info));
195 
196  // allocate buffer and do the job.
197  octave_idx_type lwork = rlwork;
198  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
199  OCTAVE_LOCAL_BUFFER (double, work, lwork);
200  F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
201  work, lwork, info));
202  }
203  }
204 }
205 
206 #ifdef HAVE_QRUPDATE
207 
208 void
210 {
211  octave_idx_type m = q.rows ();
212  octave_idx_type n = r.columns ();
213  octave_idx_type k = q.columns ();
214 
215  if (u.length () == m && v.length () == n)
216  {
217  ColumnVector utmp = u, vtmp = v;
218  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
219  F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (),
220  m, r.fortran_vec (), k,
221  utmp.fortran_vec (), vtmp.fortran_vec (), w));
222  }
223  else
224  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
225 }
226 
227 void
228 QR::update (const Matrix& u, const Matrix& v)
229 {
230  octave_idx_type m = q.rows ();
231  octave_idx_type n = r.columns ();
232  octave_idx_type k = q.columns ();
233 
234  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
235  {
236  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
237  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
238  {
239  ColumnVector utmp = u.column (i), vtmp = v.column (i);
240  F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (),
241  m, r.fortran_vec (), k,
242  utmp.fortran_vec (), vtmp.fortran_vec (),
243  w));
244  }
245  }
246  else
247  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
248 }
249 
250 void
252 {
253  octave_idx_type m = q.rows ();
254  octave_idx_type n = r.columns ();
255  octave_idx_type k = q.columns ();
256 
257  if (u.length () != m)
258  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
259  else if (j < 0 || j > n)
260  (*current_liboctave_error_handler) ("qrinsert: index out of range");
261  else
262  {
263  if (k < m)
264  {
265  q.resize (m, k+1);
266  r.resize (k+1, n+1);
267  }
268  else
269  {
270  r.resize (k, n+1);
271  }
272 
273  ColumnVector utmp = u;
274  OCTAVE_LOCAL_BUFFER (double, w, k);
275  F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), q.rows (),
276  r.fortran_vec (), r.rows (), j + 1,
277  utmp.data (), w));
278  }
279 }
280 
281 void
283 {
284  octave_idx_type m = q.rows ();
285  octave_idx_type n = r.columns ();
286  octave_idx_type k = q.columns ();
287 
289  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
290  octave_idx_type nj = js.length ();
291  bool dups = false;
292  for (octave_idx_type i = 0; i < nj - 1; i++)
293  dups = dups && js(i) == js(i+1);
294 
295  if (dups)
296  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
297  else if (u.length () != m || u.columns () != nj)
298  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
299  else if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
300  (*current_liboctave_error_handler) ("qrinsert: index out of range");
301  else if (nj > 0)
302  {
303  octave_idx_type kmax = std::min (k + nj, m);
304  if (k < m)
305  {
306  q.resize (m, kmax);
307  r.resize (kmax, n + nj);
308  }
309  else
310  {
311  r.resize (k, n + nj);
312  }
313 
314  OCTAVE_LOCAL_BUFFER (double, w, kmax);
315  for (volatile octave_idx_type i = 0; i < js.length (); i++)
316  {
317  octave_idx_type ii = i;
318  ColumnVector utmp = u.column (jsi(i));
319  F77_XFCN (dqrinc, DQRINC, (m, n + ii, std::min (kmax, k + ii),
320  q.fortran_vec (), q.rows (),
321  r.fortran_vec (), r.rows (), js(ii) + 1,
322  utmp.data (), w));
323  }
324  }
325 }
326 
327 void
329 {
330  octave_idx_type m = q.rows ();
331  octave_idx_type k = r.rows ();
332  octave_idx_type n = r.columns ();
333 
334  if (j < 0 || j > n-1)
335  (*current_liboctave_error_handler) ("qrdelete: index out of range");
336  else
337  {
338  OCTAVE_LOCAL_BUFFER (double, w, k);
339  F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
340  r.fortran_vec (), r.rows (), j + 1, w));
341 
342  if (k < m)
343  {
344  q.resize (m, k-1);
345  r.resize (k-1, n-1);
346  }
347  else
348  {
349  r.resize (k, n-1);
350  }
351  }
352 }
353 
354 void
356 {
357  octave_idx_type m = q.rows ();
358  octave_idx_type n = r.columns ();
359  octave_idx_type k = q.columns ();
360 
362  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
363  octave_idx_type nj = js.length ();
364  bool dups = false;
365  for (octave_idx_type i = 0; i < nj - 1; i++)
366  dups = dups && js(i) == js(i+1);
367 
368  if (dups)
369  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
370  else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
371  (*current_liboctave_error_handler) ("qrinsert: index out of range");
372  else if (nj > 0)
373  {
374  OCTAVE_LOCAL_BUFFER (double, w, k);
375  for (volatile octave_idx_type i = 0; i < js.length (); i++)
376  {
377  octave_idx_type ii = i;
378  F77_XFCN (dqrdec, DQRDEC, (m, n - ii, k == m ? k : k - ii,
379  q.fortran_vec (), q.rows (),
380  r.fortran_vec (), r.rows (),
381  js(ii) + 1, w));
382  }
383  if (k < m)
384  {
385  q.resize (m, k - nj);
386  r.resize (k - nj, n - nj);
387  }
388  else
389  {
390  r.resize (k, n - nj);
391  }
392 
393  }
394 }
395 
396 void
398 {
399  octave_idx_type m = r.rows ();
400  octave_idx_type n = r.columns ();
401  octave_idx_type k = std::min (m, n);
402 
403  if (! q.is_square () || u.length () != n)
404  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
405  else if (j < 0 || j > m)
406  (*current_liboctave_error_handler) ("qrinsert: index out of range");
407  else
408  {
409  q.resize (m + 1, m + 1);
410  r.resize (m + 1, n);
411  RowVector utmp = u;
412  OCTAVE_LOCAL_BUFFER (double, w, k);
413  F77_XFCN (dqrinr, DQRINR, (m, n, q.fortran_vec (), q.rows (),
414  r.fortran_vec (), r.rows (),
415  j + 1, utmp.fortran_vec (), w));
416 
417  }
418 }
419 
420 void
422 {
423  octave_idx_type m = r.rows ();
424  octave_idx_type n = r.columns ();
425 
426  if (! q.is_square ())
427  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
428  else if (j < 0 || j > m-1)
429  (*current_liboctave_error_handler) ("qrdelete: index out of range");
430  else
431  {
432  OCTAVE_LOCAL_BUFFER (double, w, 2*m);
433  F77_XFCN (dqrder, DQRDER, (m, n, q.fortran_vec (), q.rows (),
434  r.fortran_vec (), r.rows (), j + 1,
435  w));
436 
437  q.resize (m - 1, m - 1);
438  r.resize (m - 1, n);
439  }
440 }
441 
442 void
444 {
445  octave_idx_type m = q.rows ();
446  octave_idx_type k = r.rows ();
447  octave_idx_type n = r.columns ();
448 
449  if (i < 0 || i > n-1 || j < 0 || j > n-1)
450  (*current_liboctave_error_handler) ("qrshift: index out of range");
451  else
452  {
453  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
454  F77_XFCN (dqrshc, DQRSHC, (m, n, k,
455  q.fortran_vec (), q.rows (),
456  r.fortran_vec (), r.rows (),
457  i + 1, j + 1, w));
458  }
459 }
460 
461 #else
462 
463 // Replacement update methods.
464 
465 void
466 QR::update (const ColumnVector& u, const ColumnVector& v)
467 {
468  warn_qrupdate_once ();
469 
470  octave_idx_type m = q.rows ();
471  octave_idx_type n = r.columns ();
472 
473  if (u.length () == m && v.length () == n)
474  {
475  init (q*r + Matrix (u) * Matrix (v).transpose (), get_type ());
476  }
477  else
478  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
479 }
480 
481 void
482 QR::update (const Matrix& u, const Matrix& v)
483 {
484  warn_qrupdate_once ();
485 
486  octave_idx_type m = q.rows ();
487  octave_idx_type n = r.columns ();
488 
489  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
490  {
491  init (q*r + u * v.transpose (), get_type ());
492  }
493  else
494  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
495 }
496 
497 static
498 Matrix insert_col (const Matrix& a, octave_idx_type i,
499  const ColumnVector& x)
500 {
501  Matrix retval (a.rows (), a.columns () + 1);
502  retval.assign (idx_vector::colon, idx_vector (0, i),
503  a.index (idx_vector::colon, idx_vector (0, i)));
504  retval.assign (idx_vector::colon, idx_vector (i), x);
505  retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
506  a.index (idx_vector::colon, idx_vector (i, a.columns ())));
507  return retval;
508 }
509 
510 static
511 Matrix insert_row (const Matrix& a, octave_idx_type i,
512  const RowVector& x)
513 {
514  Matrix retval (a.rows () + 1, a.columns ());
515  retval.assign (idx_vector (0, i), idx_vector::colon,
516  a.index (idx_vector (0, i), idx_vector::colon));
517  retval.assign (idx_vector (i), idx_vector::colon, x);
518  retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
519  a.index (idx_vector (i, a.rows ()), idx_vector::colon));
520  return retval;
521 }
522 
523 static
524 Matrix delete_col (const Matrix& a, octave_idx_type i)
525 {
526  Matrix retval = a;
527  retval.delete_elements (1, idx_vector (i));
528  return retval;
529 }
530 
531 static
532 Matrix delete_row (const Matrix& a, octave_idx_type i)
533 {
534  Matrix retval = a;
535  retval.delete_elements (0, idx_vector (i));
536  return retval;
537 }
538 
539 static
540 Matrix shift_cols (const Matrix& a,
542 {
543  octave_idx_type n = a.columns ();
545  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
546  if (i < j)
547  {
548  for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
549  p(j) = i;
550  }
551  else if (j < i)
552  {
553  p(j) = i;
554  for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
555  }
556 
557  return a.index (idx_vector::colon, idx_vector (p));
558 }
559 
560 void
562 {
563  warn_qrupdate_once ();
564 
565  octave_idx_type m = q.rows ();
566  octave_idx_type n = r.columns ();
567 
568  if (u.length () != m)
569  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
570  else if (j < 0 || j > n)
571  (*current_liboctave_error_handler) ("qrinsert: index out of range");
572  else
573  {
574  init (::insert_col (q*r, j, u), get_type ());
575  }
576 }
577 
578 void
579 QR::insert_col (const Matrix& u, const Array<octave_idx_type>& j)
580 {
581  warn_qrupdate_once ();
582 
583  octave_idx_type m = q.rows ();
584  octave_idx_type n = r.columns ();
585 
587  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
588  octave_idx_type nj = js.length ();
589  bool dups = false;
590  for (octave_idx_type i = 0; i < nj - 1; i++)
591  dups = dups && js(i) == js(i+1);
592 
593  if (dups)
594  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
595  else if (u.length () != m || u.columns () != nj)
596  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
597  else if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
598  (*current_liboctave_error_handler) ("qrinsert: index out of range");
599  else if (nj > 0)
600  {
601  Matrix a = q*r;
602  for (octave_idx_type i = 0; i < js.length (); i++)
603  a = ::insert_col (a, js(i), u.column (i));
604  init (a, get_type ());
605  }
606 }
607 
608 void
610 {
611  warn_qrupdate_once ();
612 
613  octave_idx_type m = q.rows ();
614  octave_idx_type n = r.columns ();
615 
616  if (j < 0 || j > n-1)
617  (*current_liboctave_error_handler) ("qrdelete: index out of range");
618  else
619  {
620  init (::delete_col (q*r, j), get_type ());
621  }
622 }
623 
624 void
626 {
627  warn_qrupdate_once ();
628 
629  octave_idx_type m = q.rows ();
630  octave_idx_type n = r.columns ();
631 
633  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
634  octave_idx_type nj = js.length ();
635  bool dups = false;
636  for (octave_idx_type i = 0; i < nj - 1; i++)
637  dups = dups && js(i) == js(i+1);
638 
639  if (dups)
640  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
641  else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
642  (*current_liboctave_error_handler) ("qrinsert: index out of range");
643  else if (nj > 0)
644  {
645  Matrix a = q*r;
646  for (octave_idx_type i = 0; i < js.length (); i++)
647  a = ::delete_col (a, js(i));
648  init (a, get_type ());
649  }
650 }
651 
652 void
654 {
655  warn_qrupdate_once ();
656 
657  octave_idx_type m = r.rows ();
658  octave_idx_type n = r.columns ();
659 
660  if (! q.is_square () || u.length () != n)
661  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
662  else if (j < 0 || j > m)
663  (*current_liboctave_error_handler) ("qrinsert: index out of range");
664  else
665  {
666  init (::insert_row (q*r, j, u), get_type ());
667  }
668 }
669 
670 void
672 {
673  octave_idx_type m = r.rows ();
674  octave_idx_type n = r.columns ();
675 
676  if (! q.is_square ())
677  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
678  else if (j < 0 || j > m-1)
679  (*current_liboctave_error_handler) ("qrdelete: index out of range");
680  else
681  {
682  init (::delete_row (q*r, j), get_type ());
683  }
684 }
685 
686 void
688 {
689  warn_qrupdate_once ();
690 
691  octave_idx_type m = q.rows ();
692  octave_idx_type n = r.columns ();
693 
694  if (i < 0 || i > n-1 || j < 0 || j > n-1)
695  (*current_liboctave_error_handler) ("qrshift: index out of range");
696  else
697  {
698  init (::shift_cols (q*r, i, j), get_type ());
699  }
700 }
701 
702 void warn_qrupdate_once (void)
703 {
704  static bool warned = false;
705  if (! warned)
706  {
707  (*current_liboctave_warning_handler)
708  ("In this version of Octave, QR & Cholesky updating routines\n"
709  "simply update the matrix and recalculate factorizations.\n"
710  "To use fast algorithms, link Octave with the qrupdate library.\n"
711  "See <http://sourceforge.net/projects/qrupdate>.\n");
712  warned = true;
713  }
714 }
715 
716 #endif