Mercurial > hg > octave-nkf

comparison liboctave/dbleCHOL.cc @ 8547:d66c9b6e506a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
imported patch qrupdate.diff
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 20 Jan 2009 21:16:42 +0100
parents 25bc2d31e1bf
children a6edd5c23cb5
comparison
equal deleted inserted replaced
8546:3d8a914c580e 8547:d66c9b6e506a
1 /* 1 /*
2 2
3 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 3 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
4 John W. Eaton 4 John W. Eaton
5 Copyright (C) 2008, 2009 Jaroslav Hajek
5 6
6 This file is part of Octave. 7 This file is part of Octave.
7 8
8 Octave is free software; you can redistribute it and/or modify it 9 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the 10 under the terms of the GNU General Public License as published by the
18 You should have received a copy of the GNU General Public License 19 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see 20 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>. 21 <http://www.gnu.org/licenses/>.
21 22
22 */ 23 */
23
24 // updating/downdating by Jaroslav Hajek 2008
25 24
26 #ifdef HAVE_CONFIG_H 25 #ifdef HAVE_CONFIG_H
27 #include <config.h> 26 #include <config.h>
28 #endif 27 #endif
29 28
50 F77_RET_T 49 F77_RET_T
51 F77_FUNC (dpocon, DPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, 50 F77_FUNC (dpocon, DPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
52 double*, const octave_idx_type&, const double&, 51 double*, const octave_idx_type&, const double&,
53 double&, double*, octave_idx_type*, 52 double&, double*, octave_idx_type*,
54 octave_idx_type& F77_CHAR_ARG_LEN_DECL); 53 octave_idx_type& F77_CHAR_ARG_LEN_DECL);
55 F77_RET_T 54 #ifdef HAVE_QRUPDATE
56 F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*, double*, double*); 55
57 56 F77_RET_T
58 F77_RET_T 57 F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*, const octave_idx_type&,
59 F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, double*, double*, 58 double*, double*);
59
60 F77_RET_T
61 F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, const octave_idx_type&,
62 double*, double*, octave_idx_type&);
63
64 F77_RET_T
65 F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, double*, const octave_idx_type&,
66 const octave_idx_type&, double*, double*,
60 octave_idx_type&); 67 octave_idx_type&);
61 68
62 F77_RET_T 69 F77_RET_T
63 F77_FUNC (dqrshc, DQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 70 F77_FUNC (dchdex, DCHDEX) (const octave_idx_type&, double*, const octave_idx_type&,
64 double*, double*, const octave_idx_type&, const octave_idx_type&); 71 const octave_idx_type&, double*);
65 72
66 F77_RET_T 73 F77_RET_T
67 F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, const double*, double*, const octave_idx_type&, 74 F77_FUNC (dchshx, DCHSHX) (const octave_idx_type&, double*, const octave_idx_type&,
68 const double*, octave_idx_type&); 75 const octave_idx_type&, const octave_idx_type&,
69 76 double*);
70 F77_RET_T 77 #endif
71 F77_FUNC (dchdex, DCHDEX) (const octave_idx_type&, const double*, double*, const octave_idx_type&);
72 } 78 }
73 79
74 octave_idx_type 80 octave_idx_type
75 CHOL::init (const Matrix& a, bool calc_cond) 81 CHOL::init (const Matrix& a, bool calc_cond)
76 { 82 {
184 chol_mat = R; 190 chol_mat = R;
185 else 191 else
186 (*current_liboctave_error_handler) ("CHOL requires square matrix"); 192 (*current_liboctave_error_handler) ("CHOL requires square matrix");
187 } 193 }
188 194
195 #ifdef HAVE_QRUPDATE
196
189 void 197 void
190 CHOL::update (const Matrix& u) 198 CHOL::update (const ColumnVector& u)
191 { 199 {
192 octave_idx_type n = chol_mat.rows (); 200 octave_idx_type n = chol_mat.rows ();
193 201
194 if (u.length () == n) 202 if (u.length () == n)
195 { 203 {
196 Matrix tmp = u; 204 ColumnVector utmp = u;
197 205
198 OCTAVE_LOCAL_BUFFER (double, w, n); 206 OCTAVE_LOCAL_BUFFER (double, w, n);
199 207
200 F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), 208 F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
201 tmp.fortran_vec (), w)); 209 utmp.fortran_vec (), w));
202 } 210 }
203 else 211 else
204 (*current_liboctave_error_handler) ("CHOL update dimension mismatch"); 212 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
205 } 213 }
206 214
207 octave_idx_type 215 octave_idx_type
208 CHOL::downdate (const Matrix& u) 216 CHOL::downdate (const ColumnVector& u)
209 { 217 {
210 octave_idx_type info = -1; 218 octave_idx_type info = -1;
211 219
212 octave_idx_type n = chol_mat.rows (); 220 octave_idx_type n = chol_mat.rows ();
213 221
214 if (u.length () == n) 222 if (u.length () == n)
215 { 223 {
216 Matrix tmp = u; 224 ColumnVector utmp = u;
217 225
218 OCTAVE_LOCAL_BUFFER (double, w, n); 226 OCTAVE_LOCAL_BUFFER (double, w, n);
219 227
220 F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), 228 F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
221 tmp.fortran_vec (), w, info)); 229 utmp.fortran_vec (), w, info));
222 } 230 }
223 else 231 else
224 (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch"); 232 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
225 233
226 return info; 234 return info;
227 } 235 }
228 236
229 octave_idx_type 237 octave_idx_type
230 CHOL::insert_sym (const Matrix& u, octave_idx_type j) 238 CHOL::insert_sym (const ColumnVector& u, octave_idx_type j)
231 { 239 {
232 octave_idx_type info = -1; 240 octave_idx_type info = -1;
233 241
234 octave_idx_type n = chol_mat.rows (); 242 octave_idx_type n = chol_mat.rows ();
235 243
236 if (u.length () != n+1) 244 if (u.length () != n + 1)
237 (*current_liboctave_error_handler) ("CHOL insert dimension mismatch"); 245 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
238 else if (j < 0 || j > n) 246 else if (j < 0 || j > n)
239 (*current_liboctave_error_handler) ("CHOL insert index out of range"); 247 (*current_liboctave_error_handler) ("cholinsert: index out of range");
240 else 248 else
241 { 249 {
242 Matrix chol_mat1 (n+1, n+1); 250 ColumnVector utmp = u;
243 251
244 F77_XFCN (dchinx, DCHINX, (n, chol_mat.data (), chol_mat1.fortran_vec (), 252 OCTAVE_LOCAL_BUFFER (double, w, n);
245 j+1, u.data (), info)); 253
246 254 chol_mat.resize (n+1, n+1);
247 chol_mat = chol_mat1; 255
256 F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
257 j + 1, utmp.fortran_vec (), w, info));
248 } 258 }
249 259
250 return info; 260 return info;
251 } 261 }
252 262
254 CHOL::delete_sym (octave_idx_type j) 264 CHOL::delete_sym (octave_idx_type j)
255 { 265 {
256 octave_idx_type n = chol_mat.rows (); 266 octave_idx_type n = chol_mat.rows ();
257 267
258 if (j < 0 || j > n-1) 268 if (j < 0 || j > n-1)
259 (*current_liboctave_error_handler) ("CHOL delete index out of range"); 269 (*current_liboctave_error_handler) ("choldelete: index out of range");
260 else 270 else
261 { 271 {
262 Matrix chol_mat1 (n-1, n-1); 272 OCTAVE_LOCAL_BUFFER (double, w, n);
263 273
264 F77_XFCN (dchdex, DCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1)); 274 F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
265 275 j + 1, w));
266 chol_mat = chol_mat1; 276
277 chol_mat.resize (n-1, n-1);
267 } 278 }
268 } 279 }
269 280
270 void 281 void
271 CHOL::shift_sym (octave_idx_type i, octave_idx_type j) 282 CHOL::shift_sym (octave_idx_type i, octave_idx_type j)
272 { 283 {
273 octave_idx_type n = chol_mat.rows (); 284 octave_idx_type n = chol_mat.rows ();
274 double dummy;
275 285
276 if (i < 0 || i > n-1 || j < 0 || j > n-1) 286 if (i < 0 || i > n-1 || j < 0 || j > n-1)
277 (*current_liboctave_error_handler) ("CHOL shift index out of range"); 287 (*current_liboctave_error_handler) ("cholshift: index out of range");
278 else 288 else
279 F77_XFCN (dqrshc, DQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1)); 289 {
280 } 290 OCTAVE_LOCAL_BUFFER (double, w, 2*n);
291
292 F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
293 i + 1, j + 1, w));
294 }
295 }
296
297 #endif
281 298
282 Matrix 299 Matrix
283 chol2inv (const Matrix& r) 300 chol2inv (const Matrix& r)
284 { 301 {
285 return chol2inv_internal (r); 302 return chol2inv_internal (r);