Mercurial > jwe > octave

comparison liboctave/numeric/SparseCmplxLU.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 bd1752782e56
children
comparison
equal deleted inserted replaced
21135:95da3bc8a281 21136:7cac4e7458f2
107 0, 0, 107 0, 0,
108 &Symbolic, control, info); 108 &Symbolic, control, info);
109 109
110 if (status < 0) 110 if (status < 0)
111 { 111 {
112 UMFPACK_ZNAME (report_status) (control, status);
113 UMFPACK_ZNAME (report_info) (control, info);
114
115 UMFPACK_ZNAME (free_symbolic) (&Symbolic);
116
112 (*current_liboctave_error_handler) 117 (*current_liboctave_error_handler)
113 ("SparseComplexLU::SparseComplexLU symbolic factorization failed"); 118 ("SparseComplexLU::SparseComplexLU symbolic factorization failed");
114
115 UMFPACK_ZNAME (report_status) (control, status);
116 UMFPACK_ZNAME (report_info) (control, info);
117
118 UMFPACK_ZNAME (free_symbolic) (&Symbolic);
119 } 119 }
120 else 120 else
121 { 121 {
122 UMFPACK_ZNAME (report_symbolic) (Symbolic, control); 122 UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
123 123
130 130
131 cond = Info (UMFPACK_RCOND); 131 cond = Info (UMFPACK_RCOND);
132 132
133 if (status < 0) 133 if (status < 0)
134 { 134 {
135 UMFPACK_ZNAME (report_status) (control, status);
136 UMFPACK_ZNAME (report_info) (control, info);
137
138 UMFPACK_ZNAME (free_numeric) (&Numeric);
139
135 (*current_liboctave_error_handler) 140 (*current_liboctave_error_handler)
136 ("SparseComplexLU::SparseComplexLU numeric factorization failed"); 141 ("SparseComplexLU::SparseComplexLU numeric factorization failed");
137
138 UMFPACK_ZNAME (report_status) (control, status);
139 UMFPACK_ZNAME (report_info) (control, info);
140
141 UMFPACK_ZNAME (free_numeric) (&Numeric);
142 } 142 }
143 else 143 else
144 { 144 {
145 UMFPACK_ZNAME (report_numeric) (Numeric, control); 145 UMFPACK_ZNAME (report_numeric) (Numeric, control);
146 146
148 status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz, &ignore1, 148 status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz, &ignore1,
149 &ignore2, &ignore3, Numeric); 149 &ignore2, &ignore3, Numeric);
150 150
151 if (status < 0) 151 if (status < 0)
152 { 152 {
153 UMFPACK_ZNAME (report_status) (control, status);
154 UMFPACK_ZNAME (report_info) (control, info);
155
156 UMFPACK_ZNAME (free_numeric) (&Numeric);
157
153 (*current_liboctave_error_handler) 158 (*current_liboctave_error_handler)
154 ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); 159 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
155
156 UMFPACK_ZNAME (report_status) (control, status);
157 UMFPACK_ZNAME (report_info) (control, info);
158
159 UMFPACK_ZNAME (free_numeric) (&Numeric);
160 } 160 }
161 else 161 else
162 { 162 {
163 octave_idx_type n_inner = (nr < nc ? nr : nc); 163 octave_idx_type n_inner = (nr < nc ? nr : nc);
164 164
207 207
208 UMFPACK_ZNAME (free_numeric) (&Numeric); 208 UMFPACK_ZNAME (free_numeric) (&Numeric);
209 209
210 if (status < 0) 210 if (status < 0)
211 { 211 {
212 UMFPACK_ZNAME (report_status) (control, status);
213
212 (*current_liboctave_error_handler) 214 (*current_liboctave_error_handler)
213 ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); 215 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
214
215 UMFPACK_ZNAME (report_status) (control, status);
216 } 216 }
217 else 217 else
218 { 218 {
219 Lfact = Lfact.transpose (); 219 Lfact = Lfact.transpose ();
220 220
254 { 254 {
255 #ifdef HAVE_UMFPACK 255 #ifdef HAVE_UMFPACK
256 if (milu) 256 if (milu)
257 (*current_liboctave_error_handler) 257 (*current_liboctave_error_handler)
258 ("Modified incomplete LU not implemented"); 258 ("Modified incomplete LU not implemented");
259 else 259
260 { 260 octave_idx_type nr = a.rows ();
261 octave_idx_type nr = a.rows (); 261 octave_idx_type nc = a.cols ();
262 octave_idx_type nc = a.cols (); 262
263 263 // Setup the control parameters
264 // Setup the control parameters 264 Matrix Control (UMFPACK_CONTROL, 1);
265 Matrix Control (UMFPACK_CONTROL, 1); 265 double *control = Control.fortran_vec ();
266 double *control = Control.fortran_vec (); 266 UMFPACK_ZNAME (defaults) (control);
267 UMFPACK_ZNAME (defaults) (control); 267
268 268 double tmp = octave_sparse_params::get_key ("spumoni");
269 double tmp = octave_sparse_params::get_key ("spumoni"); 269 if (! xisnan (tmp))
270 if (! xisnan (tmp)) 270 Control (UMFPACK_PRL) = tmp;
271 Control (UMFPACK_PRL) = tmp; 271 if (piv_thres.numel () == 2)
272 if (piv_thres.numel () == 2) 272 {
273 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
274 if (! xisnan (tmp))
275 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
276 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
277 if (! xisnan (tmp))
278 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
279 }
280 else
281 {
282 tmp = octave_sparse_params::get_key ("piv_tol");
283 if (! xisnan (tmp))
284 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
285
286 tmp = octave_sparse_params::get_key ("sym_tol");
287 if (! xisnan (tmp))
288 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
289 }
290
291 if (droptol >= 0.)
292 Control (UMFPACK_DROPTOL) = droptol;
293
294 // Set whether we are allowed to modify Q or not
295 if (FixedQ)
296 Control (UMFPACK_FIXQ) = 1.0;
297 else
298 {
299 tmp = octave_sparse_params::get_key ("autoamd");
300 if (! xisnan (tmp))
301 Control (UMFPACK_FIXQ) = tmp;
302 }
303
304 // Turn-off UMFPACK scaling for LU
305 if (scale)
306 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
307 else
308 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
309
310 UMFPACK_ZNAME (report_control) (control);
311
312 const octave_idx_type *Ap = a.cidx ();
313 const octave_idx_type *Ai = a.ridx ();
314 const Complex *Ax = a.data ();
315
316 UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
317 reinterpret_cast<const double *> (Ax), 0,
318 1, control);
319
320 void *Symbolic;
321 Matrix Info (1, UMFPACK_INFO);
322 double *info = Info.fortran_vec ();
323 int status;
324
325 // Null loop so that qinit is imediately deallocated when not
326 // needed
327 do
328 {
329 OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc);
330
331 for (octave_idx_type i = 0; i < nc; i++)
332 qinit[i] = static_cast<octave_idx_type> (Qinit (i));
333
334 status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai,
335 reinterpret_cast<const double *> (Ax),
336 0, qinit, &Symbolic, control,
337 info);
338 }
339 while (0);
340
341 if (status < 0)
342 {
343 UMFPACK_ZNAME (report_status) (control, status);
344 UMFPACK_ZNAME (report_info) (control, info);
345
346 UMFPACK_ZNAME (free_symbolic) (&Symbolic);
347
348 (*current_liboctave_error_handler)
349 ("SparseComplexLU::SparseComplexLU symbolic factorization failed");
350 }
351 else
352 {
353 UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
354
355 void *Numeric;
356 status = UMFPACK_ZNAME (numeric) (Ap, Ai,
357 reinterpret_cast<const double *> (Ax), 0,
358 Symbolic, &Numeric, control, info);
359 UMFPACK_ZNAME (free_symbolic) (&Symbolic);
360
361 cond = Info (UMFPACK_RCOND);
362
363 if (status < 0)
273 { 364 {
274 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0)); 365 UMFPACK_ZNAME (report_status) (control, status);
275 if (! xisnan (tmp)) 366 UMFPACK_ZNAME (report_info) (control, info);
276 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; 367
277 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1)); 368 UMFPACK_ZNAME (free_numeric) (&Numeric);
278 if (! xisnan (tmp)) 369
279 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; 370 (*current_liboctave_error_handler)
371 ("SparseComplexLU::SparseComplexLU numeric factorization failed");
280 } 372 }
281 else 373 else
282 { 374 {
283 tmp = octave_sparse_params::get_key ("piv_tol"); 375 UMFPACK_ZNAME (report_numeric) (Numeric, control);
284 if (! xisnan (tmp)) 376
285 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; 377 octave_idx_type lnz, unz, ignore1, ignore2, ignore3;
286 378 status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz,
287 tmp = octave_sparse_params::get_key ("sym_tol"); 379 &ignore1, &ignore2, &ignore3,
288 if (! xisnan (tmp)) 380 Numeric);
289 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
290 }
291
292 if (droptol >= 0.)
293 Control (UMFPACK_DROPTOL) = droptol;
294
295 // Set whether we are allowed to modify Q or not
296 if (FixedQ)
297 Control (UMFPACK_FIXQ) = 1.0;
298 else
299 {
300 tmp = octave_sparse_params::get_key ("autoamd");
301 if (! xisnan (tmp))
302 Control (UMFPACK_FIXQ) = tmp;
303 }
304
305 // Turn-off UMFPACK scaling for LU
306 if (scale)
307 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
308 else
309 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
310
311 UMFPACK_ZNAME (report_control) (control);
312
313 const octave_idx_type *Ap = a.cidx ();
314 const octave_idx_type *Ai = a.ridx ();
315 const Complex *Ax = a.data ();
316
317 UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
318 reinterpret_cast<const double *> (Ax), 0,
319 1, control);
320
321 void *Symbolic;
322 Matrix Info (1, UMFPACK_INFO);
323 double *info = Info.fortran_vec ();
324 int status;
325
326 // Null loop so that qinit is imediately deallocated when not
327 // needed
328 do
329 {
330 OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc);
331
332 for (octave_idx_type i = 0; i < nc; i++)
333 qinit[i] = static_cast<octave_idx_type> (Qinit (i));
334
335 status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai,
336 reinterpret_cast<const double *> (Ax),
337 0, qinit, &Symbolic, control,
338 info);
339 }
340 while (0);
341
342 if (status < 0)
343 {
344 (*current_liboctave_error_handler)
345 ("SparseComplexLU::SparseComplexLU symbolic factorization failed");
346
347 UMFPACK_ZNAME (report_status) (control, status);
348 UMFPACK_ZNAME (report_info) (control, info);
349
350 UMFPACK_ZNAME (free_symbolic) (&Symbolic);
351 }
352 else
353 {
354 UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
355
356 void *Numeric;
357 status = UMFPACK_ZNAME (numeric) (Ap, Ai,
358 reinterpret_cast<const double *> (Ax), 0,
359 Symbolic, &Numeric, control, info);
360 UMFPACK_ZNAME (free_symbolic) (&Symbolic);
361
362 cond = Info (UMFPACK_RCOND);
363 381
364 if (status < 0) 382 if (status < 0)
365 { 383 {
366 (*current_liboctave_error_handler)
367 ("SparseComplexLU::SparseComplexLU numeric factorization failed");
368
369 UMFPACK_ZNAME (report_status) (control, status); 384 UMFPACK_ZNAME (report_status) (control, status);
370 UMFPACK_ZNAME (report_info) (control, info); 385 UMFPACK_ZNAME (report_info) (control, info);
371 386
372 UMFPACK_ZNAME (free_numeric) (&Numeric); 387 UMFPACK_ZNAME (free_numeric) (&Numeric);
388
389 (*current_liboctave_error_handler)
390 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
373 } 391 }
374 else 392 else
375 { 393 {
376 UMFPACK_ZNAME (report_numeric) (Numeric, control); 394 octave_idx_type n_inner = (nr < nc ? nr : nc);
377 395
378 octave_idx_type lnz, unz, ignore1, ignore2, ignore3; 396 if (lnz < 1)
379 status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz, 397 Lfact = SparseComplexMatrix (n_inner, nr,
380 &ignore1, &ignore2, &ignore3, 398 static_cast<octave_idx_type> (1));
381 Numeric); 399 else
400 Lfact = SparseComplexMatrix (n_inner, nr, lnz);
401
402 octave_idx_type *Ltp = Lfact.cidx ();
403 octave_idx_type *Ltj = Lfact.ridx ();
404 Complex *Ltx = Lfact.data ();
405
406 if (unz < 1)
407 Ufact = SparseComplexMatrix (n_inner, nc,
408 static_cast<octave_idx_type> (1));
409 else
410 Ufact = SparseComplexMatrix (n_inner, nc, unz);
411
412 octave_idx_type *Up = Ufact.cidx ();
413 octave_idx_type *Uj = Ufact.ridx ();
414 Complex *Ux = Ufact.data ();
415
416 Rfact = SparseMatrix (nr, nr, nr);
417 for (octave_idx_type i = 0; i < nr; i++)
418 {
419 Rfact.xridx (i) = i;
420 Rfact.xcidx (i) = i;
421 }
422 Rfact.xcidx (nr) = nr;
423 double *Rx = Rfact.data ();
424
425 P.resize (dim_vector (nr, 1));
426 octave_idx_type *p = P.fortran_vec ();
427
428 Q.resize (dim_vector (nc, 1));
429 octave_idx_type *q = Q.fortran_vec ();
430
431 octave_idx_type do_recip;
432 status =
433 UMFPACK_ZNAME (get_numeric) (Ltp, Ltj,
434 reinterpret_cast<double *> (Ltx),
435 0, Up, Uj,
436 reinterpret_cast<double *> (Ux),
437 0, p, q, 0, 0,
438 &do_recip, Rx, Numeric);
439
440 UMFPACK_ZNAME (free_numeric) (&Numeric);
382 441
383 if (status < 0) 442 if (status < 0)
384 { 443 {
444 UMFPACK_ZNAME (report_status) (control, status);
445
385 (*current_liboctave_error_handler) 446 (*current_liboctave_error_handler)
386 ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); 447 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
387 448 }
388 UMFPACK_ZNAME (report_status) (control, status); 449 else
389 UMFPACK_ZNAME (report_info) (control, info); 450 {
390 451 Lfact = Lfact.transpose ();
391 UMFPACK_ZNAME (free_numeric) (&Numeric); 452
392 } 453 if (do_recip)
393 else 454 for (octave_idx_type i = 0; i < nr; i++)
394 { 455 Rx[i] = 1.0 / Rx[i];
395 octave_idx_type n_inner = (nr < nc ? nr : nc); 456
396 457 UMFPACK_ZNAME (report_matrix) (nr, n_inner,
397 if (lnz < 1) 458 Lfact.cidx (),
398 Lfact = SparseComplexMatrix (n_inner, nr, 459 Lfact.ridx (),
399 static_cast<octave_idx_type> (1)); 460 reinterpret_cast<double *> (Lfact.data ()),
400 else 461 0, 1, control);
401 Lfact = SparseComplexMatrix (n_inner, nr, lnz); 462
402 463 UMFPACK_ZNAME (report_matrix) (n_inner, nc,
403 octave_idx_type *Ltp = Lfact.cidx (); 464 Ufact.cidx (),
404 octave_idx_type *Ltj = Lfact.ridx (); 465 Ufact.ridx (),
405 Complex *Ltx = Lfact.data (); 466 reinterpret_cast<double *> (Ufact.data ()),
406 467 0, 1, control);
407 if (unz < 1) 468 UMFPACK_ZNAME (report_perm) (nr, p, control);
408 Ufact = SparseComplexMatrix (n_inner, nc, 469 UMFPACK_ZNAME (report_perm) (nc, q, control);
409 static_cast<octave_idx_type> (1)); 470 }
410 else 471
411 Ufact = SparseComplexMatrix (n_inner, nc, unz); 472 UMFPACK_ZNAME (report_info) (control, info);
412
413 octave_idx_type *Up = Ufact.cidx ();
414 octave_idx_type *Uj = Ufact.ridx ();
415 Complex *Ux = Ufact.data ();
416
417 Rfact = SparseMatrix (nr, nr, nr);
418 for (octave_idx_type i = 0; i < nr; i++)
419 {
420 Rfact.xridx (i) = i;
421 Rfact.xcidx (i) = i;
422 }
423 Rfact.xcidx (nr) = nr;
424 double *Rx = Rfact.data ();
425
426 P.resize (dim_vector (nr, 1));
427 octave_idx_type *p = P.fortran_vec ();
428
429 Q.resize (dim_vector (nc, 1));
430 octave_idx_type *q = Q.fortran_vec ();
431
432 octave_idx_type do_recip;
433 status =
434 UMFPACK_ZNAME (get_numeric) (Ltp, Ltj,
435 reinterpret_cast<double *> (Ltx),
436 0, Up, Uj,
437 reinterpret_cast<double *> (Ux),
438 0, p, q, 0, 0,
439 &do_recip, Rx, Numeric);
440
441 UMFPACK_ZNAME (free_numeric) (&Numeric);
442
443 if (status < 0)
444 {
445 (*current_liboctave_error_handler)
446 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
447
448 UMFPACK_ZNAME (report_status) (control, status);
449 }
450 else
451 {
452 Lfact = Lfact.transpose ();
453
454 if (do_recip)
455 for (octave_idx_type i = 0; i < nr; i++)
456 Rx[i] = 1.0 / Rx[i];
457
458 UMFPACK_ZNAME (report_matrix) (nr, n_inner,
459 Lfact.cidx (),
460 Lfact.ridx (),
461 reinterpret_cast<double *> (Lfact.data ()),
462 0, 1, control);
463
464 UMFPACK_ZNAME (report_matrix) (n_inner, nc,
465 Ufact.cidx (),
466 Ufact.ridx (),
467 reinterpret_cast<double *> (Ufact.data ()),
468 0, 1, control);
469 UMFPACK_ZNAME (report_perm) (nr, p, control);
470 UMFPACK_ZNAME (report_perm) (nc, q, control);
471 }
472
473 UMFPACK_ZNAME (report_info) (control, info);
474 }
475 } 473 }
476 } 474 }
477 475 }
478 if (udiag) 476
479 (*current_liboctave_error_handler) 477 if (udiag)
480 ("Option udiag of incomplete LU not implemented"); 478 (*current_liboctave_error_handler)
481 } 479 ("Option udiag of incomplete LU not implemented");
482 480
483 #else 481 #else
484 (*current_liboctave_error_handler) 482 (*current_liboctave_error_handler)
485 ("support for UMFPACK was unavailable or disabled when liboctave was built"); 483 ("support for UMFPACK was unavailable or disabled when liboctave was built");
486 #endif 484 #endif