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