Mercurial > hg > octave-nkf

comparison liboctave/SparseCmplxLU.cc @ 5282:5bdc3f24cd5f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-04-14 22:17:26 by dbateman]
author dbateman
date Thu, 14 Apr 2005 22:17:27 +0000
parents 23b37da9fd5b
children 4c8a2e4e0717
comparison
equal deleted inserted replaced
5281:f3266e7dbb99 5282:5bdc3f24cd5f
221 #endif 221 #endif
222 } 222 }
223 223
224 SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a, 224 SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a,
225 const ColumnVector& Qinit, 225 const ColumnVector& Qinit,
226 double piv_thres, bool FixedQ) 226 double piv_thres, bool FixedQ,
227 double droptol, bool milu, bool udiag)
227 { 228 {
228 #ifdef HAVE_UMFPACK 229 #ifdef HAVE_UMFPACK
229 octave_idx_type nr = a.rows (); 230 if (milu)
230 octave_idx_type nc = a.cols (); 231 (*current_liboctave_error_handler)
231 232 ("Modified incomplete LU not implemented");
232 // Setup the control parameters
233 Matrix Control (UMFPACK_CONTROL, 1);
234 double *control = Control.fortran_vec ();
235 umfpack_zi_defaults (control);
236
237 double tmp = Voctave_sparse_controls.get_key ("spumoni");
238 if (!xisnan (tmp))
239 Control (UMFPACK_PRL) = tmp;
240 if (piv_thres >= 0.)
241 {
242 piv_thres = (piv_thres > 1. ? 1. : piv_thres);
243 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres;
244 Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres;
245 }
246 else 233 else
247 { 234 {
248 tmp = Voctave_sparse_controls.get_key ("piv_tol"); 235 octave_idx_type nr = a.rows ();
236 octave_idx_type nc = a.cols ();
237
238 // Setup the control parameters
239 Matrix Control (UMFPACK_CONTROL, 1);
240 double *control = Control.fortran_vec ();
241 umfpack_zi_defaults (control);
242
243 double tmp = Voctave_sparse_controls.get_key ("spumoni");
249 if (!xisnan (tmp)) 244 if (!xisnan (tmp))
250 { 245 Control (UMFPACK_PRL) = tmp;
251 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; 246 if (piv_thres >= 0.)
252 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; 247 {
253 } 248 piv_thres = (piv_thres > 1. ? 1. : piv_thres);
254 } 249 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres;
255 250 Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres;
256 // Set whether we are allowed to modify Q or not 251 }
257 if (FixedQ) 252 else
258 Control (UMFPACK_FIXQ) = 1.0; 253 {
259 else 254 tmp = Voctave_sparse_controls.get_key ("piv_tol");
260 { 255 if (!xisnan (tmp))
261 tmp = Voctave_sparse_controls.get_key ("autoamd"); 256 {
262 if (!xisnan (tmp)) 257 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
263 Control (UMFPACK_FIXQ) = tmp; 258 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
264 } 259 }
265 260 }
266 // Turn-off UMFPACK scaling for LU 261
267 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; 262 if (droptol >= 0.)
268 263 Control (UMFPACK_DROPTOL) = droptol;
269 umfpack_zi_report_control (control); 264
270 265 // Set whether we are allowed to modify Q or not
271 const octave_idx_type *Ap = a.cidx (); 266 if (FixedQ)
272 const octave_idx_type *Ai = a.ridx (); 267 Control (UMFPACK_FIXQ) = 1.0;
273 const Complex *Ax = a.data (); 268 else
274 269 {
275 umfpack_zi_report_matrix (nr, nc, Ap, Ai, X_CAST (const double *, Ax), NULL, 270 tmp = Voctave_sparse_controls.get_key ("autoamd");
276 1, control); 271 if (!xisnan (tmp))
277 272 Control (UMFPACK_FIXQ) = tmp;
278 void *Symbolic; 273 }
279 Matrix Info (1, UMFPACK_INFO); 274
280 double *info = Info.fortran_vec (); 275 // Turn-off UMFPACK scaling for LU
281 int status; 276 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
282 277
283 // Null loop so that qinit is imediately deallocated when not needed 278 umfpack_zi_report_control (control);
284 do { 279
285 OCTAVE_LOCAL_BUFFER (int, qinit, nc); 280 const octave_idx_type *Ap = a.cidx ();
286 281 const octave_idx_type *Ai = a.ridx ();
287 for (int i = 0; i < nc; i++) 282 const Complex *Ax = a.data ();
288 qinit [i] = static_cast<int> (Qinit (i)); 283
289 284 umfpack_zi_report_matrix (nr, nc, Ap, Ai,
290 status = umfpack_zi_qsymbolic (nr, nc, Ap, Ai, X_CAST (const double *, Ax), 285 X_CAST (const double *, Ax), NULL,
291 NULL, qinit, &Symbolic, control, info); 286 1, control);
292 } while (0); 287
293 288 void *Symbolic;
294 if (status < 0) 289 Matrix Info (1, UMFPACK_INFO);
295 { 290 double *info = Info.fortran_vec ();
296 (*current_liboctave_error_handler) 291 int status;
292
293 // Null loop so that qinit is imediately deallocated when not
294 // needed
295 do {
296 OCTAVE_LOCAL_BUFFER (int, qinit, nc);
297
298 for (int i = 0; i < nc; i++)
299 qinit [i] = static_cast<int> (Qinit (i));
300
301 status = umfpack_zi_qsymbolic (nr, nc, Ap, Ai,
302 X_CAST (const double *, Ax),
303 NULL, qinit, &Symbolic, control,
304 info);
305 } while (0);
306
307 if (status < 0)
308 {
309 (*current_liboctave_error_handler)
297 ("SparseComplexLU::SparseComplexLU symbolic factorization failed"); 310 ("SparseComplexLU::SparseComplexLU symbolic factorization failed");
298
299 umfpack_zi_report_status (control, status);
300 umfpack_zi_report_info (control, info);
301
302 umfpack_zi_free_symbolic (&Symbolic) ;
303 }
304 else
305 {
306 umfpack_zi_report_symbolic (Symbolic, control);
307
308 void *Numeric;
309 status = umfpack_zi_numeric (Ap, Ai, X_CAST (const double *, Ax), NULL,
310 Symbolic, &Numeric, control, info) ;
311 umfpack_zi_free_symbolic (&Symbolic) ;
312
313 cond = Info (UMFPACK_RCOND);
314
315 if (status < 0)
316 {
317 (*current_liboctave_error_handler)
318 ("SparseComplexLU::SparseComplexLU numeric factorization failed");
319 311
320 umfpack_zi_report_status (control, status); 312 umfpack_zi_report_status (control, status);
321 umfpack_zi_report_info (control, info); 313 umfpack_zi_report_info (control, info);
322 314
323 umfpack_zi_free_numeric (&Numeric); 315 umfpack_zi_free_symbolic (&Symbolic) ;
324 } 316 }
325 else 317 else
326 { 318 {
327 umfpack_zi_report_numeric (Numeric, control); 319 umfpack_zi_report_symbolic (Symbolic, control);
328 320
329 int lnz, unz, ignore1, ignore2, ignore3; 321 void *Numeric;
330 status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1, &ignore2, 322 status = umfpack_zi_numeric (Ap, Ai,
331 &ignore3, Numeric) ; 323 X_CAST (const double *, Ax), NULL,
332 324 Symbolic, &Numeric, control, info) ;
325 umfpack_zi_free_symbolic (&Symbolic) ;
326
327 cond = Info (UMFPACK_RCOND);
328
333 if (status < 0) 329 if (status < 0)
334 { 330 {
335 (*current_liboctave_error_handler) 331 (*current_liboctave_error_handler)
336 ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); 332 ("SparseComplexLU::SparseComplexLU numeric factorization failed");
337 333
338 umfpack_zi_report_status (control, status); 334 umfpack_zi_report_status (control, status);
339 umfpack_zi_report_info (control, info); 335 umfpack_zi_report_info (control, info);
340 336
341 umfpack_zi_free_numeric (&Numeric); 337 umfpack_zi_free_numeric (&Numeric);
342 } 338 }
343 else 339 else
344 { 340 {
345 int n_inner = (nr < nc ? nr : nc); 341 umfpack_zi_report_numeric (Numeric, control);
346 342
347 if (lnz < 1) 343 int lnz, unz, ignore1, ignore2, ignore3;
348 Lfact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nr, 344 status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1,
349 static_cast<octave_idx_type> (1)); 345 &ignore2, &ignore3, Numeric);
350 else 346
351 Lfact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nr, 347 if (status < 0)
352 static_cast<octave_idx_type> (lnz));
353
354 octave_idx_type *Ltp = Lfact.cidx ();
355 octave_idx_type *Ltj = Lfact.ridx ();
356 Complex *Ltx = Lfact.data ();
357
358 if (unz < 1)
359 Ufact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nc,
360 static_cast<octave_idx_type> (1));
361 else
362 Ufact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nc, unz);
363
364 octave_idx_type *Up = Ufact.cidx ();
365 octave_idx_type *Uj = Ufact.ridx ();
366 Complex *Ux = Ufact.data ();
367
368 P.resize (nr);
369 int *p = P.fortran_vec ();
370
371 Q.resize (nc);
372 int *q = Q.fortran_vec ();
373
374 int do_recip;
375 status = umfpack_zi_get_numeric (Ltp, Ltj, X_CAST (double *, Ltx),
376 NULL, Up, Uj,
377 X_CAST (double *, Ux), NULL, p,
378 q, NULL, NULL, &do_recip,
379 NULL, Numeric) ;
380
381 umfpack_zi_free_numeric (&Numeric) ;
382
383 if (status < 0 || do_recip)
384 { 348 {
385 (*current_liboctave_error_handler) 349 (*current_liboctave_error_handler)
386 ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); 350 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
387 351
388 umfpack_zi_report_status (control, status); 352 umfpack_zi_report_status (control, status);
353 umfpack_zi_report_info (control, info);
354
355 umfpack_zi_free_numeric (&Numeric);
389 } 356 }
390 else 357 else
391 { 358 {
392 Lfact = Lfact.transpose (); 359 int n_inner = (nr < nc ? nr : nc);
393 360
394 umfpack_zi_report_matrix (nr, n_inner, Lfact.cidx (), 361 if (lnz < 1)
395 Lfact.ridx (), 362 Lfact = SparseComplexMatrix
396 X_CAST (double *, Lfact.data()), 363 (static_cast<octave_idx_type> (n_inner), nr,
397 NULL, 1, control); 364 static_cast<octave_idx_type> (1));
398 365 else
399 umfpack_zi_report_matrix (n_inner, nc, Ufact.cidx (), 366 Lfact = SparseComplexMatrix
400 Ufact.ridx (), 367 (static_cast<octave_idx_type> (n_inner), nr,
401 X_CAST (double *, Ufact.data()), 368 static_cast<octave_idx_type> (lnz));
402 NULL, 1, control); 369
403 umfpack_zi_report_perm (nr, p, control); 370 octave_idx_type *Ltp = Lfact.cidx ();
404 umfpack_zi_report_perm (nc, q, control); 371 octave_idx_type *Ltj = Lfact.ridx ();
372 Complex *Ltx = Lfact.data ();
373
374 if (unz < 1)
375 Ufact = SparseComplexMatrix
376 (static_cast<octave_idx_type> (n_inner), nc,
377 static_cast<octave_idx_type> (1));
378 else
379 Ufact = SparseComplexMatrix
380 (static_cast<octave_idx_type> (n_inner), nc, unz);
381
382 octave_idx_type *Up = Ufact.cidx ();
383 octave_idx_type *Uj = Ufact.ridx ();
384 Complex *Ux = Ufact.data ();
385
386 P.resize (nr);
387 int *p = P.fortran_vec ();
388
389 Q.resize (nc);
390 int *q = Q.fortran_vec ();
391
392 int do_recip;
393 status =
394 umfpack_zi_get_numeric (Ltp, Ltj,
395 X_CAST (double *, Ltx),
396 NULL, Up, Uj,
397 X_CAST (double *, Ux),
398 NULL, p, q, NULL, NULL,
399 &do_recip, NULL, Numeric) ;
400
401 umfpack_zi_free_numeric (&Numeric) ;
402
403 if (status < 0 || do_recip)
404 {
405 (*current_liboctave_error_handler)
406 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
407
408 umfpack_zi_report_status (control, status);
409 }
410 else
411 {
412 Lfact = Lfact.transpose ();
413
414 umfpack_zi_report_matrix (nr, n_inner,
415 Lfact.cidx (),
416 Lfact.ridx (),
417 X_CAST (double *, Lfact.data()),
418 NULL, 1, control);
419
420 umfpack_zi_report_matrix (n_inner, nc,
421 Ufact.cidx (),
422 Ufact.ridx (),
423 X_CAST (double *, Ufact.data()),
424 NULL, 1, control);
425 umfpack_zi_report_perm (nr, p, control);
426 umfpack_zi_report_perm (nc, q, control);
427 }
428
429 umfpack_zi_report_info (control, info);
405 } 430 }
406
407 umfpack_zi_report_info (control, info);
408 } 431 }
409 } 432 }
433
434 if (udiag)
435 (*current_liboctave_error_handler)
436 ("Option udiag of incomplete LU not implemented");
410 } 437 }
411 #else 438 #else
412 (*current_liboctave_error_handler) ("UMFPACK not installed"); 439 (*current_liboctave_error_handler) ("UMFPACK not installed");
413 #endif 440 #endif
414 } 441 }