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
floatQR.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 "floatQR.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<FloatMatrix>;
39 
40 extern "C"
41 {
42  F77_RET_T
43  F77_FUNC (sgeqrf, SGEQRF) (const octave_idx_type&, const octave_idx_type&,
44  float*, const octave_idx_type&, float*, float*,
45  const octave_idx_type&, octave_idx_type&);
46 
47  F77_RET_T
48  F77_FUNC (sorgqr, SORGQR) (const octave_idx_type&, const octave_idx_type&,
49  const octave_idx_type&, float*,
50  const octave_idx_type&, float*, float*,
51  const octave_idx_type&, octave_idx_type&);
52 
53 #ifdef HAVE_QRUPDATE
54 
55  F77_RET_T
56  F77_FUNC (sqr1up, SQR1UP) (const octave_idx_type&, const octave_idx_type&,
57  const octave_idx_type&, float*,
58  const octave_idx_type&, float*,
59  const octave_idx_type&, float*, float*, float*);
60 
61  F77_RET_T
62  F77_FUNC (sqrinc, SQRINC) (const octave_idx_type&, const octave_idx_type&,
63  const octave_idx_type&, float*,
64  const octave_idx_type&, float*,
65  const octave_idx_type&,
66  const octave_idx_type&, const float*, float*);
67 
68  F77_RET_T
69  F77_FUNC (sqrdec, SQRDEC) (const octave_idx_type&, const octave_idx_type&,
70  const octave_idx_type&, float*,
71  const octave_idx_type&, float*,
72  const octave_idx_type&,
73  const octave_idx_type&, float*);
74 
75  F77_RET_T
76  F77_FUNC (sqrinr, SQRINR) (const octave_idx_type&, const octave_idx_type&,
77  float*, const octave_idx_type&,
78  float*, const octave_idx_type&,
79  const octave_idx_type&, const float*, float*);
80 
81  F77_RET_T
82  F77_FUNC (sqrder, SQRDER) (const octave_idx_type&, const octave_idx_type&,
83  float*, const octave_idx_type&,
84  float*, const octave_idx_type&,
85  const octave_idx_type&, float*);
86 
87  F77_RET_T
88  F77_FUNC (sqrshc, SQRSHC) (const octave_idx_type&, const octave_idx_type&,
89  const octave_idx_type&, float*,
90  const octave_idx_type&, float*,
91  const octave_idx_type&,
92  const octave_idx_type&, const octave_idx_type&,
93  float*);
94 
95 #endif
96 }
97 
99 {
100  init (a, qr_type);
101 }
102 
103 void
105 {
106  octave_idx_type m = a.rows ();
107  octave_idx_type n = a.cols ();
108 
109  octave_idx_type min_mn = m < n ? m : n;
110  OCTAVE_LOCAL_BUFFER (float, tau, min_mn);
111 
112  octave_idx_type info = 0;
113 
114  FloatMatrix afact = a;
115  if (m > n && qr_type == qr_type_std)
116  afact.resize (m, m);
117 
118  if (m > 0)
119  {
120  // workspace query.
121  float rlwork;
122  F77_XFCN (sgeqrf, SGEQRF, (m, n, afact.fortran_vec (), m, tau,
123  &rlwork, -1, info));
124 
125  // allocate buffer and do the job.
126  octave_idx_type lwork = rlwork;
127  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
128  OCTAVE_LOCAL_BUFFER (float, work, lwork);
129  F77_XFCN (sgeqrf, SGEQRF, (m, n, afact.fortran_vec (), m, tau,
130  work, lwork, info));
131  }
132 
133  form (n, afact, tau, qr_type);
134 }
135 
137  float *tau, qr_type_t qr_type)
138 {
139  octave_idx_type m = afact.rows (), min_mn = std::min (m, n);
140  octave_idx_type info;
141 
142  if (qr_type == qr_type_raw)
143  {
144  for (octave_idx_type j = 0; j < min_mn; j++)
145  {
146  octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
147  for (octave_idx_type i = limit + 1; i < m; i++)
148  afact.elem (i, j) *= tau[j];
149  }
150 
151  r = afact;
152  }
153  else
154  {
155  // Attempt to minimize copying.
156  if (m >= n)
157  {
158  // afact will become q.
159  q = afact;
160  octave_idx_type k = qr_type == qr_type_economy ? n : m;
161  r = FloatMatrix (k, n);
162  for (octave_idx_type j = 0; j < n; j++)
163  {
164  octave_idx_type i = 0;
165  for (; i <= j; i++)
166  r.xelem (i, j) = afact.xelem (i, j);
167  for (; i < k; i++)
168  r.xelem (i, j) = 0;
169  }
170  afact = FloatMatrix (); // optimize memory
171  }
172  else
173  {
174  // afact will become r.
175  q = FloatMatrix (m, m);
176  for (octave_idx_type j = 0; j < m; j++)
177  for (octave_idx_type i = j + 1; i < m; i++)
178  {
179  q.xelem (i, j) = afact.xelem (i, j);
180  afact.xelem (i, j) = 0;
181  }
182  r = afact;
183  }
184 
185 
186  if (m > 0)
187  {
188  octave_idx_type k = q.columns ();
189  // workspace query.
190  float rlwork;
191  F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
192  &rlwork, -1, info));
193 
194  // allocate buffer and do the job.
195  octave_idx_type lwork = rlwork;
196  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
197  OCTAVE_LOCAL_BUFFER (float, work, lwork);
198  F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
199  work, lwork, info));
200  }
201  }
202 }
203 
204 #ifdef HAVE_QRUPDATE
205 
206 void
208 {
209  octave_idx_type m = q.rows ();
210  octave_idx_type n = r.columns ();
211  octave_idx_type k = q.columns ();
212 
213  if (u.length () == m && v.length () == n)
214  {
215  FloatColumnVector utmp = u, vtmp = v;
216  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
217  F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (),
218  m, r.fortran_vec (), k,
219  utmp.fortran_vec (), vtmp.fortran_vec (), w));
220  }
221  else
222  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
223 }
224 
225 void
227 {
228  octave_idx_type m = q.rows ();
229  octave_idx_type n = r.columns ();
230  octave_idx_type k = q.columns ();
231 
232  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
233  {
234  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
235  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
236  {
237  FloatColumnVector utmp = u.column (i), vtmp = v.column (i);
238  F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (),
239  m, r.fortran_vec (), k,
240  utmp.fortran_vec (), vtmp.fortran_vec (),
241  w));
242  }
243  }
244  else
245  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
246 }
247 
248 void
250 {
251  octave_idx_type m = q.rows ();
252  octave_idx_type n = r.columns ();
253  octave_idx_type k = q.columns ();
254 
255  if (u.length () != m)
256  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
257  else if (j < 0 || j > n)
258  (*current_liboctave_error_handler) ("qrinsert: index out of range");
259  else
260  {
261  if (k < m)
262  {
263  q.resize (m, k+1);
264  r.resize (k+1, n+1);
265  }
266  else
267  {
268  r.resize (k, n+1);
269  }
270 
271  FloatColumnVector utmp = u;
272  OCTAVE_LOCAL_BUFFER (float, w, k);
273  F77_XFCN (sqrinc, SQRINC, (m, n, k, q.fortran_vec (), q.rows (),
274  r.fortran_vec (), r.rows (), j + 1,
275  utmp.data (), w));
276  }
277 }
278 
279 void
281 {
282  octave_idx_type m = q.rows ();
283  octave_idx_type n = r.columns ();
284  octave_idx_type k = q.columns ();
285 
287  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
288  octave_idx_type nj = js.length ();
289  bool dups = false;
290  for (octave_idx_type i = 0; i < nj - 1; i++)
291  dups = dups && js(i) == js(i+1);
292 
293  if (dups)
294  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
295  else if (u.length () != m || u.columns () != nj)
296  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
297  else if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
298  (*current_liboctave_error_handler) ("qrinsert: index out of range");
299  else if (nj > 0)
300  {
301  octave_idx_type kmax = std::min (k + nj, m);
302  if (k < m)
303  {
304  q.resize (m, kmax);
305  r.resize (kmax, n + nj);
306  }
307  else
308  {
309  r.resize (k, n + nj);
310  }
311 
312  OCTAVE_LOCAL_BUFFER (float, w, kmax);
313  for (volatile octave_idx_type i = 0; i < js.length (); i++)
314  {
315  octave_idx_type ii = i;
316  FloatColumnVector utmp = u.column (jsi(i));
317  F77_XFCN (sqrinc, SQRINC, (m, n + ii, std::min (kmax, k + ii),
318  q.fortran_vec (), q.rows (),
319  r.fortran_vec (), r.rows (), js(ii) + 1,
320  utmp.data (), w));
321  }
322  }
323 }
324 
325 void
327 {
328  octave_idx_type m = q.rows ();
329  octave_idx_type k = r.rows ();
330  octave_idx_type n = r.columns ();
331 
332  if (j < 0 || j > n-1)
333  (*current_liboctave_error_handler) ("qrdelete: index out of range");
334  else
335  {
336  OCTAVE_LOCAL_BUFFER (float, w, k);
337  F77_XFCN (sqrdec, SQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
338  r.fortran_vec (), r.rows (), j + 1, w));
339 
340  if (k < m)
341  {
342  q.resize (m, k-1);
343  r.resize (k-1, n-1);
344  }
345  else
346  {
347  r.resize (k, n-1);
348  }
349  }
350 }
351 
352 void
354 {
355  octave_idx_type m = q.rows ();
356  octave_idx_type n = r.columns ();
357  octave_idx_type k = q.columns ();
358 
360  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
361  octave_idx_type nj = js.length ();
362  bool dups = false;
363  for (octave_idx_type i = 0; i < nj - 1; i++)
364  dups = dups && js(i) == js(i+1);
365 
366  if (dups)
367  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
368  else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
369  (*current_liboctave_error_handler) ("qrinsert: index out of range");
370  else if (nj > 0)
371  {
372  OCTAVE_LOCAL_BUFFER (float, w, k);
373  for (volatile octave_idx_type i = 0; i < js.length (); i++)
374  {
375  octave_idx_type ii = i;
376  F77_XFCN (sqrdec, SQRDEC, (m, n - ii, k == m ? k : k - ii,
377  q.fortran_vec (), q.rows (),
378  r.fortran_vec (), r.rows (),
379  js(ii) + 1, w));
380  }
381  if (k < m)
382  {
383  q.resize (m, k - nj);
384  r.resize (k - nj, n - nj);
385  }
386  else
387  {
388  r.resize (k, n - nj);
389  }
390 
391  }
392 }
393 
394 void
396 {
397  octave_idx_type m = r.rows ();
398  octave_idx_type n = r.columns ();
399  octave_idx_type k = std::min (m, n);
400 
401  if (! q.is_square () || u.length () != n)
402  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
403  else if (j < 0 || j > m)
404  (*current_liboctave_error_handler) ("qrinsert: index out of range");
405  else
406  {
407  q.resize (m + 1, m + 1);
408  r.resize (m + 1, n);
409  FloatRowVector utmp = u;
410  OCTAVE_LOCAL_BUFFER (float, w, k);
411  F77_XFCN (sqrinr, SQRINR, (m, n, q.fortran_vec (), q.rows (),
412  r.fortran_vec (), r.rows (),
413  j + 1, utmp.fortran_vec (), w));
414 
415  }
416 }
417 
418 void
420 {
421  octave_idx_type m = r.rows ();
422  octave_idx_type n = r.columns ();
423 
424  if (! q.is_square ())
425  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
426  else if (j < 0 || j > m-1)
427  (*current_liboctave_error_handler) ("qrdelete: index out of range");
428  else
429  {
430  OCTAVE_LOCAL_BUFFER (float, w, 2*m);
431  F77_XFCN (sqrder, SQRDER, (m, n, q.fortran_vec (), q.rows (),
432  r.fortran_vec (), r.rows (), j + 1,
433  w));
434 
435  q.resize (m - 1, m - 1);
436  r.resize (m - 1, n);
437  }
438 }
439 
440 void
442 {
443  octave_idx_type m = q.rows ();
444  octave_idx_type k = r.rows ();
445  octave_idx_type n = r.columns ();
446 
447  if (i < 0 || i > n-1 || j < 0 || j > n-1)
448  (*current_liboctave_error_handler) ("qrshift: index out of range");
449  else
450  {
451  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
452  F77_XFCN (sqrshc, SQRSHC, (m, n, k,
453  q.fortran_vec (), q.rows (),
454  r.fortran_vec (), r.rows (),
455  i + 1, j + 1, w));
456  }
457 }
458 
459 #else
460 
461 // Replacement update methods.
462 
463 void
465 {
466  warn_qrupdate_once ();
467 
468  octave_idx_type m = q.rows ();
469  octave_idx_type n = r.columns ();
470 
471  if (u.length () == m && v.length () == n)
472  {
473  init (q*r + FloatMatrix (u) * FloatMatrix (v).transpose (), get_type ());
474  }
475  else
476  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
477 }
478 
479 void
480 FloatQR::update (const FloatMatrix& u, const FloatMatrix& v)
481 {
482  warn_qrupdate_once ();
483 
484  octave_idx_type m = q.rows ();
485  octave_idx_type n = r.columns ();
486 
487  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
488  {
489  init (q*r + u * v.transpose (), get_type ());
490  }
491  else
492  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
493 }
494 
495 static
496 FloatMatrix insert_col (const FloatMatrix& a, octave_idx_type i,
497  const FloatColumnVector& x)
498 {
499  FloatMatrix retval (a.rows (), a.columns () + 1);
500  retval.assign (idx_vector::colon, idx_vector (0, i),
501  a.index (idx_vector::colon, idx_vector (0, i)));
502  retval.assign (idx_vector::colon, idx_vector (i), x);
503  retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
504  a.index (idx_vector::colon, idx_vector (i, a.columns ())));
505  return retval;
506 }
507 
508 static
509 FloatMatrix insert_row (const FloatMatrix& a, octave_idx_type i,
510  const FloatRowVector& x)
511 {
512  FloatMatrix retval (a.rows () + 1, a.columns ());
513  retval.assign (idx_vector (0, i), idx_vector::colon,
514  a.index (idx_vector (0, i), idx_vector::colon));
515  retval.assign (idx_vector (i), idx_vector::colon, x);
516  retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
517  a.index (idx_vector (i, a.rows ()), idx_vector::colon));
518  return retval;
519 }
520 
521 static
522 FloatMatrix delete_col (const FloatMatrix& a, octave_idx_type i)
523 {
524  FloatMatrix retval = a;
525  retval.delete_elements (1, idx_vector (i));
526  return retval;
527 }
528 
529 static
530 FloatMatrix delete_row (const FloatMatrix& a, octave_idx_type i)
531 {
532  FloatMatrix retval = a;
533  retval.delete_elements (0, idx_vector (i));
534  return retval;
535 }
536 
537 static
538 FloatMatrix shift_cols (const FloatMatrix& a,
540 {
541  octave_idx_type n = a.columns ();
543  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
544  if (i < j)
545  {
546  for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
547  p(j) = i;
548  }
549  else if (j < i)
550  {
551  p(j) = i;
552  for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
553  }
554 
555  return a.index (idx_vector::colon, idx_vector (p));
556 }
557 
558 void
560 {
561  warn_qrupdate_once ();
562 
563  octave_idx_type m = q.rows ();
564  octave_idx_type n = r.columns ();
565 
566  if (u.length () != m)
567  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
568  else if (j < 0 || j > n)
569  (*current_liboctave_error_handler) ("qrinsert: index out of range");
570  else
571  {
572  init (::insert_col (q*r, j, u), get_type ());
573  }
574 }
575 
576 void
578 {
579  warn_qrupdate_once ();
580 
581  octave_idx_type m = q.rows ();
582  octave_idx_type n = r.columns ();
583 
585  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
586  octave_idx_type nj = js.length ();
587  bool dups = false;
588  for (octave_idx_type i = 0; i < nj - 1; i++)
589  dups = dups && js(i) == js(i+1);
590 
591  if (dups)
592  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
593  else if (u.length () != m || u.columns () != nj)
594  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
595  else if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
596  (*current_liboctave_error_handler) ("qrinsert: index out of range");
597  else if (nj > 0)
598  {
599  FloatMatrix a = q*r;
600  for (octave_idx_type i = 0; i < js.length (); i++)
601  a = ::insert_col (a, js(i), u.column (i));
602  init (a, get_type ());
603  }
604 }
605 
606 void
608 {
609  warn_qrupdate_once ();
610 
611  octave_idx_type m = q.rows ();
612  octave_idx_type n = r.columns ();
613 
614  if (j < 0 || j > n-1)
615  (*current_liboctave_error_handler) ("qrdelete: index out of range");
616  else
617  {
618  init (::delete_col (q*r, j), get_type ());
619  }
620 }
621 
622 void
624 {
625  warn_qrupdate_once ();
626 
627  octave_idx_type m = q.rows ();
628  octave_idx_type n = r.columns ();
629 
631  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
632  octave_idx_type nj = js.length ();
633  bool dups = false;
634  for (octave_idx_type i = 0; i < nj - 1; i++)
635  dups = dups && js(i) == js(i+1);
636 
637  if (dups)
638  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
639  else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
640  (*current_liboctave_error_handler) ("qrinsert: index out of range");
641  else if (nj > 0)
642  {
643  FloatMatrix a = q*r;
644  for (octave_idx_type i = 0; i < js.length (); i++)
645  a = ::delete_col (a, js(i));
646  init (a, get_type ());
647  }
648 }
649 
650 void
652 {
653  warn_qrupdate_once ();
654 
655  octave_idx_type m = r.rows ();
656  octave_idx_type n = r.columns ();
657 
658  if (! q.is_square () || u.length () != n)
659  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
660  else if (j < 0 || j > m)
661  (*current_liboctave_error_handler) ("qrinsert: index out of range");
662  else
663  {
664  init (::insert_row (q*r, j, u), get_type ());
665  }
666 }
667 
668 void
670 {
671  warn_qrupdate_once ();
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 #endif