Bug Summary

File:liboctave/numeric/floatQR.cc
Location:line 195, column 11
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 "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
38template class base_qr<FloatMatrix>;
39
40extern "C"
41{
42 F77_RET_Tint
43 F77_FUNC (sgeqrf, 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_Tint
48 F77_FUNC (sorgqr, 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_QRUPDATE1
54
55 F77_RET_Tint
56 F77_FUNC (sqr1up, 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_Tint
62 F77_FUNC (sqrinc, 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_Tint
69 F77_FUNC (sqrdec, 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_Tint
76 F77_FUNC (sqrinr, 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_Tint
82 F77_FUNC (sqrder, 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_Tint
88 F77_FUNC (sqrshc, 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
98FloatQR::FloatQR (const FloatMatrix& a, qr_type_t qr_type)
99{
100 init (a, qr_type);
101}
102
103void
104FloatQR::init (const FloatMatrix& a, qr_type_t qr_type)
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)octave_local_buffer<float> _buffer_tau (min_mn); float *
tau = _buffer_tau
;
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,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"
, "sgeqrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sgeqrf_ (m, n, afact.fortran_vec (), m, tau, &rlwork,
-1, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
123 &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"
, "sgeqrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sgeqrf_ (m, n, afact.fortran_vec (), m, tau, &rlwork,
-1, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
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)octave_local_buffer<float> _buffer_work (lwork); float *
work = _buffer_work
;
129 F77_XFCN (sgeqrf, SGEQRF, (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"
, "sgeqrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sgeqrf_ (m, n, afact.fortran_vec (), m, tau, work, lwork,
info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
130 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"
, "sgeqrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sgeqrf_ (m, n, afact.fortran_vec (), m, tau, work, lwork,
info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
131 }
132
133 form (n, afact, tau, qr_type);
134}
135
136void FloatQR::form (octave_idx_type n, FloatMatrix& afact,
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)
1
Assuming 'qr_type' is not equal to qr_type_raw
2
Taking false branch
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)
3
Taking false branch
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++)
4
Assuming 'j' is < 'm'
5
Loop condition is true. Entering loop body
8
Loop condition is false. Execution continues on line 182
177 for (octave_idx_type i = j + 1; i < m; i++)
6
Assuming 'i' is >= 'm'
7
Loop condition is false. Execution continues on line 176
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)
9
Taking true branch
187 {
188 octave_idx_type k = q.columns ();
189 // workspace query.
190 float rlwork;
10
'rlwork' declared without an initial value
191 F77_XFCN (sorgqr, SORGQR, (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"
, "sorgqr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sorgqr_ (m, k, min_mn, q.fortran_vec (), m, tau, &rlwork
, -1, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
192 &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"
, "sorgqr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sorgqr_ (m, k, min_mn, q.fortran_vec (), m, tau, &rlwork
, -1, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
193
194 // allocate buffer and do the job.
195 octave_idx_type lwork = rlwork;
11
Assigned value is garbage or undefined
196 lwork = std::max (lwork, static_cast<octave_idx_type> (1));
197 OCTAVE_LOCAL_BUFFER (float, work, lwork)octave_local_buffer<float> _buffer_work (lwork); float *
work = _buffer_work
;
198 F77_XFCN (sorgqr, SORGQR, (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"
, "sorgqr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sorgqr_ (m, k, min_mn, q.fortran_vec (), m, tau, work, lwork
, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
199 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"
, "sorgqr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sorgqr_ (m, k, min_mn, q.fortran_vec (), m, tau, work, lwork
, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
200 }
201 }
202}
203
204#ifdef HAVE_QRUPDATE1
205
206void
207FloatQR::update (const FloatColumnVector& u, const FloatColumnVector& v)
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)octave_local_buffer<float> _buffer_w (2*k); float *w = _buffer_w;
217 F77_XFCN (sqr1up, SQR1UP, (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"
, "sqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqr1up_ (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)
218 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"
, "sqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqr1up_ (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)
219 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"
, "sqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqr1up_ (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 }
221 else
222 (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
223}
224
225void
226FloatQR::update (const FloatMatrix& u, const FloatMatrix& v)
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)octave_local_buffer<float> _buffer_w (2*k); float *w = _buffer_w;
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 (),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"
, "sqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqr1up_ (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)
239 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"
, "sqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqr1up_ (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)
240 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"
, "sqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqr1up_ (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 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"
, "sqr1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqr1up_ (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 }
243 }
244 else
245 (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
246}
247
248void
249FloatQR::insert_col (const FloatColumnVector& u, octave_idx_type j)
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)octave_local_buffer<float> _buffer_w (k); float *w = _buffer_w;
273 F77_XFCN (sqrinc, SQRINC, (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"
, "sqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrinc_ (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)
274 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"
, "sqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrinc_ (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)
275 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"
, "sqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrinc_ (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 }
277}
278
279void
280FloatQR::insert_col (const FloatMatrix& u, const Array<octave_idx_type>& j)
281{
282 octave_idx_type m = q.rows ();
283 octave_idx_type n = r.columns ();
284 octave_idx_type k = q.columns ();
285
286 Array<octave_idx_type> jsi;
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)octave_local_buffer<float> _buffer_w (kmax); float *w =
_buffer_w
;
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),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"
, "sqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrinc_ (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)
318 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"
, "sqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrinc_ (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)
319 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"
, "sqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrinc_ (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 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"
, "sqrinc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrinc_ (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 }
322 }
323}
324
325void
326FloatQR::delete_col (octave_idx_type j)
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)octave_local_buffer<float> _buffer_w (k); float *w = _buffer_w;
337 F77_XFCN (sqrdec, SQRDEC, (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"
, "sqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrdec_ (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)
338 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"
, "sqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrdec_ (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)
;
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
352void
353FloatQR::delete_col (const Array<octave_idx_type>& j)
354{
355 octave_idx_type m = q.rows ();
356 octave_idx_type n = r.columns ();
357 octave_idx_type k = q.columns ();
358
359 Array<octave_idx_type> jsi;
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)octave_local_buffer<float> _buffer_w (k); float *w = _buffer_w;
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,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"
, "sqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrdec_ (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)
377 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"
, "sqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrdec_ (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)
378 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"
, "sqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrdec_ (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 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"
, "sqrdec_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrdec_ (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 }
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
394void
395FloatQR::insert_row (const FloatRowVector& u, octave_idx_type j)
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)octave_local_buffer<float> _buffer_w (k); float *w = _buffer_w;
411 F77_XFCN (sqrinr, SQRINR, (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"
, "sqrinr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrinr_ (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)
412 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"
, "sqrinr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrinr_ (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)
413 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"
, "sqrinr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrinr_ (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
415 }
416}
417
418void
419FloatQR::delete_row (octave_idx_type j)
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)octave_local_buffer<float> _buffer_w (2*m); float *w = _buffer_w;
431 F77_XFCN (sqrder, SQRDER, (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"
, "sqrder_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrder_ (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)
432 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"
, "sqrder_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrder_ (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)
433 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"
, "sqrder_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrder_ (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
435 q.resize (m - 1, m - 1);
436 r.resize (m - 1, n);
437 }
438}
439
440void
441FloatQR::shift_cols (octave_idx_type i, octave_idx_type j)
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)octave_local_buffer<float> _buffer_w (2*k); float *w = _buffer_w;
452 F77_XFCN (sqrshc, SQRSHC, (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"
, "sqrshc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrshc_ (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)
453 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"
, "sqrshc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrshc_ (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)
454 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"
, "sqrshc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrshc_ (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 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"
, "sqrshc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; sqrshc_ (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 }
457}
458
459#else
460
461// Replacement update methods.
462
463void
464FloatQR::update (const FloatColumnVector& u, const FloatColumnVector& v)
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
479void
480FloatQR::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
495static
496FloatMatrix 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
508static
509FloatMatrix 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
521static
522FloatMatrix 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
529static
530FloatMatrix 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
537static
538FloatMatrix shift_cols (const FloatMatrix& a,
539 octave_idx_type i, octave_idx_type j)
540{
541 octave_idx_type n = a.columns ();
542 Array<octave_idx_type> p (n);
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
558void
559FloatQR::insert_col (const FloatColumnVector& u, octave_idx_type j)
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
576void
577FloatQR::insert_col (const FloatMatrix& u, const Array<octave_idx_type>& j)
578{
579 warn_qrupdate_once ();
580
581 octave_idx_type m = q.rows ();
582 octave_idx_type n = r.columns ();
583
584 Array<octave_idx_type> jsi;
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
606void
607FloatQR::delete_col (octave_idx_type j)
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
622void
623FloatQR::delete_col (const Array<octave_idx_type>& j)
624{
625 warn_qrupdate_once ();
626
627 octave_idx_type m = q.rows ();
628 octave_idx_type n = r.columns ();
629
630 Array<octave_idx_type> jsi;
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
650void
651FloatQR::insert_row (const FloatRowVector& u, octave_idx_type j)
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
668void
669FloatQR::delete_row (octave_idx_type j)
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
686void
687FloatQR::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
702#endif