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