25 #if defined (HAVE_CONFIG_H) 74 double *Rs,
void *Numeric);
80 void *Symbolic,
void **Numeric,
81 const double *Control,
double *Info);
89 const double *Control,
double *Info);
107 template <
typename T>
111 template <
typename T>
114 const double *Control);
116 template <
typename T>
120 template <
typename T>
124 #if defined (HAVE_UMFPACK) 158 &ignore1, &ignore2, &ignore3, Numeric);
183 const double *Ax,
void *Symbolic,
void **Numeric,
184 const double *Control,
double *Info)
188 Ax, Symbolic, Numeric, Control, Info);
197 const double *Control,
double *Info)
203 Symbolic, Control, Info);
225 const double *Control)
294 &ignore1, &ignore2, &ignore3, Numeric);
307 reinterpret_cast<double *
> (Lz),
310 reinterpret_cast<double *
> (Uz),
313 reinterpret_cast<double *
> (Dz),
322 const Complex *Az,
void *Symbolic,
void **Numeric,
323 const double *Control,
double *Info)
327 reinterpret_cast<const double *
> (Az),
328 nullptr, Symbolic, Numeric, Control, Info);
337 void **Symbolic,
const double *Control,
double *Info)
342 reinterpret_cast<const double *
> (Az),
344 Symbolic, Control, Info);
371 reinterpret_cast<const double *
> (Az),
372 nullptr, col_form, Control);
406 template <
typename lu_type>
409 : Lfact (), Ufact (), Rfact (), cond (0), P (),
Q ()
411 #if defined (HAVE_UMFPACK) 416 Matrix Control (UMFPACK_CONTROL, 1);
418 umfpack_defaults<lu_elt_type> (control);
422 Control (UMFPACK_PRL) =
tmp;
424 if (piv_thres.
numel () == 2)
426 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
428 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
430 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
432 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
438 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
442 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
448 Control (UMFPACK_FIXQ) =
tmp;
452 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
454 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
456 umfpack_report_control<lu_elt_type> (control);
462 umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax,
467 Matrix Info (1, UMFPACK_INFO);
469 int status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax,
nullptr,
470 &Symbolic, control, info);
474 umfpack_report_status<lu_elt_type> (control, status);
475 umfpack_report_info<lu_elt_type> (control, info);
477 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
479 (*current_liboctave_error_handler)
480 (
"sparse_lu: symbolic factorization failed");
484 umfpack_report_symbolic<lu_elt_type> (Symbolic, control);
487 status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic,
488 &Numeric, control, info);
489 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
491 cond = Info (UMFPACK_RCOND);
495 umfpack_report_status<lu_elt_type> (control, status);
496 umfpack_report_info<lu_elt_type> (control, info);
498 umfpack_free_numeric<lu_elt_type> (&Numeric);
500 (*current_liboctave_error_handler)
501 (
"sparse_lu: numeric factorization failed");
505 umfpack_report_numeric<lu_elt_type> (Numeric, control);
508 status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);
512 umfpack_report_status<lu_elt_type> (control, status);
513 umfpack_report_info<lu_elt_type> (control, info);
515 umfpack_free_numeric<lu_elt_type> (&Numeric);
517 (*current_liboctave_error_handler)
518 (
"sparse_lu: extracting LU factors failed");
525 Lfact = lu_type (n_inner, nr,
526 static_cast<octave_idx_type> (1));
528 Lfact = lu_type (n_inner, nr, lnz);
535 Ufact = lu_type (n_inner, nc,
536 static_cast<octave_idx_type> (1));
538 Ufact = lu_type (n_inner, nc, unz);
560 status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx,
566 umfpack_free_numeric<lu_elt_type> (&Numeric);
570 umfpack_report_status<lu_elt_type> (control, status);
572 (*current_liboctave_error_handler)
573 (
"sparse_lu: extracting LU factors failed");
583 umfpack_report_matrix<lu_elt_type> (nr, n_inner,
589 umfpack_report_matrix<lu_elt_type> (n_inner, nc,
595 umfpack_report_perm<lu_elt_type> (nr,
p, control);
596 umfpack_report_perm<lu_elt_type> (nc, q, control);
599 umfpack_report_info<lu_elt_type> (control, info);
606 octave_unused_parameter (
a);
607 octave_unused_parameter (piv_thres);
608 octave_unused_parameter (
scale);
610 (*current_liboctave_error_handler)
611 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
616 template <
typename lu_type>
620 bool FixedQ,
double droptol,
621 bool milu,
bool udiag)
622 : Lfact (), Ufact (), Rfact (), cond (0), P (),
Q ()
624 #if defined (HAVE_UMFPACK) 627 (*current_liboctave_error_handler)
628 (
"Modified incomplete LU not implemented");
634 Matrix Control (UMFPACK_CONTROL, 1);
636 umfpack_defaults<lu_elt_type> (control);
640 Control (UMFPACK_PRL) =
tmp;
642 if (piv_thres.
numel () == 2)
644 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
646 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
647 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
649 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
655 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
659 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
663 Control (UMFPACK_DROPTOL) = droptol;
667 Control (UMFPACK_FIXQ) = 1.0;
672 Control (UMFPACK_FIXQ) =
tmp;
677 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
679 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
681 umfpack_report_control<lu_elt_type> (control);
687 umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax,
692 Matrix Info (1, UMFPACK_INFO);
702 qinit[
i] = static_cast<octave_idx_type> (Qinit (
i));
704 status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax,
705 qinit, &Symbolic, control,
712 umfpack_report_status<lu_elt_type> (control, status);
713 umfpack_report_info<lu_elt_type> (control, info);
715 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
717 (*current_liboctave_error_handler)
718 (
"sparse_lu: symbolic factorization failed");
722 umfpack_report_symbolic<lu_elt_type> (Symbolic, control);
725 status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic,
726 &Numeric, control, info);
727 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
729 cond = Info (UMFPACK_RCOND);
733 umfpack_report_status<lu_elt_type> (control, status);
734 umfpack_report_info<lu_elt_type> (control, info);
736 umfpack_free_numeric<lu_elt_type> (&Numeric);
738 (*current_liboctave_error_handler)
739 (
"sparse_lu: numeric factorization failed");
743 umfpack_report_numeric<lu_elt_type> (Numeric, control);
746 status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);
750 umfpack_report_status<lu_elt_type> (control, status);
751 umfpack_report_info<lu_elt_type> (control, info);
753 umfpack_free_numeric<lu_elt_type> (&Numeric);
755 (*current_liboctave_error_handler)
756 (
"sparse_lu: extracting LU factors failed");
763 Lfact = lu_type (n_inner, nr,
764 static_cast<octave_idx_type> (1));
766 Lfact = lu_type (n_inner, nr, lnz);
773 Ufact = lu_type (n_inner, nc,
774 static_cast<octave_idx_type> (1));
776 Ufact = lu_type (n_inner, nc, unz);
798 status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx,
804 umfpack_free_numeric<lu_elt_type> (&Numeric);
808 umfpack_report_status<lu_elt_type> (control, status);
810 (*current_liboctave_error_handler)
811 (
"sparse_lu: extracting LU factors failed");
821 umfpack_report_matrix<lu_elt_type> (nr, n_inner,
827 umfpack_report_matrix<lu_elt_type> (n_inner, nc,
833 umfpack_report_perm<lu_elt_type> (nr,
p, control);
834 umfpack_report_perm<lu_elt_type> (nc, q, control);
837 umfpack_report_info<lu_elt_type> (control, info);
843 (*current_liboctave_error_handler)
844 (
"Option udiag of incomplete LU not implemented");
848 octave_unused_parameter (
a);
849 octave_unused_parameter (Qinit);
850 octave_unused_parameter (piv_thres);
851 octave_unused_parameter (
scale);
852 octave_unused_parameter (FixedQ);
853 octave_unused_parameter (droptol);
854 octave_unused_parameter (milu);
855 octave_unused_parameter (udiag);
857 (*current_liboctave_error_handler)
858 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
863 template <
typename lu_type>
871 lu_type Yout (nr, nc, Lfact.nnz () + Ufact.nnz () - (nr<nz?nr:nz));
879 Yout.xridx (ii) = Ufact.ridx (
i);
880 Yout.xdata (ii++) = Ufact.data (
i);
887 i < Lfact.cidx (j +1);
i++)
889 Yout.xridx (ii) = Lfact.ridx (
i);
890 Yout.xdata (ii++) = Lfact.data (
i);
894 Yout.xcidx (j + 1) = ii;
900 template <
typename lu_type>
920 template <
typename lu_type>
929 Pout.
xelem (
i) =
static_cast<double> (P(
i) + 1);
934 template <
typename lu_type>
941 template <
typename lu_type>
961 template <
typename lu_type>
970 Pout.
xelem (
i) =
static_cast<double> (
Q(
i) + 1);
975 template <
typename lu_type>
void umfpack_report_perm< Complex >(octave_idx_type np, const octave_idx_type *Perm, const double *Control)
lu_type::element_type lu_elt_type
MArray< octave_idx_type > Q
ColumnVector Pc_vec(void) const
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)
void umfpack_report_status< double >(double *Control, octave_idx_type status)
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)
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 * cidx(void)
const T * fortran_vec(void) const
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)
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)
SparseMatrix Pr(void) const
#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)
PermMatrix Pr_mat(void) const
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)
SparseMatrix Pc(void) const
void umfpack_free_symbolic< Complex >(void **Symbolic)
void umfpack_report_numeric< double >(void *Numeric, const double *Control)
void umfpack_free_symbolic(void **Symbolic)
PermMatrix Pc_mat(void) const
void resize(const dim_vector &dv, const T &rfv)
Resizing (with fill).
octave_idx_type * ridx(void)
octave_idx_type * xridx(void)
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)
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
void umfpack_report_info(const double *Control, const double *Info)
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)
void scale(Matrix &m, double x, double y, double z)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
ColumnVector Pr_vec(void) const
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
octave_idx_type * xcidx(void)
octave_idx_type numel(void) const
Number of elements in the array.
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)
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
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)