Mercurial > jwe > octave

comparison liboctave/numeric/dbleCHOL.cc @ 21136:7cac4e7458f2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: clean up code around calls to current_liboctave_error_handler. Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code. * Array-util.cc, Array.cc, CColVector.cc, CDiagMatrix.cc, CMatrix.cc, CNDArray.cc, CRowVector.cc, CSparse.cc, DiagArray2.cc, MArray.cc, PermMatrix.cc, Sparse.cc, Sparse.h, chMatrix.cc, chNDArray.cc, dColVector.cc, dDiagMatrix.cc, dMatrix.cc, dNDArray.cc, dRowVector.cc, dSparse.cc, fCColVector.cc, fCDiagMatrix.cc, fCMatrix.cc, fCNDArray.cc, fCRowVector.cc, fColVector.cc, fDiagMatrix.cc, fMatrix.cc, fNDArray.cc, fRowVector.cc, idx-vector.cc, CmplxAEPBAL.cc, CmplxCHOL.cc, CmplxGEPBAL.cc, CmplxHESS.cc, CmplxLU.cc, CmplxQR.cc, CmplxSCHUR.cc, CmplxSVD.cc, DASPK.cc, EIG.cc, LSODE.cc, Quad.cc, SparseCmplxCHOL.cc, SparseCmplxLU.cc, SparseCmplxQR.cc, SparseQR.cc, SparsedbleCHOL.cc, SparsedbleLU.cc, base-lu.cc, bsxfun-defs.cc, dbleAEPBAL.cc, dbleCHOL.cc, dbleGEPBAL.cc, dbleHESS.cc, dbleLU.cc, dbleQR.cc, dbleSCHUR.cc, dbleSVD.cc, eigs-base.cc, fCmplxAEPBAL.cc, fCmplxCHOL.cc, fCmplxLU.cc, fCmplxQR.cc, fCmplxSCHUR.cc, fEIG.cc, floatAEPBAL.cc, floatCHOL.cc, floatGEPBAL.cc, floatHESS.cc, floatLU.cc, floatQR.cc, floatSCHUR.cc, floatSVD.cc, lo-specfun.cc, oct-fftw.cc, oct-rand.cc, oct-spparms.cc, sparse-base-chol.cc, sparse-dmsolve.cc, file-ops.cc, lo-sysdep.cc, mach-info.cc, oct-env.cc, oct-syscalls.cc, cmd-edit.cc, cmd-hist.cc, data-conv.cc, lo-ieee.cc, lo-regexp.cc, oct-base64.cc, oct-shlib.cc, pathsearch.cc, singleton-cleanup.cc, sparse-util.cc, unwind-prot.cc: Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code.
author Rik <rik@octave.org>
date Sat, 23 Jan 2016 13:52:03 -0800
parents ab705b42cfd8
children f7121e111991
comparison
equal deleted inserted replaced
21135:95da3bc8a281 21136:7cac4e7458f2
91 { 91 {
92 octave_idx_type a_nr = a.rows (); 92 octave_idx_type a_nr = a.rows ();
93 octave_idx_type a_nc = a.cols (); 93 octave_idx_type a_nc = a.cols ();
94 94
95 if (a_nr != a_nc) 95 if (a_nr != a_nc)
96 { 96 (*current_liboctave_error_handler) ("CHOL requires square matrix");
97 (*current_liboctave_error_handler) ("CHOL requires square matrix");
98 return -1;
99 }
100 97
101 octave_idx_type n = a_nc; 98 octave_idx_type n = a_nc;
102 octave_idx_type info; 99 octave_idx_type info;
103 100
104 is_upper = upper; 101 is_upper = upper;
168 Matrix retval; 165 Matrix retval;
169 166
170 octave_idx_type r_nr = r.rows (); 167 octave_idx_type r_nr = r.rows ();
171 octave_idx_type r_nc = r.cols (); 168 octave_idx_type r_nc = r.cols ();
172 169
173 if (r_nr == r_nc) 170 if (r_nr != r_nc)
174 { 171 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
175 octave_idx_type n = r_nc; 172
176 octave_idx_type info = 0; 173 octave_idx_type n = r_nc;
177 174 octave_idx_type info = 0;
178 Matrix tmp = r; 175
179 double *v = tmp.fortran_vec (); 176 Matrix tmp = r;
180 177 double *v = tmp.fortran_vec ();
181 if (info == 0) 178
179 if (info == 0)
180 {
181 if (is_upper)
182 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
183 v, n, info
184 F77_CHAR_ARG_LEN (1)));
185 else
186 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
187 v, n, info
188 F77_CHAR_ARG_LEN (1)));
189
190 // If someone thinks of a more graceful way of doing this (or
191 // faster for that matter :-)), please let me know!
192
193 if (n > 1)
182 { 194 {
183 if (is_upper) 195 if (is_upper)
184 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, 196 for (octave_idx_type j = 0; j < r_nc; j++)
185 v, n, info 197 for (octave_idx_type i = j+1; i < r_nr; i++)
186 F77_CHAR_ARG_LEN (1))); 198 tmp.xelem (i, j) = tmp.xelem (j, i);
187 else 199 else
188 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n, 200 for (octave_idx_type j = 0; j < r_nc; j++)
189 v, n, info 201 for (octave_idx_type i = j+1; i < r_nr; i++)
190 F77_CHAR_ARG_LEN (1))); 202 tmp.xelem (j, i) = tmp.xelem (i, j);
191
192 // If someone thinks of a more graceful way of doing this (or
193 // faster for that matter :-)), please let me know!
194
195 if (n > 1)
196 {
197 if (is_upper)
198 for (octave_idx_type j = 0; j < r_nc; j++)
199 for (octave_idx_type i = j+1; i < r_nr; i++)
200 tmp.xelem (i, j) = tmp.xelem (j, i);
201 else
202 for (octave_idx_type j = 0; j < r_nc; j++)
203 for (octave_idx_type i = j+1; i < r_nr; i++)
204 tmp.xelem (j, i) = tmp.xelem (i, j);
205 }
206
207 retval = tmp;
208 } 203 }
209 } 204
210 else 205 retval = tmp;
211 (*current_liboctave_error_handler) ("chol2inv requires square matrix"); 206 }
212 207
213 return retval; 208 return retval;
214 } 209 }
215 210
216 // Compute the inverse of a matrix using the Cholesky factorization. 211 // Compute the inverse of a matrix using the Cholesky factorization.
221 } 216 }
222 217
223 void 218 void
224 CHOL::set (const Matrix& R) 219 CHOL::set (const Matrix& R)
225 { 220 {
226 if (R.is_square ()) 221 if (! R.is_square ())
227 chol_mat = R;
228 else
229 (*current_liboctave_error_handler) ("CHOL requires square matrix"); 222 (*current_liboctave_error_handler) ("CHOL requires square matrix");
223
224 chol_mat = R;
230 } 225 }
231 226
232 #ifdef HAVE_QRUPDATE 227 #ifdef HAVE_QRUPDATE
233 228
234 void 229 void
235 CHOL::update (const ColumnVector& u) 230 CHOL::update (const ColumnVector& u)
236 { 231 {
237 octave_idx_type n = chol_mat.rows (); 232 octave_idx_type n = chol_mat.rows ();
238 233
239 if (u.numel () == n) 234 if (u.numel () != n)
240 {
241 ColumnVector utmp = u;
242
243 OCTAVE_LOCAL_BUFFER (double, w, n);
244
245 F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
246 utmp.fortran_vec (), w));
247 }
248 else
249 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); 235 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
236
237 ColumnVector utmp = u;
238
239 OCTAVE_LOCAL_BUFFER (double, w, n);
240
241 F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
242 utmp.fortran_vec (), w));
250 } 243 }
251 244
252 octave_idx_type 245 octave_idx_type
253 CHOL::downdate (const ColumnVector& u) 246 CHOL::downdate (const ColumnVector& u)
254 { 247 {
255 octave_idx_type info = -1; 248 octave_idx_type info = -1;
256 249
257 octave_idx_type n = chol_mat.rows (); 250 octave_idx_type n = chol_mat.rows ();
258 251
259 if (u.numel () == n) 252 if (u.numel () != n)
260 {
261 ColumnVector utmp = u;
262
263 OCTAVE_LOCAL_BUFFER (double, w, n);
264
265 F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
266 utmp.fortran_vec (), w, info));
267 }
268 else
269 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); 253 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
254
255 ColumnVector utmp = u;
256
257 OCTAVE_LOCAL_BUFFER (double, w, n);
258
259 F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
260 utmp.fortran_vec (), w, info));
270 261
271 return info; 262 return info;
272 } 263 }
273 264
274 octave_idx_type 265 octave_idx_type
278 269
279 octave_idx_type n = chol_mat.rows (); 270 octave_idx_type n = chol_mat.rows ();
280 271
281 if (u.numel () != n + 1) 272 if (u.numel () != n + 1)
282 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); 273 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
283 else if (j < 0 || j > n) 274 if (j < 0 || j > n)
284 (*current_liboctave_error_handler) ("cholinsert: index out of range"); 275 (*current_liboctave_error_handler) ("cholinsert: index out of range");
285 else 276
286 { 277 ColumnVector utmp = u;
287 ColumnVector utmp = u; 278
288 279 OCTAVE_LOCAL_BUFFER (double, w, n);
289 OCTAVE_LOCAL_BUFFER (double, w, n); 280
290 281 chol_mat.resize (n+1, n+1);
291 chol_mat.resize (n+1, n+1); 282
292 283 F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
293 F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), 284 j + 1, utmp.fortran_vec (), w, info));
294 j + 1, utmp.fortran_vec (), w, info));
295 }
296 285
297 return info; 286 return info;
298 } 287 }
299 288
300 void 289 void
302 { 291 {
303 octave_idx_type n = chol_mat.rows (); 292 octave_idx_type n = chol_mat.rows ();
304 293
305 if (j < 0 || j > n-1) 294 if (j < 0 || j > n-1)
306 (*current_liboctave_error_handler) ("choldelete: index out of range"); 295 (*current_liboctave_error_handler) ("choldelete: index out of range");
307 else 296
308 { 297 OCTAVE_LOCAL_BUFFER (double, w, n);
309 OCTAVE_LOCAL_BUFFER (double, w, n); 298
310 299 F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
311 F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), 300 j + 1, w));
312 j + 1, w)); 301
313 302 chol_mat.resize (n-1, n-1);
314 chol_mat.resize (n-1, n-1);
315 }
316 } 303 }
317 304
318 void 305 void
319 CHOL::shift_sym (octave_idx_type i, octave_idx_type j) 306 CHOL::shift_sym (octave_idx_type i, octave_idx_type j)
320 { 307 {
321 octave_idx_type n = chol_mat.rows (); 308 octave_idx_type n = chol_mat.rows ();
322 309
323 if (i < 0 || i > n-1 || j < 0 || j > n-1) 310 if (i < 0 || i > n-1 || j < 0 || j > n-1)
324 (*current_liboctave_error_handler) ("cholshift: index out of range"); 311 (*current_liboctave_error_handler) ("cholshift: index out of range");
325 else 312
326 { 313 OCTAVE_LOCAL_BUFFER (double, w, 2*n);
327 OCTAVE_LOCAL_BUFFER (double, w, 2*n); 314
328 315 F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
329 F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), 316 i + 1, j + 1, w));
330 i + 1, j + 1, w));
331 }
332 } 317 }
333 318
334 #else 319 #else
335 320
336 void 321 void
338 { 323 {
339 warn_qrupdate_once (); 324 warn_qrupdate_once ();
340 325
341 octave_idx_type n = chol_mat.rows (); 326 octave_idx_type n = chol_mat.rows ();
342 327
343 if (u.numel () == n) 328 if (u.numel () != n)
344 {
345 init (chol_mat.transpose () * chol_mat
346 + Matrix (u) * Matrix (u).transpose (), true, false);
347 }
348 else
349 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); 329 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
330
331 init (chol_mat.transpose () * chol_mat
332 + Matrix (u) * Matrix (u).transpose (), true, false);
350 } 333 }
351 334
352 static bool 335 static bool
353 singular (const Matrix& a) 336 singular (const Matrix& a)
354 { 337 {
364 347
365 octave_idx_type info = -1; 348 octave_idx_type info = -1;
366 349
367 octave_idx_type n = chol_mat.rows (); 350 octave_idx_type n = chol_mat.rows ();
368 351
369 if (u.numel () == n) 352 if (u.numel () != n)
370 { 353 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
371 if (singular (chol_mat)) 354
372 info = 2; 355 if (singular (chol_mat))
373 else 356 info = 2;
374 {
375 info = init (chol_mat.transpose () * chol_mat
376 - Matrix (u) * Matrix (u).transpose (), true, false);
377 if (info) info = 1;
378 }
379 }
380 else 357 else
381 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); 358 {
359 info = init (chol_mat.transpose () * chol_mat
360 - Matrix (u) * Matrix (u).transpose (), true, false);
361 if (info) info = 1;
362 }
382 363
383 return info; 364 return info;
384 } 365 }
385 366
386 octave_idx_type 367 octave_idx_type
392 373
393 octave_idx_type n = chol_mat.rows (); 374 octave_idx_type n = chol_mat.rows ();
394 375
395 if (u.numel () != n + 1) 376 if (u.numel () != n + 1)
396 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); 377 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
397 else if (j < 0 || j > n) 378 if (j < 0 || j > n)
398 (*current_liboctave_error_handler) ("cholinsert: index out of range"); 379 (*current_liboctave_error_handler) ("cholinsert: index out of range");
380
381 if (singular (chol_mat))
382 info = 2;
399 else 383 else
400 { 384 {
401 if (singular (chol_mat)) 385 Matrix a = chol_mat.transpose () * chol_mat;
402 info = 2; 386 Matrix a1 (n+1, n+1);
403 else 387 for (octave_idx_type k = 0; k < n+1; k++)
404 { 388 for (octave_idx_type l = 0; l < n+1; l++)
405 Matrix a = chol_mat.transpose () * chol_mat; 389 {
406 Matrix a1 (n+1, n+1); 390 if (l == j)
407 for (octave_idx_type k = 0; k < n+1; k++) 391 a1(k, l) = u(k);
408 for (octave_idx_type l = 0; l < n+1; l++) 392 else if (k == j)
409 { 393 a1(k, l) = u(l);
410 if (l == j) 394 else
411 a1(k, l) = u(k); 395 a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
412 else if (k == j) 396 }
413 a1(k, l) = u(l); 397 info = init (a1, true, false);
414 else 398 if (info) info = 1;
415 a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
416 }
417 info = init (a1, true, false);
418 if (info) info = 1;
419 }
420 } 399 }
421 400
422 return info; 401 return info;
423 } 402 }
424 403
429 408
430 octave_idx_type n = chol_mat.rows (); 409 octave_idx_type n = chol_mat.rows ();
431 410
432 if (j < 0 || j > n-1) 411 if (j < 0 || j > n-1)
433 (*current_liboctave_error_handler) ("choldelete: index out of range"); 412 (*current_liboctave_error_handler) ("choldelete: index out of range");
434 else 413
435 { 414 Matrix a = chol_mat.transpose () * chol_mat;
436 Matrix a = chol_mat.transpose () * chol_mat; 415 a.delete_elements (1, idx_vector (j));
437 a.delete_elements (1, idx_vector (j)); 416 a.delete_elements (0, idx_vector (j));
438 a.delete_elements (0, idx_vector (j)); 417 init (a, true, false);
439 init (a, true, false);
440 }
441 } 418 }
442 419
443 void 420 void
444 CHOL::shift_sym (octave_idx_type i, octave_idx_type j) 421 CHOL::shift_sym (octave_idx_type i, octave_idx_type j)
445 { 422 {
447 424
448 octave_idx_type n = chol_mat.rows (); 425 octave_idx_type n = chol_mat.rows ();
449 426
450 if (i < 0 || i > n-1 || j < 0 || j > n-1) 427 if (i < 0 || i > n-1 || j < 0 || j > n-1)
451 (*current_liboctave_error_handler) ("cholshift: index out of range"); 428 (*current_liboctave_error_handler) ("cholshift: index out of range");
452 else 429
453 { 430 Matrix a = chol_mat.transpose () * chol_mat;
454 Matrix a = chol_mat.transpose () * chol_mat; 431 Array<octave_idx_type> p (dim_vector (n, 1));
455 Array<octave_idx_type> p (dim_vector (n, 1)); 432 for (octave_idx_type k = 0; k < n; k++) p(k) = k;
456 for (octave_idx_type k = 0; k < n; k++) p(k) = k; 433 if (i < j)
457 if (i < j) 434 {
458 { 435 for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
459 for (octave_idx_type k = i; k < j; k++) p(k) = k+1; 436 p(j) = i;
460 p(j) = i; 437 }
461 } 438 else if (j < i)
462 else if (j < i) 439 {
463 { 440 p(j) = i;
464 p(j) = i; 441 for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
465 for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; 442 }
466 } 443
467 444 init (a.index (idx_vector (p), idx_vector (p)), true, false);
468 init (a.index (idx_vector (p), idx_vector (p)), true, false);
469 }
470 } 445 }
471 446
472 #endif 447 #endif
473 448
474 Matrix 449 Matrix