Mercurial > octave-nkf

comparison liboctave/SparsedbleLU.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
215 (*current_liboctave_error_handler) ("UMFPACK not installed"); 215 (*current_liboctave_error_handler) ("UMFPACK not installed");
216 #endif 216 #endif
217 } 217 }
218 218
219 SparseLU::SparseLU (const SparseMatrix& a, const ColumnVector& Qinit, 219 SparseLU::SparseLU (const SparseMatrix& a, const ColumnVector& Qinit,
220 double piv_thres, bool FixedQ) 220 double piv_thres, bool FixedQ, double droptol,
221 bool milu, bool udiag)
221 { 222 {
222 #ifdef HAVE_UMFPACK 223 #ifdef HAVE_UMFPACK
223 octave_idx_type nr = a.rows (); 224 if (milu)
224 octave_idx_type nc = a.cols (); 225 (*current_liboctave_error_handler)
225 226 ("Modified incomplete LU not implemented");
226 // Setup the control parameters
227 Matrix Control (UMFPACK_CONTROL, 1);
228 double *control = Control.fortran_vec ();
229 umfpack_di_defaults (control);
230
231 double tmp = Voctave_sparse_controls.get_key ("spumoni");
232 if (!xisnan (tmp))
233 Control (UMFPACK_PRL) = tmp;
234 if (piv_thres >= 0.)
235 {
236 piv_thres = (piv_thres > 1. ? 1. : piv_thres);
237 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres;
238 Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres;
239 }
240 else 227 else
241 { 228 {
242 tmp = Voctave_sparse_controls.get_key ("piv_tol"); 229 octave_idx_type nr = a.rows ();
230 octave_idx_type nc = a.cols ();
231
232 // Setup the control parameters
233 Matrix Control (UMFPACK_CONTROL, 1);
234 double *control = Control.fortran_vec ();
235 umfpack_di_defaults (control);
236
237 double tmp = Voctave_sparse_controls.get_key ("spumoni");
243 if (!xisnan (tmp)) 238 if (!xisnan (tmp))
244 { 239 Control (UMFPACK_PRL) = tmp;
245 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; 240 if (piv_thres >= 0.)
246 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; 241 {
247 } 242 piv_thres = (piv_thres > 1. ? 1. : piv_thres);
248 } 243 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres;
249 244 Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres;
250 // Set whether we are allowed to modify Q or not 245 }
251 if (FixedQ) 246 else
252 Control (UMFPACK_FIXQ) = 1.0; 247 {
253 else 248 tmp = Voctave_sparse_controls.get_key ("piv_tol");
254 { 249 if (!xisnan (tmp))
255 tmp = Voctave_sparse_controls.get_key ("autoamd"); 250 {
256 if (!xisnan (tmp)) 251 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
257 Control (UMFPACK_FIXQ) = tmp; 252 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
258 } 253 }
259 254 }
260 // Turn-off UMFPACK scaling for LU 255
261 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; 256 if (droptol >= 0.)
262 257 Control (UMFPACK_DROPTOL) = droptol;
263 umfpack_di_report_control (control); 258
264 259
265 const octave_idx_type *Ap = a.cidx (); 260 // Set whether we are allowed to modify Q or not
266 const octave_idx_type *Ai = a.ridx (); 261 if (FixedQ)
267 const double *Ax = a.data (); 262 Control (UMFPACK_FIXQ) = 1.0;
268 263 else
269 umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control); 264 {
270 265 tmp = Voctave_sparse_controls.get_key ("autoamd");
271 void *Symbolic; 266 if (!xisnan (tmp))
272 Matrix Info (1, UMFPACK_INFO); 267 Control (UMFPACK_FIXQ) = tmp;
273 double *info = Info.fortran_vec (); 268 }
274 int status; 269
275 270 // Turn-off UMFPACK scaling for LU
276 // Null loop so that qinit is imediately deallocated when not needed 271 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
277 do { 272
278 OCTAVE_LOCAL_BUFFER (int, qinit, nc); 273 umfpack_di_report_control (control);
279 274
280 for (int i = 0; i < nc; i++) 275 const octave_idx_type *Ap = a.cidx ();
281 qinit [i] = static_cast<int> (Qinit (i)); 276 const octave_idx_type *Ai = a.ridx ();
282 277 const double *Ax = a.data ();
283 status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, qinit, 278
284 &Symbolic, control, info); 279 umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control);
285 } while (0); 280
286 281 void *Symbolic;
287 if (status < 0) 282 Matrix Info (1, UMFPACK_INFO);
288 { 283 double *info = Info.fortran_vec ();
289 (*current_liboctave_error_handler) 284 int status;
285
286 // Null loop so that qinit is imediately deallocated when not needed
287 do {
288 OCTAVE_LOCAL_BUFFER (int, qinit, nc);
289
290 for (int i = 0; i < nc; i++)
291 qinit [i] = static_cast<int> (Qinit (i));
292
293 status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, qinit,
294 &Symbolic, control, info);
295 } while (0);
296
297 if (status < 0)
298 {
299 (*current_liboctave_error_handler)
290 ("SparseLU::SparseLU symbolic factorization failed"); 300 ("SparseLU::SparseLU symbolic factorization failed");
291
292 umfpack_di_report_status (control, status);
293 umfpack_di_report_info (control, info);
294
295 umfpack_di_free_symbolic (&Symbolic) ;
296 }
297 else
298 {
299 umfpack_di_report_symbolic (Symbolic, control);
300
301 void *Numeric;
302 status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric,
303 control, info) ;
304 umfpack_di_free_symbolic (&Symbolic) ;
305
306 cond = Info (UMFPACK_RCOND);
307
308 if (status < 0)
309 {
310 (*current_liboctave_error_handler)
311 ("SparseLU::SparseLU numeric factorization failed");
312 301
313 umfpack_di_report_status (control, status); 302 umfpack_di_report_status (control, status);
314 umfpack_di_report_info (control, info); 303 umfpack_di_report_info (control, info);
315 304
316 umfpack_di_free_numeric (&Numeric); 305 umfpack_di_free_symbolic (&Symbolic) ;
317 } 306 }
318 else 307 else
319 { 308 {
320 umfpack_di_report_numeric (Numeric, control); 309 umfpack_di_report_symbolic (Symbolic, control);
321 310
322 int lnz, unz, ignore1, ignore2, ignore3; 311 void *Numeric;
323 status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2, 312 status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric,
324 &ignore3, Numeric) ; 313 control, info) ;
325 314 umfpack_di_free_symbolic (&Symbolic) ;
315
316 cond = Info (UMFPACK_RCOND);
317
326 if (status < 0) 318 if (status < 0)
327 { 319 {
328 (*current_liboctave_error_handler) 320 (*current_liboctave_error_handler)
329 ("SparseLU::SparseLU extracting LU factors failed"); 321 ("SparseLU::SparseLU numeric factorization failed");
330 322
331 umfpack_di_report_status (control, status); 323 umfpack_di_report_status (control, status);
332 umfpack_di_report_info (control, info); 324 umfpack_di_report_info (control, info);
333 325
334 umfpack_di_free_numeric (&Numeric); 326 umfpack_di_free_numeric (&Numeric);
335 } 327 }
336 else 328 else
337 { 329 {
338 int n_inner = (nr < nc ? nr : nc); 330 umfpack_di_report_numeric (Numeric, control);
339 331
340 if (lnz < 1) 332 int lnz, unz, ignore1, ignore2, ignore3;
341 Lfact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nr, 333 status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2,
342 static_cast<octave_idx_type> (1)); 334 &ignore3, Numeric) ;
343 else 335
344 Lfact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nr, 336 if (status < 0)
345 static_cast<octave_idx_type> (lnz));
346
347 octave_idx_type *Ltp = Lfact.cidx ();
348 octave_idx_type *Ltj = Lfact.ridx ();
349 double *Ltx = Lfact.data ();
350
351 if (unz < 1)
352 Ufact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nc,
353 static_cast<octave_idx_type> (1));
354 else
355 Ufact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nc,
356 static_cast<octave_idx_type> (unz));
357
358 octave_idx_type *Up = Ufact.cidx ();
359 octave_idx_type *Uj = Ufact.ridx ();
360 double *Ux = Ufact.data ();
361
362 P.resize (nr);
363 int *p = P.fortran_vec ();
364
365 Q.resize (nc);
366 int *q = Q.fortran_vec ();
367
368 int do_recip;
369 status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Uj,
370 Ux, p, q, (double *) NULL,
371 &do_recip, (double *) NULL,
372 Numeric) ;
373
374 umfpack_di_free_numeric (&Numeric) ;
375
376 if (status < 0 || do_recip)
377 { 337 {
378 (*current_liboctave_error_handler) 338 (*current_liboctave_error_handler)
379 ("SparseLU::SparseLU extracting LU factors failed"); 339 ("SparseLU::SparseLU extracting LU factors failed");
380 340
381 umfpack_di_report_status (control, status); 341 umfpack_di_report_status (control, status);
342 umfpack_di_report_info (control, info);
343
344 umfpack_di_free_numeric (&Numeric);
382 } 345 }
383 else 346 else
384 { 347 {
385 Lfact = Lfact.transpose (); 348 int n_inner = (nr < nc ? nr : nc);
386 umfpack_di_report_matrix (nr, n_inner, Lfact.cidx (), 349
387 Lfact.ridx (), Lfact.data (), 350 if (lnz < 1)
388 1, control); 351 Lfact = SparseMatrix
389 umfpack_di_report_matrix (n_inner, nc, Ufact.cidx (), 352 (static_cast<octave_idx_type> (n_inner), nr,
390 Ufact.ridx (), Ufact.data (), 353 static_cast<octave_idx_type> (1));
391 1, control); 354 else
392 umfpack_di_report_perm (nr, p, control); 355 Lfact = SparseMatrix
393 umfpack_di_report_perm (nc, q, control); 356 (static_cast<octave_idx_type> (n_inner), nr,
357 static_cast<octave_idx_type> (lnz));
358
359 octave_idx_type *Ltp = Lfact.cidx ();
360 octave_idx_type *Ltj = Lfact.ridx ();
361 double *Ltx = Lfact.data ();
362
363 if (unz < 1)
364 Ufact = SparseMatrix
365 (static_cast<octave_idx_type> (n_inner), nc,
366 static_cast<octave_idx_type> (1));
367 else
368 Ufact = SparseMatrix
369 (static_cast<octave_idx_type> (n_inner), nc,
370 static_cast<octave_idx_type> (unz));
371
372 octave_idx_type *Up = Ufact.cidx ();
373 octave_idx_type *Uj = Ufact.ridx ();
374 double *Ux = Ufact.data ();
375
376 P.resize (nr);
377 int *p = P.fortran_vec ();
378
379 Q.resize (nc);
380 int *q = Q.fortran_vec ();
381
382 int do_recip;
383 status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Uj,
384 Ux, p, q,
385 (double *) NULL,
386 &do_recip,
387 (double *) NULL,
388 Numeric) ;
389
390 umfpack_di_free_numeric (&Numeric) ;
391
392 if (status < 0 || do_recip)
393 {
394 (*current_liboctave_error_handler)
395 ("SparseLU::SparseLU extracting LU factors failed");
396
397 umfpack_di_report_status (control, status);
398 }
399 else
400 {
401 Lfact = Lfact.transpose ();
402 umfpack_di_report_matrix (nr, n_inner,
403 Lfact.cidx (),
404 Lfact.ridx (),
405 Lfact.data (),
406 1, control);
407 umfpack_di_report_matrix (n_inner, nc,
408 Ufact.cidx (),
409 Ufact.ridx (),
410 Ufact.data (),
411 1, control);
412 umfpack_di_report_perm (nr, p, control);
413 umfpack_di_report_perm (nc, q, control);
414 }
415
416 umfpack_di_report_info (control, info);
394 } 417 }
395
396 umfpack_di_report_info (control, info);
397 } 418 }
398 } 419 }
420
421 if (udiag)
422 (*current_liboctave_error_handler)
423 ("Option udiag of incomplete LU not implemented");
399 } 424 }
400 #else 425 #else
401 (*current_liboctave_error_handler) ("UMFPACK not installed"); 426 (*current_liboctave_error_handler) ("UMFPACK not installed");
402 #endif 427 #endif
403 } 428 }