25 #if defined (HAVE_CONFIG_H)
73 double *Rs,
void *Numeric);
79 void *Symbolic,
void **Numeric,
80 const double *Control,
double *Info);
88 const double *Control,
double *Info);
106 template <
typename T>
110 template <
typename T>
113 const double *Control);
115 template <
typename T>
119 template <
typename T>
123 #if defined (HAVE_UMFPACK)
155 return UMFPACK_DNAME (get_lunz) (lnz, unz, &ignore1, &ignore2,
167 return UMFPACK_DNAME (get_numeric) (Lp, Lj, Lx, Up, Ui, Ux,
p, q, Dx,
168 do_recip, Rs, Numeric);
175 const double *Ax,
void *Symbolic,
void **Numeric,
176 const double *Control,
double *Info)
178 return UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic, Numeric, Control,
188 const double *Control,
double *Info)
190 return UMFPACK_DNAME (qsymbolic) (n_row, n_col, Ap, Ai, Ax, Qinit,
191 Symbolic, Control, Info);
213 const double *Control)
278 return UMFPACK_ZNAME (get_lunz) (lnz, unz, &ignore1, &ignore2,
291 reinterpret_cast<double *
> (Lz),
293 reinterpret_cast<double *> (Uz),
295 reinterpret_cast<double *
> (Dz),
296 0, do_recip, Rs, Numeric);
303 const Complex *Az,
void *Symbolic,
void **Numeric,
304 const double *Control,
double *Info)
307 reinterpret_cast<const double *
> (Az),
308 0, Symbolic, Numeric, Control, Info);
317 void **Symbolic,
const double *Control,
double *Info)
320 reinterpret_cast<const double *
> (Az),
321 0, Qinit, Symbolic, Control, Info);
346 reinterpret_cast<const double *
> (Az),
347 0, col_form, Control);
381 template <
typename lu_type>
384 : Lfact (), Ufact (), Rfact (), cond (0),
P (),
Q ()
386 #if defined (HAVE_UMFPACK)
391 Matrix Control (UMFPACK_CONTROL, 1);
393 umfpack_defaults<lu_elt_type> (control);
397 Control (UMFPACK_PRL) =
tmp;
399 if (piv_thres.
numel () == 2)
401 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
403 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
405 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
407 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
413 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
417 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
423 Control (UMFPACK_FIXQ) =
tmp;
427 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
429 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
431 umfpack_report_control<lu_elt_type> (control);
437 umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax,
442 Matrix Info (1, UMFPACK_INFO);
444 int status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax, 0,
445 &Symbolic, control, info);
449 umfpack_report_status<lu_elt_type> (control, status);
450 umfpack_report_info<lu_elt_type> (control, info);
452 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
454 (*current_liboctave_error_handler)
455 (
"sparse_lu: symbolic factorization failed");
459 umfpack_report_symbolic<lu_elt_type> (Symbolic, control);
462 status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic,
463 &Numeric, control, info);
464 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
466 cond = Info (UMFPACK_RCOND);
470 umfpack_report_status<lu_elt_type> (control, status);
471 umfpack_report_info<lu_elt_type> (control, info);
473 umfpack_free_numeric<lu_elt_type> (&Numeric);
475 (*current_liboctave_error_handler)
476 (
"sparse_lu: numeric factorization failed");
480 umfpack_report_numeric<lu_elt_type> (Numeric, control);
483 status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);
487 umfpack_report_status<lu_elt_type> (control, status);
488 umfpack_report_info<lu_elt_type> (control, info);
490 umfpack_free_numeric<lu_elt_type> (&Numeric);
492 (*current_liboctave_error_handler)
493 (
"sparse_lu: extracting LU factors failed");
500 Lfact = lu_type (n_inner, nr,
501 static_cast<octave_idx_type> (1));
503 Lfact = lu_type (n_inner, nr, lnz);
510 Ufact = lu_type (n_inner, nc,
511 static_cast<octave_idx_type> (1));
513 Ufact = lu_type (n_inner, nc, unz);
535 status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx,
541 umfpack_free_numeric<lu_elt_type> (&Numeric);
545 umfpack_report_status<lu_elt_type> (control, status);
547 (*current_liboctave_error_handler)
548 (
"sparse_lu: extracting LU factors failed");
558 umfpack_report_matrix<lu_elt_type> (nr, n_inner,
564 umfpack_report_matrix<lu_elt_type> (n_inner, nc,
570 umfpack_report_perm<lu_elt_type> (nr,
p, control);
571 umfpack_report_perm<lu_elt_type> (nc, q, control);
574 umfpack_report_info<lu_elt_type> (control, info);
581 octave_unused_parameter (a);
582 octave_unused_parameter (piv_thres);
583 octave_unused_parameter (scale);
585 (*current_liboctave_error_handler)
586 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
591 template <
typename lu_type>
595 bool FixedQ,
double droptol,
596 bool milu,
bool udiag)
597 : Lfact (), Ufact (), Rfact (), cond (0),
P (),
Q ()
599 #if defined (HAVE_UMFPACK)
602 (*current_liboctave_error_handler)
603 (
"Modified incomplete LU not implemented");
609 Matrix Control (UMFPACK_CONTROL, 1);
611 umfpack_defaults<lu_elt_type> (control);
615 Control (UMFPACK_PRL) =
tmp;
617 if (piv_thres.
numel () == 2)
619 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
621 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
622 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
624 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
630 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
634 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
638 Control (UMFPACK_DROPTOL) = droptol;
642 Control (UMFPACK_FIXQ) = 1.0;
647 Control (UMFPACK_FIXQ) =
tmp;
652 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
654 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
656 umfpack_report_control<lu_elt_type> (control);
662 umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax,
667 Matrix Info (1, UMFPACK_INFO);
677 qinit[
i] = static_cast<octave_idx_type> (Qinit (
i));
679 status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax,
680 qinit, &Symbolic, control,
687 umfpack_report_status<lu_elt_type> (control, status);
688 umfpack_report_info<lu_elt_type> (control, info);
690 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
692 (*current_liboctave_error_handler)
693 (
"sparse_lu: symbolic factorization failed");
697 umfpack_report_symbolic<lu_elt_type> (Symbolic, control);
700 status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic,
701 &Numeric, control, info);
702 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
704 cond = Info (UMFPACK_RCOND);
708 umfpack_report_status<lu_elt_type> (control, status);
709 umfpack_report_info<lu_elt_type> (control, info);
711 umfpack_free_numeric<lu_elt_type> (&Numeric);
713 (*current_liboctave_error_handler)
714 (
"sparse_lu: numeric factorization failed");
718 umfpack_report_numeric<lu_elt_type> (Numeric, control);
721 status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);
725 umfpack_report_status<lu_elt_type> (control, status);
726 umfpack_report_info<lu_elt_type> (control, info);
728 umfpack_free_numeric<lu_elt_type> (&Numeric);
730 (*current_liboctave_error_handler)
731 (
"sparse_lu: extracting LU factors failed");
738 Lfact = lu_type (n_inner, nr,
739 static_cast<octave_idx_type> (1));
741 Lfact = lu_type (n_inner, nr, lnz);
748 Ufact = lu_type (n_inner, nc,
749 static_cast<octave_idx_type> (1));
751 Ufact = lu_type (n_inner, nc, unz);
773 status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx,
779 umfpack_free_numeric<lu_elt_type> (&Numeric);
783 umfpack_report_status<lu_elt_type> (control, status);
785 (*current_liboctave_error_handler)
786 (
"sparse_lu: extracting LU factors failed");
796 umfpack_report_matrix<lu_elt_type> (nr, n_inner,
802 umfpack_report_matrix<lu_elt_type> (n_inner, nc,
808 umfpack_report_perm<lu_elt_type> (nr,
p, control);
809 umfpack_report_perm<lu_elt_type> (nc, q, control);
812 umfpack_report_info<lu_elt_type> (control, info);
818 (*current_liboctave_error_handler)
819 (
"Option udiag of incomplete LU not implemented");
823 octave_unused_parameter (a);
824 octave_unused_parameter (Qinit);
825 octave_unused_parameter (piv_thres);
826 octave_unused_parameter (scale);
827 octave_unused_parameter (FixedQ);
828 octave_unused_parameter (droptol);
829 octave_unused_parameter (milu);
830 octave_unused_parameter (udiag);
832 (*current_liboctave_error_handler)
833 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
838 template <
typename lu_type>
846 lu_type Yout (nr, nc, Lfact.nnz () + Ufact.nnz () - (nr<nz?nr:nz));
854 Yout.xridx (ii) = Ufact.ridx (
i);
855 Yout.xdata (ii++) = Ufact.data (
i);
862 i < Lfact.cidx (j +1);
i++)
864 Yout.xridx (ii) = Lfact.ridx (
i);
865 Yout.xdata (ii++) = Lfact.data (
i);
869 Yout.xcidx (j + 1) = ii;
875 template <
typename lu_type>
895 template <
typename lu_type>
904 Pout.
xelem (
i) =
static_cast<double> (
P(
i) + 1);
909 template <
typename lu_type>
916 template <
typename lu_type>
936 template <
typename lu_type>
945 Pout.
xelem (
i) =
static_cast<double> (
Q(
i) + 1);
950 template <
typename lu_type>
void umfpack_report_perm< Complex >(octave_idx_type np, const octave_idx_type *Perm, const double *Control)
ColumnVector Pc_vec(void) const
octave_idx_type * xridx(void)
lu_type::element_type lu_elt_type
MArray< octave_idx_type > Q
void umfpack_defaults(double *Control)
octave_idx_type umfpack_numeric< double >(const octave_idx_type *Ap, const octave_idx_type *Ai, const double *Ax, void *Symbolic, void **Numeric, const double *Control, double *Info)
Octave interface to the compression and uncompression libraries.
void umfpack_report_status< double >(double *Control, octave_idx_type status)
SparseMatrix Pc(void) const
void umfpack_report_symbolic< double >(void *Symbolic, const double *Control)
octave_idx_type umfpack_qsymbolic(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const T *Ax, const octave_idx_type *Qinit, void **Symbolic, const double *Control, double *Info)
void umfpack_report_status< Complex >(double *Control, octave_idx_type status)
octave_idx_type numel(void) const
Number of elements in the array.
SparseMatrix Pr(void) const
PermMatrix Pr_mat(void) const
void umfpack_report_info< double >(const double *Control, const double *Info)
void umfpack_defaults< double >(double *Control)
void umfpack_report_numeric< Complex >(void *Numeric, const double *Control)
octave_idx_type * xcidx(void)
octave_idx_type umfpack_qsymbolic< double >(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const double *Ax, const octave_idx_type *Qinit, void **Symbolic, const double *Control, double *Info)
void umfpack_defaults< Complex >(double *Control)
octave_idx_type umfpack_numeric(const octave_idx_type *Ap, const octave_idx_type *Ai, const T *Ax, void *Symbolic, void **Numeric, const double *Control, double *Info)
#define UMFPACK_DNAME(name)
ColumnVector Pr_vec(void) const
octave_idx_type umfpack_get_numeric< double >(octave_idx_type *Lp, octave_idx_type *Lj, double *Lx, octave_idx_type *Up, octave_idx_type *Ui, double *Ux, octave_idx_type *p, octave_idx_type *q, double *Dx, octave_idx_type *do_recip, double *Rs, void *Numeric)
octave_idx_type * cidx(void)
#define UMFPACK_ZNAME(name)
static double get_key(const std::string &key)
void umfpack_report_matrix(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const T *Ax, octave_idx_type col_form, const double *Control)
void umfpack_free_numeric< Complex >(void **Numeric)
void umfpack_report_control(const double *Control)
void umfpack_free_symbolic< double >(void **Symbolic)
void umfpack_report_symbolic(void *Symbolic, const double *Control)
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
void umfpack_report_numeric(void *Numeric, const double *Control)
void umfpack_free_numeric(void **Numeric)
octave_idx_type umfpack_get_lunz< Complex >(octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
octave_idx_type umfpack_qsymbolic< Complex >(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const Complex *Az, const octave_idx_type *Qinit, void **Symbolic, const double *Control, double *Info)
void umfpack_free_symbolic< Complex >(void **Symbolic)
void umfpack_report_numeric< double >(void *Numeric, const double *Control)
void umfpack_free_symbolic(void **Symbolic)
void resize(const dim_vector &dv, const T &rfv)
octave_idx_type umfpack_get_numeric< Complex >(octave_idx_type *Lp, octave_idx_type *Lj, Complex *Lz, octave_idx_type *Up, octave_idx_type *Ui, Complex *Uz, octave_idx_type *p, octave_idx_type *q, double *Dz, octave_idx_type *do_recip, double *Rs, void *Numeric)
void umfpack_report_control< Complex >(const double *Control)
T & xelem(octave_idx_type n)
void umfpack_report_info(const double *Control, const double *Info)
octave_idx_type * ridx(void)
=val(i)}if ode{val(i)}occurs in table i
void umfpack_report_perm< double >(octave_idx_type np, const octave_idx_type *Perm, const double *Control)
void umfpack_report_matrix< Complex >(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const Complex *Az, octave_idx_type col_form, const double *Control)
PermMatrix Pc_mat(void) const
void scale(Matrix &m, double x, double y, double z)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
octave_idx_type umfpack_numeric< Complex >(const octave_idx_type *Ap, const octave_idx_type *Ai, const Complex *Az, void *Symbolic, void **Numeric, const double *Control, double *Info)
void umfpack_report_symbolic< Complex >(void *Symbolic, const double *Control)
void umfpack_report_perm(octave_idx_type np, const octave_idx_type *Perm, const double *Control)
MArray< octave_idx_type > P
std::complex< double > Complex
const T * fortran_vec(void) const
void umfpack_report_control< double >(const double *Control)
Vector representing the dimensions (size) of an Array.
void umfpack_free_numeric< double >(void **Numeric)
octave_idx_type umfpack_get_lunz< double >(octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
void umfpack_report_info< Complex >(const double *Control, const double *Info)
void umfpack_report_status(double *Control, octave_idx_type status)
void umfpack_report_matrix< double >(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const double *Ax, octave_idx_type col_form, const double *Control)
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Q
octave_idx_type umfpack_get_numeric(octave_idx_type *Lp, octave_idx_type *Lj, T *Lx, octave_idx_type *Up, octave_idx_type *Ui, T *Ux, octave_idx_type *p, octave_idx_type *q, double *Dz_packed, octave_idx_type *do_recip, double *Rs, void *Numeric)
octave_idx_type umfpack_get_lunz(octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)