File: | liboctave/numeric/floatQR.cc |
Location: | line 195, column 11 |
Description: | Assigned value is garbage or undefined |
1 | /* | |||
2 | ||||
3 | Copyright (C) 1994-2013 John W. Eaton | |||
4 | Copyright (C) 2008-2009 Jaroslav Hajek | |||
5 | Copyright (C) 2009 VZLU Prague | |||
6 | ||||
7 | This file is part of Octave. | |||
8 | ||||
9 | Octave is free software; you can redistribute it and/or modify it | |||
10 | under the terms of the GNU General Public License as published by the | |||
11 | Free Software Foundation; either version 3 of the License, or (at your | |||
12 | option) any later version. | |||
13 | ||||
14 | Octave is distributed in the hope that it will be useful, but WITHOUT | |||
15 | ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |||
16 | FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |||
17 | for more details. | |||
18 | ||||
19 | You should have received a copy of the GNU General Public License | |||
20 | along with Octave; see the file COPYING. If not, see | |||
21 | <http://www.gnu.org/licenses/>. | |||
22 | ||||
23 | */ | |||
24 | ||||
25 | #ifdef HAVE_CONFIG_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 | ||||
38 | template class base_qr<FloatMatrix>; | |||
39 | ||||
40 | extern "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 | ||||
98 | FloatQR::FloatQR (const FloatMatrix& a, qr_type_t qr_type) | |||
99 | { | |||
100 | init (a, qr_type); | |||
101 | } | |||
102 | ||||
103 | void | |||
104 | FloatQR::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 | ||||
136 | void 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) | |||
| ||||
143 | { | |||
144 | for (octave_idx_type j = 0; j < min_mn; j++) | |||
145 | { | |||
146 | octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; | |||
147 | for (octave_idx_type i = limit + 1; i < m; i++) | |||
148 | afact.elem (i, j) *= tau[j]; | |||
149 | } | |||
150 | ||||
151 | r = afact; | |||
152 | } | |||
153 | else | |||
154 | { | |||
155 | // Attempt to minimize copying. | |||
156 | if (m >= n) | |||
157 | { | |||
158 | // afact will become q. | |||
159 | q = afact; | |||
160 | octave_idx_type k = qr_type == qr_type_economy ? n : m; | |||
161 | r = FloatMatrix (k, n); | |||
162 | for (octave_idx_type j = 0; j < n; j++) | |||
163 | { | |||
164 | octave_idx_type i = 0; | |||
165 | for (; i <= j; i++) | |||
166 | r.xelem (i, j) = afact.xelem (i, j); | |||
167 | for (; i < k; i++) | |||
168 | r.xelem (i, j) = 0; | |||
169 | } | |||
170 | afact = FloatMatrix (); // optimize memory | |||
171 | } | |||
172 | else | |||
173 | { | |||
174 | // afact will become r. | |||
175 | q = FloatMatrix (m, m); | |||
176 | for (octave_idx_type j = 0; j < m; j++) | |||
177 | for (octave_idx_type i = j + 1; i < m; i++) | |||
178 | { | |||
179 | q.xelem (i, j) = afact.xelem (i, j); | |||
180 | afact.xelem (i, j) = 0; | |||
181 | } | |||
182 | r = afact; | |||
183 | } | |||
184 | ||||
185 | ||||
186 | if (m > 0) | |||
187 | { | |||
188 | octave_idx_type k = q.columns (); | |||
189 | // workspace query. | |||
190 | float rlwork; | |||
191 | F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,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; | |||
| ||||
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 | ||||
206 | void | |||
207 | FloatQR::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 | ||||
225 | void | |||
226 | FloatQR::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 | ||||
248 | void | |||
249 | FloatQR::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 | ||||
279 | void | |||
280 | FloatQR::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 | ||||
325 | void | |||
326 | FloatQR::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 | ||||
352 | void | |||
353 | FloatQR::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 | ||||
394 | void | |||
395 | FloatQR::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 | ||||
418 | void | |||
419 | FloatQR::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 | ||||
440 | void | |||
441 | FloatQR::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 | ||||
463 | void | |||
464 | FloatQR::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 | ||||
479 | void | |||
480 | FloatQR::update (const FloatMatrix& u, const FloatMatrix& v) | |||
481 | { | |||
482 | warn_qrupdate_once (); | |||
483 | ||||
484 | octave_idx_type m = q.rows (); | |||
485 | octave_idx_type n = r.columns (); | |||
486 | ||||
487 | if (u.rows () == m && v.rows () == n && u.cols () == v.cols ()) | |||
488 | { | |||
489 | init (q*r + u * v.transpose (), get_type ()); | |||
490 | } | |||
491 | else | |||
492 | (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch"); | |||
493 | } | |||
494 | ||||
495 | static | |||
496 | FloatMatrix insert_col (const FloatMatrix& a, octave_idx_type i, | |||
497 | const FloatColumnVector& x) | |||
498 | { | |||
499 | FloatMatrix retval (a.rows (), a.columns () + 1); | |||
500 | retval.assign (idx_vector::colon, idx_vector (0, i), | |||
501 | a.index (idx_vector::colon, idx_vector (0, i))); | |||
502 | retval.assign (idx_vector::colon, idx_vector (i), x); | |||
503 | retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()), | |||
504 | a.index (idx_vector::colon, idx_vector (i, a.columns ()))); | |||
505 | return retval; | |||
506 | } | |||
507 | ||||
508 | static | |||
509 | FloatMatrix insert_row (const FloatMatrix& a, octave_idx_type i, | |||
510 | const FloatRowVector& x) | |||
511 | { | |||
512 | FloatMatrix retval (a.rows () + 1, a.columns ()); | |||
513 | retval.assign (idx_vector (0, i), idx_vector::colon, | |||
514 | a.index (idx_vector (0, i), idx_vector::colon)); | |||
515 | retval.assign (idx_vector (i), idx_vector::colon, x); | |||
516 | retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon, | |||
517 | a.index (idx_vector (i, a.rows ()), idx_vector::colon)); | |||
518 | return retval; | |||
519 | } | |||
520 | ||||
521 | static | |||
522 | FloatMatrix delete_col (const FloatMatrix& a, octave_idx_type i) | |||
523 | { | |||
524 | FloatMatrix retval = a; | |||
525 | retval.delete_elements (1, idx_vector (i)); | |||
526 | return retval; | |||
527 | } | |||
528 | ||||
529 | static | |||
530 | FloatMatrix delete_row (const FloatMatrix& a, octave_idx_type i) | |||
531 | { | |||
532 | FloatMatrix retval = a; | |||
533 | retval.delete_elements (0, idx_vector (i)); | |||
534 | return retval; | |||
535 | } | |||
536 | ||||
537 | static | |||
538 | FloatMatrix shift_cols (const FloatMatrix& a, | |||
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 | ||||
558 | void | |||
559 | FloatQR::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 | ||||
576 | void | |||
577 | FloatQR::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 | ||||
606 | void | |||
607 | FloatQR::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 | ||||
622 | void | |||
623 | FloatQR::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 | ||||
650 | void | |||
651 | FloatQR::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 | ||||
668 | void | |||
669 | FloatQR::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 | ||||
686 | void | |||
687 | FloatQR::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 |