Bug Summary

File:liboctave/numeric/dbleQR.cc
Location:line 128, column 7
Description:Assigned value is garbage or undefined

Annotated Source Code

1/*
2
3Copyright (C) 1994-2013 John W. Eaton
4Copyright (C) 2008-2009 Jaroslav Hajek
5Copyright (C) 2009 VZLU Prague
6
7This file is part of Octave.
8
9Octave is free software; you can redistribute it and/or modify it
10under the terms of the GNU General Public License as published by the
11Free Software Foundation; either version 3 of the License, or (at your
12option) any later version.
13
14Octave is distributed in the hope that it will be useful, but WITHOUT
15ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17for more details.
18
19You should have received a copy of the GNU General Public License
20along with Octave; see the file COPYING. If not, see
21<http://www.gnu.org/licenses/>.
22
23*/
24
25#ifdef HAVE_CONFIG_H1
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
38template class base_qr<Matrix>;
39
40extern "C"
41{
42 F77_RET_Tint
43 F77_FUNC (dgeqrf, 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_Tint
49 F77_FUNC (dorgqr, 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_QRUPDATE1
55
56 F77_RET_Tint
57 F77_FUNC (dqr1up, 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_Tint
63 F77_FUNC (dqrinc, 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_Tint
70 F77_FUNC (dqrdec, 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_Tint
77 F77_FUNC (dqrinr, 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_Tint
83 F77_FUNC (dqrder, 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_Tint
89 F77_FUNC (dqrshc, 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
98const QR::type QR::raw, QR::std, QR::economy;
99
100QR::QR (const Matrix& a, qr_type_t qr_type)
101{
102 init (a, qr_type);
103}
104
105void
106QR::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;
1
'?' condition is false
112 OCTAVE_LOCAL_BUFFER (double, tau, min_mn)octave_local_buffer<double> _buffer_tau (min_mn); double
*tau = _buffer_tau
;
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)
2
Assuming 'm' is > 0
3
Taking true branch
121 {
122 // workspace query.
123 double rlwork;
4
'rlwork' declared without an initial value
124 F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dgeqrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dgeqrf_ (m, n, afact.fortran_vec (), m, tau, &rlwork,
-1, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
125 &rlwork, -1, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dgeqrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dgeqrf_ (m, n, afact.fortran_vec (), m, tau, &rlwork,
-1, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
126
127 // allocate buffer and do the job.
128 octave_idx_type lwork = rlwork;
5
Assigned value is garbage or undefined
129 lwork = std::max (lwork, static_cast<octave_idx_type> (1));
130 OCTAVE_LOCAL_BUFFER (double, work, lwork)octave_local_buffer<double> _buffer_work (lwork); double
*work = _buffer_work
;
131 F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dgeqrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dgeqrf_ (m, n, afact.fortran_vec (), m, tau, work, lwork,
info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
132 work, lwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dgeqrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dgeqrf_ (m, n, afact.fortran_vec (), m, tau, work, lwork,
info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
133 }
134
135 form (n, afact, tau, qr_type);
136}
137
138void QR::form (octave_idx_type n, Matrix& afact,
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,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dorgqr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dorgqr_ (m, k, min_mn, q.fortran_vec (), m, tau, &rlwork
, -1, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
194 &rlwork, -1, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dorgqr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dorgqr_ (m, k, min_mn, q.fortran_vec (), m, tau, &rlwork
, -1, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
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)octave_local_buffer<double> _buffer_work (lwork); double
*work = _buffer_work
;
200 F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dorgqr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dorgqr_ (m, k, min_mn, q.fortran_vec (), m, tau, work, lwork
, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
201 work, lwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dorgqr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dorgqr_ (m, k, min_mn, q.fortran_vec (), m, tau, work, lwork
, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
202 }
203 }
204}
205
206#ifdef HAVE_QRUPDATE1
207
208void
209QR::update (const ColumnVector& u, const ColumnVector& v)
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)octave_local_buffer<double> _buffer_w (2*k); double *w =
_buffer_w
;
219 F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqr1up_ (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k
, utmp.fortran_vec (), vtmp.fortran_vec (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
220 m, r.fortran_vec (), k,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqr1up_ (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k
, utmp.fortran_vec (), vtmp.fortran_vec (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
221 utmp.fortran_vec (), vtmp.fortran_vec (), w))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqr1up_ (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k
, utmp.fortran_vec (), vtmp.fortran_vec (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
222 }
223 else
224 (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
225}
226
227void
228QR::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)octave_local_buffer<double> _buffer_w (2*k); double *w =
_buffer_w
;
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 (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqr1up_ (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k
, utmp.fortran_vec (), vtmp.fortran_vec (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
241 m, r.fortran_vec (), k,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqr1up_ (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k
, utmp.fortran_vec (), vtmp.fortran_vec (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
242 utmp.fortran_vec (), vtmp.fortran_vec (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqr1up_ (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k
, utmp.fortran_vec (), vtmp.fortran_vec (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
243 w))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqr1up_ (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k
, utmp.fortran_vec (), vtmp.fortran_vec (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
244 }
245 }
246 else
247 (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
248}
249
250void
251QR::insert_col (const ColumnVector& u, octave_idx_type j)
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)octave_local_buffer<double> _buffer_w (k); double *w = _buffer_w;
275 F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), q.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrinc_ (m, n, k, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), j + 1, utmp.data (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
276 r.fortran_vec (), r.rows (), j + 1,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrinc_ (m, n, k, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), j + 1, utmp.data (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
277 utmp.data (), w))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrinc_ (m, n, k, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), j + 1, utmp.data (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
278 }
279}
280
281void
282QR::insert_col (const Matrix& u, const Array<octave_idx_type>& j)
283{
284 octave_idx_type m = q.rows ();
285 octave_idx_type n = r.columns ();
286 octave_idx_type k = q.columns ();
287
288 Array<octave_idx_type> jsi;
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)octave_local_buffer<double> _buffer_w (kmax); double *w
= _buffer_w
;
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),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrinc_ (m, n + ii, std::min (kmax, k + ii), q.fortran_vec
(), q.rows (), r.fortran_vec (), r.rows (), js(ii) + 1, utmp
.data (), w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
320 q.fortran_vec (), q.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrinc_ (m, n + ii, std::min (kmax, k + ii), q.fortran_vec
(), q.rows (), r.fortran_vec (), r.rows (), js(ii) + 1, utmp
.data (), w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
321 r.fortran_vec (), r.rows (), js(ii) + 1,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrinc_ (m, n + ii, std::min (kmax, k + ii), q.fortran_vec
(), q.rows (), r.fortran_vec (), r.rows (), js(ii) + 1, utmp
.data (), w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
322 utmp.data (), w))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrinc_ (m, n + ii, std::min (kmax, k + ii), q.fortran_vec
(), q.rows (), r.fortran_vec (), r.rows (), js(ii) + 1, utmp
.data (), w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
323 }
324 }
325}
326
327void
328QR::delete_col (octave_idx_type j)
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)octave_local_buffer<double> _buffer_w (k); double *w = _buffer_w;
339 F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), q.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrdec_ (m, n, k, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), j + 1, w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
340 r.fortran_vec (), r.rows (), j + 1, w))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrdec_ (m, n, k, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), j + 1, w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
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
354void
355QR::delete_col (const Array<octave_idx_type>& j)
356{
357 octave_idx_type m = q.rows ();
358 octave_idx_type n = r.columns ();
359 octave_idx_type k = q.columns ();
360
361 Array<octave_idx_type> jsi;
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)octave_local_buffer<double> _buffer_w (k); double *w = _buffer_w;
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,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrdec_ (m, n - ii, k == m ? k : k - ii, q.fortran_vec ()
, q.rows (), r.fortran_vec (), r.rows (), js(ii) + 1, w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
379 q.fortran_vec (), q.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrdec_ (m, n - ii, k == m ? k : k - ii, q.fortran_vec ()
, q.rows (), r.fortran_vec (), r.rows (), js(ii) + 1, w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
380 r.fortran_vec (), r.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrdec_ (m, n - ii, k == m ? k : k - ii, q.fortran_vec ()
, q.rows (), r.fortran_vec (), r.rows (), js(ii) + 1, w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
381 js(ii) + 1, w))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrdec_ (m, n - ii, k == m ? k : k - ii, q.fortran_vec ()
, q.rows (), r.fortran_vec (), r.rows (), js(ii) + 1, w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
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
396void
397QR::insert_row (const RowVector& u, octave_idx_type j)
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)octave_local_buffer<double> _buffer_w (k); double *w = _buffer_w;
413 F77_XFCN (dqrinr, DQRINR, (m, n, q.fortran_vec (), q.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrinr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrinr_ (m, n, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), j + 1, utmp.fortran_vec (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
414 r.fortran_vec (), r.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrinr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrinr_ (m, n, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), j + 1, utmp.fortran_vec (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
415 j + 1, utmp.fortran_vec (), w))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrinr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrinr_ (m, n, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), j + 1, utmp.fortran_vec (), w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
416
417 }
418}
419
420void
421QR::delete_row (octave_idx_type j)
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)octave_local_buffer<double> _buffer_w (2*m); double *w =
_buffer_w
;
433 F77_XFCN (dqrder, DQRDER, (m, n, q.fortran_vec (), q.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrder_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrder_ (m, n, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), j + 1, w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
434 r.fortran_vec (), r.rows (), j + 1,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrder_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrder_ (m, n, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), j + 1, w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
435 w))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrder_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrder_ (m, n, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), j + 1, w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
436
437 q.resize (m - 1, m - 1);
438 r.resize (m - 1, n);
439 }
440}
441
442void
443QR::shift_cols (octave_idx_type i, octave_idx_type j)
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)octave_local_buffer<double> _buffer_w (2*k); double *w =
_buffer_w
;
454 F77_XFCN (dqrshc, DQRSHC, (m, n, k,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrshc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrshc_ (m, n, k, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), i + 1, j + 1, w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
455 q.fortran_vec (), q.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrshc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrshc_ (m, n, k, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), i + 1, j + 1, w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
456 r.fortran_vec (), r.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrshc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrshc_ (m, n, k, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), i + 1, j + 1, w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
457 i + 1, j + 1, w))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dqrshc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dqrshc_ (m, n, k, q.fortran_vec (), q.rows (), r.fortran_vec
(), r.rows (), i + 1, j + 1, w); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
458 }
459}
460
461#else
462
463// Replacement update methods.
464
465void
466QR::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
481void
482QR::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
497static
498Matrix 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
510static
511Matrix 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
523static
524Matrix 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
531static
532Matrix 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
539static
540Matrix shift_cols (const Matrix& a,
541 octave_idx_type i, octave_idx_type j)
542{
543 octave_idx_type n = a.columns ();
544 Array<octave_idx_type> p (n);
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
560void
561QR::insert_col (const ColumnVector& u, octave_idx_type j)
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
578void
579QR::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
586 Array<octave_idx_type> jsi;
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
608void
609QR::delete_col (octave_idx_type j)
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
624void
625QR::delete_col (const Array<octave_idx_type>& j)
626{
627 warn_qrupdate_once ();
628
629 octave_idx_type m = q.rows ();
630 octave_idx_type n = r.columns ();
631
632 Array<octave_idx_type> jsi;
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
652void
653QR::insert_row (const RowVector& u, octave_idx_type j)
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
670void
671QR::delete_row (octave_idx_type j)
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
686void
687QR::shift_cols (octave_idx_type i, octave_idx_type j)
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
702void 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