25 #if defined (HAVE_CONFIG_H) 60 #if ! defined (USE_QRSOLVE) 109 if (nr != nr_a || nc != nc_a || nz != nz_a)
126 return !(*
this ==
a);
135 if (nr == nc && nr > 0)
188 return max (dummy_idx, dim);
212 if (nr == 0 || nc == 0 || dim >=
dv.
ndims ())
222 if (
ridx (
i) != idx_j)
256 double tmp =
elem (idx_arg(j), j);
270 if (nr == 0 || nc == 0 || dim >=
dv.
ndims ())
297 else if (ix == -1 ||
tmp >
elem (ir, ix))
298 idx_arg.
elem (ir) = j;
304 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
314 if (idx_arg(j) == -1)
322 double tmp =
elem (j, idx_arg(j));
339 return min (dummy_idx, dim);
363 if (nr == 0 || nc == 0 || dim >=
dv.
ndims ())
373 if (
ridx (
i) != idx_j)
407 double tmp =
elem (idx_arg(j), j);
421 if (nr == 0 || nc == 0 || dim >=
dv.
ndims ())
448 else if (ix == -1 ||
tmp <
elem (ir, ix))
449 idx_arg.
elem (ir) = j;
455 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
465 if (idx_arg(j) == -1)
473 double tmp =
elem (j, idx_arg(j));
537 if (rb.
rows () > 0 && rb.
cols () > 0)
547 if (rb.
rows () > 0 && rb.
cols () > 0)
607 return inverse (mattype, info, rcond, 0, 0);
615 return inverse (mattype, info, rcond, 0, 0);
622 return inverse (mattype, info, rcond, 0, 0);
627 double& rcond,
const bool,
628 const bool calccond)
const 636 if (nr == 0 || nc == 0 || nr != nc)
637 (*current_liboctave_error_handler) (
"inverse requires square matrix");
640 int typ = mattype.
type ();
644 (*current_liboctave_error_handler) (
"incorrect matrix type");
652 double *v =
retval.data ();
660 double tmp = fabs (v[
i]);
677 double& rcond,
const bool,
678 const bool calccond)
const 686 if (nr == 0 || nc == 0 || nr != nc)
687 (*current_liboctave_error_handler) (
"inverse requires square matrix");
690 int typ = mattype.
type ();
695 (*current_liboctave_error_handler) (
"incorrect matrix type");
698 double ainvnorm = 0.;
707 atmp += fabs (
data (
i));
729 retval.change_capacity (nz2);
747 (*current_liboctave_error_handler) (
"division by zero");
752 rpX =
retval.xridx (colXp);
766 while (rpX < j && rpU < j && colXp < cx && colUp < nz);
770 colUp =
cidx (j+1) - 1;
773 double pivot =
data (colUp);
774 if (pivot == 0. ||
ridx (colUp) != j)
775 (*current_liboctave_error_handler) (
"division by zero");
782 retval.change_capacity (nz2);
786 retval.xdata (cx) = v / pivot;
794 colUp =
cidx (
i+1) - 1;
797 double pivot =
data (colUp);
798 if (pivot == 0. ||
ridx (colUp) !=
i)
799 (*current_liboctave_error_handler) (
"division by zero");
803 retval.xdata (j) /= pivot;
863 (*current_liboctave_error_handler) (
"division by zero");
871 colUp =
cidx (perm[iidx]+1) - 1;
873 colUp =
cidx (perm[iidx]);
875 double pivot =
data (colUp);
877 (*current_liboctave_error_handler) (
"division by zero");
890 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
891 retval.change_capacity (nz2);
899 retval.xdata (cx++) = work[j];
915 atmp += fabs (
retval.data (
i));
920 rcond = 1. / ainvnorm / anorm;
928 double& rcond,
bool,
bool calc_cond)
const 930 int typ = mattype.
type (
false);
934 typ = mattype.
type (*
this);
937 ret =
dinverse (mattype, info, rcond,
true, calc_cond);
951 rcond = fact.
rcond ();
957 info, rcond2,
true,
false);
958 ret =
Q * InvL.
transpose () * InvL *
Q.transpose ();
978 rcond = fact.
rcond ();
981 info, rcond2,
true,
false);
1011 #if defined (HAVE_UMFPACK) 1016 if (nr == 0 || nc == 0 || nr != nc)
1025 Matrix Control (UMFPACK_CONTROL, 1);
1031 Control (UMFPACK_PRL) =
tmp;
1036 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
1037 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
1043 Control (UMFPACK_FIXQ) =
tmp;
1046 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1052 const double *Ax =
data ();
1060 Matrix Info (1, UMFPACK_INFO);
1065 Ax,
nullptr, &Symbolic, control, info);
1074 (*current_liboctave_error_handler)
1075 (
"SparseMatrix::determinant symbolic factorization failed");
1085 &Numeric, control, info);
1088 rcond = Info (UMFPACK_RCOND);
1096 (*current_liboctave_error_handler)
1097 (
"SparseMatrix::determinant numeric factorization failed");
1105 status =
UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric,
1113 (*current_liboctave_error_handler)
1114 (
"SparseMatrix::determinant error calculating determinant");
1126 octave_unused_parameter (
err);
1127 octave_unused_parameter (rcond);
1129 (*current_liboctave_error_handler)
1130 (
"support for UMFPACK was unavailable or disabled " 1131 "when liboctave was built");
1141 double& rcond, solve_singularity_handler,
1142 bool calc_cond)
const 1151 if (nr !=
b.rows ())
1153 (
"matrix dimension mismatch solution of linear equations");
1155 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1160 int typ = mattype.
type ();
1164 (*current_liboctave_error_handler) (
"incorrect matrix type");
1189 rcond = dmin / dmax;
1201 solve_singularity_handler,
bool calc_cond)
const 1210 if (nr !=
b.rows ())
1212 (
"matrix dimension mismatch solution of linear equations");
1214 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1219 int typ = mattype.
type ();
1223 (*current_liboctave_error_handler) (
"incorrect matrix type");
1236 if (
b.ridx (
i) >=
nm)
1251 for (
k =
b.cidx (j);
k <
b.cidx (j+1);
k++)
1278 rcond = dmin / dmax;
1290 solve_singularity_handler,
bool calc_cond)
const 1299 if (nr !=
b.rows ())
1301 (
"matrix dimension mismatch solution of linear equations");
1303 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1308 int typ = mattype.
type ();
1312 (*current_liboctave_error_handler) (
"incorrect matrix type");
1337 rcond = dmin / dmax;
1349 solve_singularity_handler,
bool calc_cond)
const 1358 if (nr !=
b.rows ())
1360 (
"matrix dimension mismatch solution of linear equations");
1362 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1367 int typ = mattype.
type ();
1371 (*current_liboctave_error_handler) (
"incorrect matrix type");
1384 if (
b.ridx (
i) >=
nm)
1399 for (
k =
b.cidx (j);
k <
b.cidx (j+1);
k++)
1426 rcond = dmin / dmax;
1438 solve_singularity_handler sing_handler,
1439 bool calc_cond)
const 1448 if (nr !=
b.rows ())
1450 (
"matrix dimension mismatch solution of linear equations");
1452 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1457 int typ = mattype.
type ();
1461 (*current_liboctave_error_handler) (
"incorrect matrix type");
1464 double ainvnorm = 0.;
1475 atmp += fabs (
data (
i));
1504 goto triangular_error;
1510 i <
cidx (kidx+1)-1;
i++)
1513 work[iidx] = work[iidx] -
tmp *
data (
i);
1519 retval.xelem (perm[
i], j) = work[
i];
1541 i <
cidx (iidx+1)-1;
i++)
1544 work[idx2] = work[idx2] -
tmp *
data (
i);
1551 atmp += fabs (work[
i]);
1554 if (atmp > ainvnorm)
1557 rcond = 1. / ainvnorm / anorm;
1580 goto triangular_error;
1588 work[iidx] = work[iidx] -
tmp *
data (
i);
1616 work[iidx] = work[iidx] -
tmp *
data (
i);
1623 atmp += fabs (work[
i]);
1626 if (atmp > ainvnorm)
1629 rcond = 1. / ainvnorm / anorm;
1638 sing_handler (rcond);
1645 volatile double rcond_plus_one = rcond + 1.0;
1653 sing_handler (rcond);
1667 solve_singularity_handler sing_handler,
1668 bool calc_cond)
const 1677 if (nr !=
b.rows ())
1679 (
"matrix dimension mismatch solution of linear equations");
1681 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1686 int typ = mattype.
type ();
1690 (*current_liboctave_error_handler) (
"incorrect matrix type");
1693 double ainvnorm = 0.;
1703 atmp += fabs (
data (
i));
1730 work[
b.ridx (
i)] =
b.data (
i);
1742 goto triangular_error;
1748 i <
cidx (kidx+1)-1;
i++)
1751 work[iidx] = work[iidx] -
tmp *
data (
i);
1763 if (ii + new_nnz > x_nz)
1772 if (work[rperm[
i]] != 0.)
1775 retval.xdata (ii++) = work[rperm[
i]];
1780 retval.maybe_compress ();
1801 i <
cidx (iidx+1)-1;
i++)
1804 work[idx2] = work[idx2] -
tmp *
data (
i);
1811 atmp += fabs (work[
i]);
1814 if (atmp > ainvnorm)
1817 rcond = 1. / ainvnorm / anorm;
1829 work[
b.ridx (
i)] =
b.data (
i);
1839 goto triangular_error;
1847 work[iidx] = work[iidx] -
tmp *
data (
i);
1859 if (ii + new_nnz > x_nz)
1871 retval.xdata (ii++) = work[
i];
1876 retval.maybe_compress ();
1898 work[iidx] = work[iidx] -
tmp *
data (
i);
1905 atmp += fabs (work[
i]);
1908 if (atmp > ainvnorm)
1911 rcond = 1. / ainvnorm / anorm;
1920 sing_handler (rcond);
1927 volatile double rcond_plus_one = rcond + 1.0;
1935 sing_handler (rcond);
1948 solve_singularity_handler sing_handler,
1949 bool calc_cond)
const 1958 if (nr !=
b.rows ())
1960 (
"matrix dimension mismatch solution of linear equations");
1962 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1967 int typ = mattype.
type ();
1971 (*current_liboctave_error_handler) (
"incorrect matrix type");
1974 double ainvnorm = 0.;
1985 atmp += fabs (
data (
i));
2014 goto triangular_error;
2020 i <
cidx (kidx+1)-1;
i++)
2023 cwork[iidx] = cwork[iidx] -
tmp *
data (
i);
2029 retval.xelem (perm[
i], j) = cwork[
i];
2052 i <
cidx (iidx+1)-1;
i++)
2055 work[idx2] = work[idx2] -
tmp *
data (
i);
2062 atmp += fabs (work[
i]);
2065 if (atmp > ainvnorm)
2068 rcond = 1. / ainvnorm / anorm;
2091 goto triangular_error;
2099 cwork[iidx] = cwork[iidx] -
tmp *
data (
i);
2129 work[iidx] = work[iidx] -
tmp *
data (
i);
2136 atmp += fabs (work[
i]);
2139 if (atmp > ainvnorm)
2142 rcond = 1. / ainvnorm / anorm;
2151 sing_handler (rcond);
2158 volatile double rcond_plus_one = rcond + 1.0;
2166 sing_handler (rcond);
2180 solve_singularity_handler sing_handler,
2181 bool calc_cond)
const 2190 if (nr !=
b.rows ())
2192 (
"matrix dimension mismatch solution of linear equations");
2194 if (nr == 0 || nc == 0 ||
b.cols () == 0)
2199 int typ = mattype.
type ();
2203 (*current_liboctave_error_handler) (
"incorrect matrix type");
2206 double ainvnorm = 0.;
2216 atmp += fabs (
data (
i));
2243 cwork[
b.ridx (
i)] =
b.data (
i);
2255 goto triangular_error;
2261 i <
cidx (kidx+1)-1;
i++)
2264 cwork[iidx] = cwork[iidx] -
tmp *
data (
i);
2276 if (ii + new_nnz > x_nz)
2285 if (cwork[rperm[
i]] != 0.)
2288 retval.xdata (ii++) = cwork[rperm[
i]];
2293 retval.maybe_compress ();
2315 i <
cidx (iidx+1)-1;
i++)
2318 work[idx2] = work[idx2] -
tmp *
data (
i);
2325 atmp += fabs (work[
i]);
2328 if (atmp > ainvnorm)
2331 rcond = 1. / ainvnorm / anorm;
2343 cwork[
b.ridx (
i)] =
b.data (
i);
2353 goto triangular_error;
2361 cwork[iidx] = cwork[iidx] -
tmp *
data (
i);
2373 if (ii + new_nnz > x_nz)
2385 retval.xdata (ii++) = cwork[
i];
2390 retval.maybe_compress ();
2413 work[iidx] = work[iidx] -
tmp *
data (
i);
2420 atmp += fabs (work[
i]);
2423 if (atmp > ainvnorm)
2426 rcond = 1. / ainvnorm / anorm;
2435 sing_handler (rcond);
2442 volatile double rcond_plus_one = rcond + 1.0;
2450 sing_handler (rcond);
2464 solve_singularity_handler sing_handler,
2465 bool calc_cond)
const 2474 if (nr !=
b.rows ())
2476 (
"matrix dimension mismatch solution of linear equations");
2478 if (nr == 0 || nc == 0 ||
b.cols () == 0)
2483 int typ = mattype.
type ();
2487 (*current_liboctave_error_handler) (
"incorrect matrix type");
2490 double ainvnorm = 0.;
2501 atmp += fabs (
data (
i));
2519 work[perm[
i]] =
b(
i,j);
2529 if (perm[
ridx (
i)] < minr)
2531 minr = perm[
ridx (
i)];
2535 if (minr !=
k ||
data (mini) == 0)
2538 goto triangular_error;
2541 double tmp = work[
k] /
data (mini);
2549 work[iidx] = work[iidx] -
tmp *
data (
i);
2577 if (perm[
ridx (
i)] < minr)
2579 minr = perm[
ridx (
i)];
2583 double tmp = work[
k] /
data (mini);
2592 work[iidx] = work[iidx] -
tmp *
data (
i);
2600 atmp += fabs (work[
i]);
2603 if (atmp > ainvnorm)
2606 rcond = 1. / ainvnorm / anorm;
2627 goto triangular_error;
2636 work[iidx] = work[iidx] -
tmp *
data (
i);
2666 work[iidx] = work[iidx] -
tmp *
data (
i);
2673 atmp += fabs (work[
i]);
2676 if (atmp > ainvnorm)
2679 rcond = 1. / ainvnorm / anorm;
2688 sing_handler (rcond);
2695 volatile double rcond_plus_one = rcond + 1.0;
2703 sing_handler (rcond);
2717 solve_singularity_handler sing_handler,
2718 bool calc_cond)
const 2727 if (nr !=
b.rows ())
2729 (
"matrix dimension mismatch solution of linear equations");
2731 if (nr == 0 || nc == 0 ||
b.cols () == 0)
2736 int typ = mattype.
type ();
2740 (*current_liboctave_error_handler) (
"incorrect matrix type");
2743 double ainvnorm = 0.;
2753 atmp += fabs (
data (
i));
2776 work[perm[
b.ridx (
i)]] =
b.data (
i);
2786 if (perm[
ridx (
i)] < minr)
2788 minr = perm[
ridx (
i)];
2792 if (minr !=
k ||
data (mini) == 0)
2795 goto triangular_error;
2798 double tmp = work[
k] /
data (mini);
2806 work[iidx] = work[iidx] -
tmp *
data (
i);
2818 if (ii + new_nnz > x_nz)
2830 retval.xdata (ii++) = work[
i];
2835 retval.maybe_compress ();
2856 if (perm[
ridx (
i)] < minr)
2858 minr = perm[
ridx (
i)];
2862 double tmp = work[
k] /
data (mini);
2871 work[iidx] = work[iidx] -
tmp *
data (
i);
2879 atmp += fabs (work[
i]);
2882 if (atmp > ainvnorm)
2885 rcond = 1. / ainvnorm / anorm;
2897 work[
b.ridx (
i)] =
b.data (
i);
2906 goto triangular_error;
2914 work[iidx] = work[iidx] -
tmp *
data (
i);
2926 if (ii + new_nnz > x_nz)
2938 retval.xdata (ii++) = work[
i];
2943 retval.maybe_compress ();
2966 work[iidx] = work[iidx] -
tmp *
data (
i);
2973 atmp += fabs (work[
i]);
2976 if (atmp > ainvnorm)
2979 rcond = 1. / ainvnorm / anorm;
2988 sing_handler (rcond);
2995 volatile double rcond_plus_one = rcond + 1.0;
3003 sing_handler (rcond);
3017 solve_singularity_handler sing_handler,
3018 bool calc_cond)
const 3027 if (nr !=
b.rows ())
3029 (
"matrix dimension mismatch solution of linear equations");
3031 if (nr == 0 || nc == 0 ||
b.cols () == 0)
3036 int typ = mattype.
type ();
3040 (*current_liboctave_error_handler) (
"incorrect matrix type");
3043 double ainvnorm = 0.;
3054 atmp += fabs (
data (
i));
3071 cwork[perm[
i]] =
b(
i,j);
3081 if (perm[
ridx (
i)] < minr)
3083 minr = perm[
ridx (
i)];
3087 if (minr !=
k ||
data (mini) == 0)
3090 goto triangular_error;
3101 cwork[iidx] = cwork[iidx] -
tmp *
data (
i);
3130 if (perm[
ridx (
i)] < minr)
3132 minr = perm[
ridx (
i)];
3136 double tmp = work[
k] /
data (mini);
3145 work[iidx] = work[iidx] -
tmp *
data (
i);
3153 atmp += fabs (work[
i]);
3156 if (atmp > ainvnorm)
3159 rcond = 1. / ainvnorm / anorm;
3181 goto triangular_error;
3189 cwork[iidx] = cwork[iidx] -
tmp *
data (
i);
3220 work[iidx] = work[iidx] -
tmp *
data (
i);
3227 atmp += fabs (work[
i]);
3230 if (atmp > ainvnorm)
3233 rcond = 1. / ainvnorm / anorm;
3242 sing_handler (rcond);
3249 volatile double rcond_plus_one = rcond + 1.0;
3257 sing_handler (rcond);
3271 solve_singularity_handler sing_handler,
3272 bool calc_cond)
const 3281 if (nr !=
b.rows ())
3283 (
"matrix dimension mismatch solution of linear equations");
3285 if (nr == 0 || nc == 0 ||
b.cols () == 0)
3290 int typ = mattype.
type ();
3294 (*current_liboctave_error_handler) (
"incorrect matrix type");
3297 double ainvnorm = 0.;
3307 atmp += fabs (
data (
i));
3330 cwork[perm[
b.ridx (
i)]] =
b.data (
i);
3340 if (perm[
ridx (
i)] < minr)
3342 minr = perm[
ridx (
i)];
3346 if (minr !=
k ||
data (mini) == 0)
3349 goto triangular_error;
3360 cwork[iidx] = cwork[iidx] -
tmp *
data (
i);
3372 if (ii + new_nnz > x_nz)
3384 retval.xdata (ii++) = cwork[
i];
3389 retval.maybe_compress ();
3411 if (perm[
ridx (
i)] < minr)
3413 minr = perm[
ridx (
i)];
3417 double tmp = work[
k] /
data (mini);
3426 work[iidx] = work[iidx] -
tmp *
data (
i);
3434 atmp += fabs (work[
i]);
3437 if (atmp > ainvnorm)
3440 rcond = 1. / ainvnorm / anorm;
3452 cwork[
b.ridx (
i)] =
b.data (
i);
3461 goto triangular_error;
3469 cwork[iidx] = cwork[iidx] -
tmp *
data (
i);
3481 if (ii + new_nnz > x_nz)
3493 retval.xdata (ii++) = cwork[
i];
3498 retval.maybe_compress ();
3522 work[iidx] = work[iidx] -
tmp *
data (
i);
3529 atmp += fabs (work[
i]);
3532 if (atmp > ainvnorm)
3535 rcond = 1. / ainvnorm / anorm;
3544 sing_handler (rcond);
3551 volatile double rcond_plus_one = rcond + 1.0;
3559 sing_handler (rcond);
3573 solve_singularity_handler sing_handler,
3574 bool calc_cond)
const 3582 if (nr != nc || nr !=
b.rows ())
3584 (
"matrix dimension mismatch solution of linear equations");
3586 if (nr == 0 ||
b.cols () == 0)
3589 (*current_liboctave_error_handler)
3590 (
"calculation of condition number not implemented");
3594 volatile int typ = mattype.
type ();
3612 D[nc-1] =
data (ii);
3628 else if (
ridx (
i) == j + 1)
3633 F77_INT tmp_nr = octave::to_f77_int (nr);
3671 DL[j] =
data (ii++);
3672 DU[j] =
data (ii++);
3674 D[nc-1] =
data (ii);
3691 else if (
ridx (
i) == j + 1)
3693 else if (
ridx (
i) == j - 1)
3698 F77_INT tmp_nr = octave::to_f77_int (nr);
3720 sing_handler (rcond);
3730 (*current_liboctave_error_handler) (
"incorrect matrix type");
3739 solve_singularity_handler sing_handler,
3740 bool calc_cond)
const 3748 if (nr != nc || nr !=
b.rows ())
3750 (
"matrix dimension mismatch solution of linear equations");
3752 if (nr == 0 ||
b.cols () == 0)
3755 (*current_liboctave_error_handler)
3756 (
"calculation of condition number not implemented");
3760 int typ = mattype.
type ();
3781 DL[j] =
data (ii++);
3782 DU[j] =
data (ii++);
3784 D[nc-1] =
data (ii);
3801 else if (
ridx (
i) == j + 1)
3803 else if (
ridx (
i) == j - 1)
3808 F77_INT tmp_nr = octave::to_f77_int (nr);
3812 F77_XFCN (dgttrf, DGTTRF, (tmp_nr, DL, D, DU, DU2, pipvt, tmp_err));
3821 sing_handler (rcond);
3844 work[
b.ridx (
i)] =
b.data (
i);
3849 (F77_CONST_CHAR_ARG2 (&job, 1),
3850 tmp_nr, 1, DL, D, DU, DU2, pipvt,
3852 F77_CHAR_ARG_LEN (1)));
3863 if (ii + new_nnz > x_nz)
3875 retval.xdata (ii++) = work[
i];
3880 retval.maybe_compress ();
3884 (*current_liboctave_error_handler) (
"incorrect matrix type");
3893 solve_singularity_handler sing_handler,
3894 bool calc_cond)
const 3902 if (nr != nc || nr !=
b.rows ())
3904 (
"matrix dimension mismatch solution of linear equations");
3906 if (nr == 0 ||
b.cols () == 0)
3909 (*current_liboctave_error_handler)
3910 (
"calculation of condition number not implemented");
3914 volatile int typ = mattype.
type ();
3932 D[nc-1] =
data (ii);
3948 else if (
ridx (
i) == j + 1)
3953 F77_INT tmp_nr = octave::to_f77_int (nr);
3992 DL[j] =
data (ii++);
3993 DU[j] =
data (ii++);
3995 D[nc-1] =
data (ii);
4012 else if (
ridx (
i) == j + 1)
4014 else if (
ridx (
i) == j - 1)
4019 F77_INT tmp_nr = octave::to_f77_int (nr);
4046 sing_handler (rcond);
4054 (*current_liboctave_error_handler) (
"incorrect matrix type");
4063 solve_singularity_handler sing_handler,
4064 bool calc_cond)
const 4072 if (nr != nc || nr !=
b.rows ())
4074 (
"matrix dimension mismatch solution of linear equations");
4076 if (nr == 0 ||
b.cols () == 0)
4079 (*current_liboctave_error_handler)
4080 (
"calculation of condition number not implemented");
4084 int typ = mattype.
type ();
4105 DL[j] =
data (ii++);
4106 DU[j] =
data (ii++);
4108 D[nc-1] =
data (ii);
4125 else if (
ridx (
i) == j + 1)
4127 else if (
ridx (
i) == j - 1)
4132 F77_INT tmp_nr = octave::to_f77_int (nr);
4136 F77_XFCN (dgttrf, DGTTRF, (tmp_nr, DL, D, DU, DU2, pipvt, tmp_err));
4147 sing_handler (rcond);
4180 (F77_CONST_CHAR_ARG2 (&job, 1),
4181 tmp_nr, 1, DL, D, DU, DU2, pipvt,
4183 F77_CHAR_ARG_LEN (1)));
4190 (*current_liboctave_error_handler)
4191 (
"SparseMatrix::solve solve failed");
4198 (F77_CONST_CHAR_ARG2 (&job, 1),
4199 tmp_nr, 1, DL, D, DU, DU2, pipvt,
4201 F77_CHAR_ARG_LEN (1)));
4208 (*current_liboctave_error_handler)
4209 (
"SparseMatrix::solve solve failed");
4219 if (Bx[
i] != 0. || Bz[
i] != 0.)
4222 if (ii + new_nnz > x_nz)
4231 if (Bx[
i] != 0. || Bz[
i] != 0.)
4241 retval.maybe_compress ();
4245 (*current_liboctave_error_handler) (
"incorrect matrix type");
4254 solve_singularity_handler sing_handler,
4255 bool calc_cond)
const 4263 if (nr != nc || nr !=
b.rows ())
4265 (
"matrix dimension mismatch solution of linear equations");
4267 if (nr == 0 ||
b.cols () == 0)
4272 volatile int typ = mattype.
type ();
4288 tmp_data[ii++] = 0.;
4296 m_band(ri - j, j) =
data (
i);
4304 F77_INT tmp_nr = octave::to_f77_int (nr);
4309 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4310 tmp_nr, n_lower, tmp_data, ldm, tmp_err
4311 F77_CHAR_ARG_LEN (1)));
4334 (F77_CONST_CHAR_ARG2 (&job, 1),
4335 tmp_nr, n_lower, tmp_data, ldm,
4336 anorm, rcond, pz, piz, tmp_err
4337 F77_CHAR_ARG_LEN (1)));
4344 volatile double rcond_plus_one = rcond + 1.0;
4352 sing_handler (rcond);
4371 (F77_CONST_CHAR_ARG2 (&job, 1),
4372 tmp_nr, n_lower,
b_nc, tmp_data,
4374 F77_CHAR_ARG_LEN (1)));
4381 (*current_liboctave_error_handler)
4382 (
"SparseMatrix::solve solve failed");
4394 F77_INT ldm = n_upper + 2 * n_lower + 1;
4403 for (
F77_INT j = 0; j < ldm; j++)
4405 tmp_data[ii++] = 0.;
4410 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
4420 atmp += fabs (
data (
i));
4426 F77_INT tmp_nr = octave::to_f77_int (nr);
4433 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4434 tmp_data, ldm, pipvt, tmp_err));
4447 sing_handler (rcond);
4463 F77_INT tmp_nc = octave::to_f77_int (nc);
4466 (F77_CONST_CHAR_ARG2 (&job, 1),
4467 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4468 anorm, rcond, pz, piz, tmp_err
4469 F77_CHAR_ARG_LEN (1)));
4476 volatile double rcond_plus_one = rcond + 1.0;
4484 sing_handler (rcond);
4504 (F77_CONST_CHAR_ARG2 (&job, 1),
4505 tmp_nr, n_lower, n_upper,
b_nc, tmp_data,
4507 F77_CHAR_ARG_LEN (1)));
4514 (*current_liboctave_error_handler) (
"incorrect matrix type");
4523 solve_singularity_handler sing_handler,
4524 bool calc_cond)
const 4532 if (nr != nc || nr !=
b.rows ())
4534 (
"matrix dimension mismatch solution of linear equations");
4536 if (nr == 0 ||
b.cols () == 0)
4541 volatile int typ = mattype.
type ();
4547 F77_INT ldm = octave::to_f77_int (n_lower + 1);
4556 for (
F77_INT j = 0; j < ldm; j++)
4558 tmp_data[ii++] = 0.;
4566 m_band(ri - j, j) =
data (
i);
4574 F77_INT tmp_nr = octave::to_f77_int (nr);
4579 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4580 tmp_nr, n_lower, tmp_data, ldm, tmp_err
4581 F77_CHAR_ARG_LEN (1)));
4602 (F77_CONST_CHAR_ARG2 (&job, 1),
4603 tmp_nr, n_lower, tmp_data, ldm,
4604 anorm, rcond, pz, piz, tmp_err
4605 F77_CHAR_ARG_LEN (1)));
4612 volatile double rcond_plus_one = rcond + 1.0;
4620 sing_handler (rcond);
4646 Bx[
i] =
b.elem (
i, j);
4649 (F77_CONST_CHAR_ARG2 (&job, 1),
4650 tmp_nr, n_lower, 1, tmp_data,
4651 ldm, Bx,
b_nr, tmp_err
4652 F77_CHAR_ARG_LEN (1)));
4659 (*current_liboctave_error_handler)
4660 (
"SparseMatrix::solve solve failed");
4676 sz = x_nz + (
sz > 100 ?
sz : 100);
4687 retval.maybe_compress ();
4697 F77_INT ldm = octave::to_f77_int (n_upper + 2 * n_lower + 1);
4708 tmp_data[ii++] = 0.;
4713 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
4723 atmp += fabs (
data (
i));
4729 F77_INT tmp_nr = octave::to_f77_int (nr);
4736 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4737 tmp_data, ldm, pipvt, tmp_err));
4748 sing_handler (rcond);
4764 F77_INT tmp_nc = octave::to_f77_int (nc);
4767 (F77_CONST_CHAR_ARG2 (&job, 1),
4768 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4769 anorm, rcond, pz, piz, tmp_err
4770 F77_CHAR_ARG_LEN (1)));
4777 volatile double rcond_plus_one = rcond + 1.0;
4785 sing_handler (rcond);
4811 i <
b.cidx (j+1);
i++)
4812 work[
b.ridx (
i)] =
b.data (
i);
4817 (F77_CONST_CHAR_ARG2 (&job, 1),
4818 tmp_nr, n_lower, n_upper, 1, tmp_data,
4819 ldm, pipvt, work,
b_nr, tmp_err
4820 F77_CHAR_ARG_LEN (1)));
4831 if (ii + new_nnz > x_nz)
4843 retval.xdata (ii++) = work[
i];
4848 retval.maybe_compress ();
4853 (*current_liboctave_error_handler) (
"incorrect matrix type");
4862 solve_singularity_handler sing_handler,
4863 bool calc_cond)
const 4871 if (nr != nc || nr !=
b.rows ())
4873 (
"matrix dimension mismatch solution of linear equations");
4875 if (nr == 0 ||
b.cols () == 0)
4880 volatile int typ = mattype.
type ();
4895 for (
F77_INT j = 0; j < ldm; j++)
4897 tmp_data[ii++] = 0.;
4905 m_band(ri - j, j) =
data (
i);
4913 F77_INT tmp_nr = octave::to_f77_int (nr);
4918 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4919 tmp_nr, n_lower, tmp_data, ldm, tmp_err
4920 F77_CHAR_ARG_LEN (1)));
4943 (F77_CONST_CHAR_ARG2 (&job, 1),
4944 tmp_nr, n_lower, tmp_data, ldm,
4945 anorm, rcond, pz, piz, tmp_err
4946 F77_CHAR_ARG_LEN (1)));
4953 volatile double rcond_plus_one = rcond + 1.0;
4961 sing_handler (rcond);
4991 (F77_CONST_CHAR_ARG2 (&job, 1),
4992 tmp_nr, n_lower, 1, tmp_data,
4993 ldm, Bx,
b_nr, tmp_err
4994 F77_CHAR_ARG_LEN (1)));
5001 (*current_liboctave_error_handler)
5002 (
"SparseMatrix::solve solve failed");
5008 (F77_CONST_CHAR_ARG2 (&job, 1),
5009 tmp_nr, n_lower, 1, tmp_data,
5010 ldm, Bz,
b_nr, tmp_err
5011 F77_CHAR_ARG_LEN (1)));
5018 (*current_liboctave_error_handler)
5019 (
"SparseMatrix::solve solve failed");
5036 F77_INT ldm = n_upper + 2 * n_lower + 1;
5045 for (
F77_INT j = 0; j < ldm; j++)
5047 tmp_data[ii++] = 0.;
5052 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
5062 atmp += fabs (
data (
i));
5068 F77_INT tmp_nr = octave::to_f77_int (nr);
5075 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5076 tmp_data, ldm, pipvt, tmp_err));
5087 sing_handler (rcond);
5103 F77_INT tmp_nc = octave::to_f77_int (nc);
5106 (F77_CONST_CHAR_ARG2 (&job, 1),
5107 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5108 anorm, rcond, pz, piz, tmp_err
5109 F77_CHAR_ARG_LEN (1)));
5116 volatile double rcond_plus_one = rcond + 1.0;
5124 sing_handler (rcond);
5154 (F77_CONST_CHAR_ARG2 (&job, 1),
5155 tmp_nr, n_lower, n_upper, 1, tmp_data,
5156 ldm, pipvt, Bx,
b_nr, tmp_err
5157 F77_CHAR_ARG_LEN (1)));
5162 (F77_CONST_CHAR_ARG2 (&job, 1),
5163 tmp_nr, n_lower, n_upper, 1, tmp_data,
5164 ldm, pipvt, Bz,
b_nr, tmp_err
5165 F77_CHAR_ARG_LEN (1)));
5176 (*current_liboctave_error_handler) (
"incorrect matrix type");
5185 solve_singularity_handler sing_handler,
5186 bool calc_cond)
const 5194 if (nr != nc || nr !=
b.rows ())
5196 (
"matrix dimension mismatch solution of linear equations");
5198 if (nr == 0 ||
b.cols () == 0)
5203 volatile int typ = mattype.
type ();
5218 for (
F77_INT j = 0; j < ldm; j++)
5220 tmp_data[ii++] = 0.;
5228 m_band(ri - j, j) =
data (
i);
5236 F77_INT tmp_nr = octave::to_f77_int (nr);
5241 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5242 tmp_nr, n_lower, tmp_data, ldm, tmp_err
5243 F77_CHAR_ARG_LEN (1)));
5267 (F77_CONST_CHAR_ARG2 (&job, 1),
5268 tmp_nr, n_lower, tmp_data, ldm,
5269 anorm, rcond, pz, piz, tmp_err
5270 F77_CHAR_ARG_LEN (1)));
5277 volatile double rcond_plus_one = rcond + 1.0;
5285 sing_handler (rcond);
5320 (F77_CONST_CHAR_ARG2 (&job, 1),
5321 tmp_nr, n_lower, 1, tmp_data,
5322 ldm, Bx,
b_nr, tmp_err
5323 F77_CHAR_ARG_LEN (1)));
5330 (*current_liboctave_error_handler)
5331 (
"SparseMatrix::solve solve failed");
5337 (F77_CONST_CHAR_ARG2 (&job, 1),
5338 tmp_nr, n_lower, 1, tmp_data,
5339 ldm, Bz,
b_nr, tmp_err
5340 F77_CHAR_ARG_LEN (1)));
5347 (*current_liboctave_error_handler)
5348 (
"SparseMatrix::solve solve failed");
5358 if (Bx[
i] != 0. || Bz[
i] != 0.)
5361 if (ii + new_nnz > x_nz)
5370 if (Bx[
i] != 0. || Bz[
i] != 0.)
5380 retval.maybe_compress ();
5390 F77_INT ldm = n_upper + 2 * n_lower + 1;
5399 for (
F77_INT j = 0; j < ldm; j++)
5401 tmp_data[ii++] = 0.;
5406 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
5416 atmp += fabs (
data (
i));
5422 F77_INT tmp_nr = octave::to_f77_int (nr);
5429 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5430 tmp_data, ldm, pipvt, tmp_err));
5441 sing_handler (rcond);
5457 F77_INT tmp_nc = octave::to_f77_int (nc);
5460 (F77_CONST_CHAR_ARG2 (&job, 1),
5461 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5462 anorm, rcond, pz, piz, tmp_err
5463 F77_CHAR_ARG_LEN (1)));
5470 volatile double rcond_plus_one = rcond + 1.0;
5478 sing_handler (rcond);
5509 i <
b.cidx (j+1);
i++)
5512 Bx[
b.ridx (
i)] =
c.real ();
5513 Bz[
b.ridx (
i)] =
c.imag ();
5517 (F77_CONST_CHAR_ARG2 (&job, 1),
5518 tmp_nr, n_lower, n_upper, 1, tmp_data,
5519 ldm, pipvt, Bx,
b_nr, tmp_err
5520 F77_CHAR_ARG_LEN (1)));
5525 (F77_CONST_CHAR_ARG2 (&job, 1),
5526 tmp_nr, n_lower, n_upper, 1, tmp_data,
5527 ldm, pipvt, Bz,
b_nr, tmp_err
5528 F77_CHAR_ARG_LEN (1)));
5536 if (Bx[
i] != 0. || Bz[
i] != 0.)
5539 if (ii + new_nnz > x_nz)
5548 if (Bx[
i] != 0. || Bz[
i] != 0.)
5557 retval.maybe_compress ();
5562 (*current_liboctave_error_handler) (
"incorrect matrix type");
5570 Matrix& Info, solve_singularity_handler sing_handler,
5571 bool calc_cond)
const 5574 void *Numeric =
nullptr;
5577 #if defined (HAVE_UMFPACK) 5580 Control =
Matrix (UMFPACK_CONTROL, 1);
5586 Control (UMFPACK_PRL) =
tmp;
5590 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
5591 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
5597 Control (UMFPACK_FIXQ) =
tmp;
5603 const double *Ax =
data ();
5613 Info =
Matrix (1, UMFPACK_INFO);
5618 Ax,
nullptr, &Symbolic, control, info);
5628 (*current_liboctave_error_handler)
5629 (
"SparseMatrix::solve symbolic factorization failed");
5638 Ax, Symbolic, &Numeric, control, info);
5642 rcond = Info (UMFPACK_RCOND);
5645 volatile double rcond_plus_one = rcond + 1.0;
5647 if (status == UMFPACK_WARNING_singular_matrix
5655 sing_handler (rcond);
5659 else if (status < 0)
5665 (*current_liboctave_error_handler)
5666 (
"SparseMatrix::solve numeric factorization failed");
5681 octave_unused_parameter (rcond);
5682 octave_unused_parameter (Control);
5683 octave_unused_parameter (Info);
5684 octave_unused_parameter (sing_handler);
5685 octave_unused_parameter (calc_cond);
5687 (*current_liboctave_error_handler)
5688 (
"support for UMFPACK was unavailable or disabled " 5689 "when liboctave was built");
5699 solve_singularity_handler sing_handler,
5700 bool calc_cond)
const 5708 if (nr != nc || nr !=
b.rows ())
5710 (
"matrix dimension mismatch solution of linear equations");
5712 if (nr == 0 ||
b.cols () == 0)
5717 volatile int typ = mattype.
type ();
5722 #if defined (HAVE_CHOLMOD) 5723 cholmod_common Common;
5724 cholmod_common *cm = &Common;
5728 cm->prefer_zomplex =
false;
5734 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
5738 cm->print =
static_cast<int> (spu) + 2;
5739 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5743 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5744 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5746 cm->final_ll =
true;
5748 cholmod_sparse Astore;
5749 cholmod_sparse *
A = &Astore;
5760 #if defined (OCTAVE_ENABLE_64) 5761 A->itype = CHOLMOD_LONG;
5763 A->itype = CHOLMOD_INT;
5765 A->dtype = CHOLMOD_DOUBLE;
5767 A->xtype = CHOLMOD_REAL;
5770 if (
A->x ==
nullptr)
5773 cholmod_dense Bstore;
5774 cholmod_dense *
B = &Bstore;
5775 B->nrow =
b.rows ();
5776 B->ncol =
b.cols ();
5778 B->nzmax =
B->nrow *
B->ncol;
5779 B->dtype = CHOLMOD_DOUBLE;
5780 B->xtype = CHOLMOD_REAL;
5782 B->x =
const_cast<double *
>(
b.fortran_vec ());
5783 if (
B->x ==
nullptr)
5805 volatile double rcond_plus_one = rcond + 1.0;
5813 sing_handler (rcond);
5832 retval.xelem (
i,j) =
static_cast<double *
>(X->x)[jr +
i];
5839 static char blank_name[] =
" ";
5844 (*current_liboctave_warning_with_id_handler)
5845 (
"Octave:missing-dependency",
5846 "support for CHOLMOD was unavailable or disabled " 5847 "when liboctave was built");
5856 #if defined (HAVE_UMFPACK) 5859 factorize (
err, rcond, Control, Info, sing_handler, calc_cond);
5864 Control (UMFPACK_IRSTEP) = 1;
5865 const double *Bx =
b.fortran_vec ();
5875 const double *Ax =
data ();
5882 Ax, &
result[iidx], &Bx[iidx],
5883 Numeric, control, info);
5889 (*current_liboctave_error_handler)
5890 (
"SparseMatrix::solve solve failed");
5905 octave_unused_parameter (rcond);
5906 octave_unused_parameter (sing_handler);
5907 octave_unused_parameter (calc_cond);
5909 (*current_liboctave_error_handler)
5910 (
"support for UMFPACK was unavailable or disabled " 5911 "when liboctave was built");
5915 (*current_liboctave_error_handler) (
"incorrect matrix type");
5924 solve_singularity_handler sing_handler,
5925 bool calc_cond)
const 5933 if (nr != nc || nr !=
b.rows ())
5935 (
"matrix dimension mismatch solution of linear equations");
5937 if (nr == 0 ||
b.cols () == 0)
5942 volatile int typ = mattype.
type ();
5947 #if defined (HAVE_CHOLMOD) 5948 cholmod_common Common;
5949 cholmod_common *cm = &Common;
5953 cm->prefer_zomplex =
false;
5959 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
5963 cm->print =
static_cast<int> (spu) + 2;
5964 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5968 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5969 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5971 cm->final_ll =
true;
5973 cholmod_sparse Astore;
5974 cholmod_sparse *
A = &Astore;
5985 #if defined (OCTAVE_ENABLE_64) 5986 A->itype = CHOLMOD_LONG;
5988 A->itype = CHOLMOD_INT;
5990 A->dtype = CHOLMOD_DOUBLE;
5992 A->xtype = CHOLMOD_REAL;
5995 if (
A->x ==
nullptr)
5998 cholmod_sparse Bstore;
5999 cholmod_sparse *
B = &Bstore;
6000 B->nrow =
b.rows ();
6001 B->ncol =
b.cols ();
6004 B->nzmax =
b.nnz ();
6008 #if defined (OCTAVE_ENABLE_64) 6009 B->itype = CHOLMOD_LONG;
6011 B->itype = CHOLMOD_INT;
6013 B->dtype = CHOLMOD_DOUBLE;
6015 B->xtype = CHOLMOD_REAL;
6018 if (
B->x ==
nullptr)
6039 volatile double rcond_plus_one = rcond + 1.0;
6047 sing_handler (rcond);
6062 static_cast<octave_idx_type>(X->ncol),
6063 static_cast<octave_idx_type>(X->nzmax));
6065 j <= static_cast<octave_idx_type>(X->ncol); j++)
6068 j < static_cast<octave_idx_type>(X->nzmax); j++)
6071 retval.xdata (j) =
static_cast<double *
>(X->x)[j];
6078 static char blank_name[] =
" ";
6083 (*current_liboctave_warning_with_id_handler)
6084 (
"Octave:missing-dependency",
6085 "support for CHOLMOD was unavailable or disabled " 6086 "when liboctave was built");
6095 #if defined (HAVE_UMFPACK) 6098 sing_handler, calc_cond);
6103 Control (UMFPACK_IRSTEP) = 1;
6111 const double *Ax =
data ();
6127 Bx[
i] =
b.elem (
i, j);
6132 Ax, Xx, Bx, Numeric,
6139 (*current_liboctave_error_handler)
6140 (
"SparseMatrix::solve solve failed");
6157 sz = x_nz + (
sz > 100 ?
sz : 100);
6168 retval.maybe_compress ();
6178 octave_unused_parameter (rcond);
6179 octave_unused_parameter (sing_handler);
6180 octave_unused_parameter (calc_cond);
6182 (*current_liboctave_error_handler)
6183 (
"support for UMFPACK was unavailable or disabled " 6184 "when liboctave was built");
6188 (*current_liboctave_error_handler) (
"incorrect matrix type");
6197 solve_singularity_handler sing_handler,
6198 bool calc_cond)
const 6206 if (nr != nc || nr !=
b.rows ())
6208 (
"matrix dimension mismatch solution of linear equations");
6210 if (nr == 0 ||
b.cols () == 0)
6215 volatile int typ = mattype.
type ();
6220 #if defined (HAVE_CHOLMOD) 6221 cholmod_common Common;
6222 cholmod_common *cm = &Common;
6226 cm->prefer_zomplex =
false;
6232 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6236 cm->print =
static_cast<int> (spu) + 2;
6237 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6241 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6242 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6244 cm->final_ll =
true;
6246 cholmod_sparse Astore;
6247 cholmod_sparse *
A = &Astore;
6258 #if defined (OCTAVE_ENABLE_64) 6259 A->itype = CHOLMOD_LONG;
6261 A->itype = CHOLMOD_INT;
6263 A->dtype = CHOLMOD_DOUBLE;
6265 A->xtype = CHOLMOD_REAL;
6268 if (
A->x ==
nullptr)
6271 cholmod_dense Bstore;
6272 cholmod_dense *
B = &Bstore;
6273 B->nrow =
b.rows ();
6274 B->ncol =
b.cols ();
6276 B->nzmax =
B->nrow *
B->ncol;
6277 B->dtype = CHOLMOD_DOUBLE;
6278 B->xtype = CHOLMOD_COMPLEX;
6280 B->x =
const_cast<Complex *
>(
b.fortran_vec ());
6281 if (
B->x ==
nullptr)
6302 volatile double rcond_plus_one = rcond + 1.0;
6310 sing_handler (rcond);
6336 static char blank_name[] =
" ";
6341 (*current_liboctave_warning_with_id_handler)
6342 (
"Octave:missing-dependency",
6343 "support for CHOLMOD was unavailable or disabled " 6344 "when liboctave was built");
6353 #if defined (HAVE_UMFPACK) 6356 sing_handler, calc_cond);
6361 Control (UMFPACK_IRSTEP) = 1;
6369 const double *Ax =
data ();
6391 Ax, Xx, Bx, Numeric,
6397 Numeric, control, info);
6399 if (status < 0 || status2 < 0)
6404 (*current_liboctave_error_handler)
6405 (
"SparseMatrix::solve solve failed");
6423 octave_unused_parameter (rcond);
6424 octave_unused_parameter (sing_handler);
6425 octave_unused_parameter (calc_cond);
6427 (*current_liboctave_error_handler)
6428 (
"support for UMFPACK was unavailable or disabled " 6429 "when liboctave was built");
6433 (*current_liboctave_error_handler) (
"incorrect matrix type");
6442 solve_singularity_handler sing_handler,
6443 bool calc_cond)
const 6451 if (nr != nc || nr !=
b.rows ())
6453 (
"matrix dimension mismatch solution of linear equations");
6455 if (nr == 0 ||
b.cols () == 0)
6460 volatile int typ = mattype.
type ();
6465 #if defined (HAVE_CHOLMOD) 6466 cholmod_common Common;
6467 cholmod_common *cm = &Common;
6471 cm->prefer_zomplex =
false;
6477 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6481 cm->print =
static_cast<int> (spu) + 2;
6482 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6486 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6487 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6489 cm->final_ll =
true;
6491 cholmod_sparse Astore;
6492 cholmod_sparse *
A = &Astore;
6503 #if defined (OCTAVE_ENABLE_64) 6504 A->itype = CHOLMOD_LONG;
6506 A->itype = CHOLMOD_INT;
6508 A->dtype = CHOLMOD_DOUBLE;
6510 A->xtype = CHOLMOD_REAL;
6513 if (
A->x ==
nullptr)
6516 cholmod_sparse Bstore;
6517 cholmod_sparse *
B = &Bstore;
6518 B->nrow =
b.rows ();
6519 B->ncol =
b.cols ();
6522 B->nzmax =
b.nnz ();
6526 #if defined (OCTAVE_ENABLE_64) 6527 B->itype = CHOLMOD_LONG;
6529 B->itype = CHOLMOD_INT;
6531 B->dtype = CHOLMOD_DOUBLE;
6533 B->xtype = CHOLMOD_COMPLEX;
6536 if (
B->x ==
nullptr)
6557 volatile double rcond_plus_one = rcond + 1.0;
6565 sing_handler (rcond);
6580 (static_cast<octave_idx_type>(X->nrow),
6581 static_cast<octave_idx_type>(X->ncol),
6582 static_cast<octave_idx_type>(X->nzmax));
6584 j <= static_cast<octave_idx_type>(X->ncol); j++)
6587 j < static_cast<octave_idx_type>(X->nzmax); j++)
6597 static char blank_name[] =
" ";
6602 (*current_liboctave_warning_with_id_handler)
6603 (
"Octave:missing-dependency",
6604 "support for CHOLMOD was unavailable or disabled " 6605 "when liboctave was built");
6614 #if defined (HAVE_UMFPACK) 6617 sing_handler, calc_cond);
6622 Control (UMFPACK_IRSTEP) = 1;
6630 const double *Ax =
data ();
6657 Ax, Xx, Bx, Numeric,
6663 Numeric, control, info);
6665 if (status < 0 || status2 < 0)
6670 (*current_liboctave_error_handler)
6671 (
"SparseMatrix::solve solve failed");
6688 sz = x_nz + (
sz > 100 ?
sz : 100);
6699 retval.maybe_compress ();
6708 octave_unused_parameter (rcond);
6709 octave_unused_parameter (sing_handler);
6710 octave_unused_parameter (calc_cond);
6712 (*current_liboctave_error_handler)
6713 (
"support for UMFPACK was unavailable or disabled " 6714 "when liboctave was built");
6718 (*current_liboctave_error_handler) (
"incorrect matrix type");
6729 return solve (mattype,
b, info, rcond,
nullptr);
6737 return solve (mattype,
b, info, rcond,
nullptr);
6744 return solve (mattype,
b, info, rcond,
nullptr);
6749 double& rcond, solve_singularity_handler sing_handler,
6750 bool singular_fallback)
const 6753 int typ = mattype.
type (
false);
6756 typ = mattype.
type (*
this);
6773 (*current_liboctave_error_handler) (
"unknown matrix type");
6779 #if defined (USE_QRSOLVE) 6794 return solve (mattype,
b, info, rcond,
nullptr);
6802 return solve (mattype,
b, info, rcond,
nullptr);
6809 return solve (mattype,
b, info, rcond,
nullptr);
6815 solve_singularity_handler sing_handler,
6816 bool singular_fallback)
const 6819 int typ = mattype.
type (
false);
6822 typ = mattype.
type (*
this);
6838 (*current_liboctave_error_handler) (
"unknown matrix type");
6843 #if defined (USE_QRSOLVE) 6859 return solve (mattype,
b, info, rcond,
nullptr);
6867 return solve (mattype,
b, info, rcond,
nullptr);
6874 return solve (mattype,
b, info, rcond,
nullptr);
6880 solve_singularity_handler sing_handler,
6881 bool singular_fallback)
const 6884 int typ = mattype.
type (
false);
6887 typ = mattype.
type (*
this);
6903 (*current_liboctave_error_handler) (
"unknown matrix type");
6908 #if defined (USE_QRSOLVE) 6924 return solve (mattype,
b, info, rcond,
nullptr);
6932 return solve (mattype,
b, info, rcond,
nullptr);
6939 return solve (mattype,
b, info, rcond,
nullptr);
6945 solve_singularity_handler sing_handler,
6946 bool singular_fallback)
const 6949 int typ = mattype.
type (
false);
6952 typ = mattype.
type (*
this);
6968 (*current_liboctave_error_handler) (
"unknown matrix type");
6973 #if defined (USE_QRSOLVE) 6988 return solve (mattype,
b, info, rcond);
6996 return solve (mattype,
b, info, rcond);
7003 return solve (mattype,
b, info, rcond,
nullptr);
7009 solve_singularity_handler sing_handler)
const 7012 return solve (mattype,
tmp, info, rcond,
7013 sing_handler).
column (static_cast<octave_idx_type> (0));
7021 return solve (mattype,
b, info, rcond,
nullptr);
7029 return solve (mattype,
b, info, rcond,
nullptr);
7035 double& rcond)
const 7037 return solve (mattype,
b, info, rcond,
nullptr);
7043 solve_singularity_handler sing_handler)
const 7046 return solve (mattype,
tmp, info, rcond,
7047 sing_handler).
column (static_cast<octave_idx_type> (0));
7055 return solve (
b, info, rcond,
nullptr);
7062 return solve (
b, info, rcond,
nullptr);
7067 double& rcond)
const 7069 return solve (
b, info, rcond,
nullptr);
7074 solve_singularity_handler sing_handler)
const 7077 return solve (mattype,
b,
err, rcond, sing_handler);
7085 return solve (
b, info, rcond,
nullptr);
7093 return solve (
b, info, rcond,
nullptr);
7100 return solve (
b, info, rcond,
nullptr);
7105 solve_singularity_handler sing_handler)
const 7108 return solve (mattype,
b,
err, rcond, sing_handler);
7115 return solve (
b, info, rcond,
nullptr);
7120 double& rcond)
const 7122 return solve (
b, info, rcond,
nullptr);
7128 solve_singularity_handler sing_handler)
const 7131 return solve (mattype,
b,
err, rcond, sing_handler);
7139 return solve (
b, info, rcond,
nullptr);
7146 return solve (
b, info, rcond,
nullptr);
7151 double& rcond)
const 7153 return solve (
b, info, rcond,
nullptr);
7159 solve_singularity_handler sing_handler)
const 7162 return solve (mattype,
b,
err, rcond, sing_handler);
7169 return solve (
b, info, rcond);
7176 return solve (
b, info, rcond);
7181 double& rcond)
const 7183 return solve (
b, info, rcond,
nullptr);
7189 solve_singularity_handler sing_handler)
const 7193 sing_handler).
column (static_cast<octave_idx_type> (0));
7201 return solve (
b, info, rcond,
nullptr);
7208 return solve (
b, info, rcond,
nullptr);
7213 double& rcond)
const 7215 return solve (
b, info, rcond,
nullptr);
7221 solve_singularity_handler sing_handler)
const 7225 sing_handler).
column (static_cast<octave_idx_type> (0));
7289 if (
val != 0.0 &&
val != 1.0)
7425 if ((
rows () == 1 && dim == -1) || dim == 1)
7430 (
cidx (j+1) -
cidx (j) < nr ? 0.0 : 1.0), 1.0);
7444 double d = data (i); \ 7445 tmp[ridx (i)] += d * d 7448 double d = data (i); \ 7495 os <<
a.ridx (
i) + 1 <<
' ' << j + 1 <<
' ';
7509 return read_sparse_matrix<elt_type> (
is,
a, octave_read_value<double>);
7573 return do_mul_dm_sm<SparseMatrix> (
d,
a);
7579 return do_mul_sm_dm<SparseMatrix> (
a,
d);
7585 return do_add_dm_sm<SparseMatrix> (
d,
a);
7591 return do_sub_dm_sm<SparseMatrix> (
d,
a);
7597 return do_add_sm_dm<SparseMatrix> (
a,
d);
7603 return do_sub_sm_dm<SparseMatrix> (
a,
d);
7622 #define EMPTY_RETURN_CHECK(T) \ 7623 if (nr == 0 || nc == 0) \ 7709 bool ja_lt_max = ja < ja_max;
7713 bool jb_lt_max = jb < jb_max;
7715 while (ja_lt_max || jb_lt_max)
7718 if ((! jb_lt_max) || (ja_lt_max && (
a.ridx (ja) <
b.ridx (jb))))
7723 r.
ridx (jx) =
a.ridx (ja);
7728 ja_lt_max= ja < ja_max;
7730 else if ((! ja_lt_max)
7731 || (jb_lt_max && (
b.ridx (jb) <
a.ridx (ja))))
7736 r.
ridx (jx) =
b.ridx (jb);
7741 jb_lt_max= jb < jb_max;
7749 r.
ridx (jx) =
a.ridx (ja);
7753 ja_lt_max= ja < ja_max;
7755 jb_lt_max= jb < jb_max;
7859 bool ja_lt_max = ja < ja_max;
7863 bool jb_lt_max = jb < jb_max;
7865 while (ja_lt_max || jb_lt_max)
7868 if ((! jb_lt_max) || (ja_lt_max && (
a.ridx (ja) <
b.ridx (jb))))
7873 r.
ridx (jx) =
a.ridx (ja);
7878 ja_lt_max= ja < ja_max;
7880 else if ((! ja_lt_max)
7881 || (jb_lt_max && (
b.ridx (jb) <
a.ridx (ja))))
7886 r.
ridx (jx) =
b.ridx (jb);
7891 jb_lt_max= jb < jb_max;
7899 r.
ridx (jx) =
a.ridx (ja);
7903 ja_lt_max= ja < ja_max;
7905 jb_lt_max= jb < jb_max;
template Matrix dmsolve< Matrix, SparseMatrix, Matrix >(const SparseMatrix &, const Matrix &, octave_idx_type &)
bool any_element_is_nan(void) const
void octave_write_double(std::ostream &os, double d)
octave_idx_type rows(void) const
void resize(octave_idx_type r, octave_idx_type c)
SparseBoolMatrix any(int dim=-1) const
is already an absolute the name is checked against the file system instead of Octave s loadpath In this if otherwise an empty string is returned If the first argument is a cell array of search each directory of the loadpath for element of the cell array and return the first that matches If the second optional argument return a cell array containing the list of all files that have the same name in the path If no files are found
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
SparseMatrix prod(int dim=-1) const
SparseMatrix concat(const SparseMatrix &rb, const Array< octave_idx_type > &ra_idx)
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
bool too_large_for_float(void) const
MSparse< T > & insert(const Sparse< T > &a, octave_idx_type r, octave_idx_type c)
#define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN)
const octave_base_value const Array< octave_idx_type > & ra_idx
SparseMatrix max(double d, const SparseMatrix &m)
SparseMatrix sumsq(int dim=-1) const
RowVector row(octave_idx_type i) const
#define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, INIT_VAL, MT_RESULT)
RowVector row(octave_idx_type i) const
template SparseMatrix dmsolve< SparseMatrix, SparseMatrix, SparseMatrix >(const SparseMatrix &, const SparseMatrix &, octave_idx_type &)
identity matrix If supplied two scalar respectively For allows like xample val
void SparseCholError(int status, char *file, int line, char *message)
bool ishermitian(void) const
void * factorize(octave_idx_type &err, double &rcond, Matrix &Control, Matrix &Info, solve_singularity_handler sing_handler, bool calc_cond=false) const
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
octave_idx_type columns(void) const
#define F77_DBLE_CMPLX_ARG(x)
SparseMatrix sum(int dim=-1) const
SM octinternal_do_mul_pm_sm(const PermMatrix &p, const SM &a)
octave_idx_type * cidx(void)
const T * fortran_vec(void) const
Matrix trisolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
Matrix mul_trans(const Matrix &m, const SparseMatrix &a)
octave_idx_type * triangular_perm(void) const
#define SPARSE_SMSM_BOOL_OPS(M1, M2, ZERO)
SparseMatrix transpose(void) const
template SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseMatrix, SparseComplexMatrix >(const SparseMatrix &, const SparseComplexMatrix &, octave_idx_type &)
SparseMatrix & insert(const SparseMatrix &a, octave_idx_type r, octave_idx_type c)
DET determinant(void) const
#define lo_ieee_signbit(x)
SM octinternal_do_mul_sm_pm(const SM &a, const PermMatrix &p)
#define SPARSE_FULL_TRANS_MUL(RET_TYPE, EL_TYPE, ZERO, CONJ_OP)
#define UMFPACK_DNAME(name)
bool test_any(F fcn) const
Matrix sum(int dim=-1) const
SparseMatrix cumsum(int dim=-1) const
T & elem(octave_idx_type n)
#define SPARSE_REDUCTION_OP(RET_TYPE, EL_TYPE, OP, INIT_VAL, MT_RESULT)
std::istream & operator>>(std::istream &is, SparseMatrix &a)
MSparse< T > diag(octave_idx_type k=0) const
bool any_element_not_one_or_zero(void) const
#define FULL_SPARSE_MUL(RET_TYPE, EL_TYPE, ZERO)
SparseMatrix Pr(void) const
nd example oindent opens the file binary numeric values will be read assuming they are stored in IEEE format with the least significant bit and then converted to the native representation Opening a file that is already open simply opens it again and returns a separate file id It is not an error to open a file several though writing to the same file through several different file ids may produce unexpected results The possible values of text mode reading and writing automatically converts linefeeds to the appropriate line end character for the you may append a you must also open the file in binary mode The parameter conversions are currently only supported for and permissions will be set to and then everything is written in a single operation This is very efficient and improves performance c
#define F77_XFCN(f, F, args)
SparseMatrix max(int dim=-1) const
void err_nan_to_logical_conversion(void)
static double get_key(const std::string &key)
MSparse< T > reshape(const dim_vector &new_dims) const
Matrix bsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
#define SPARSE_SMS_CMP_OPS(M, MZ, CM, S, SZ, CS)
#define SPARSE_FULL_MUL(RET_TYPE, EL_TYPE, ZERO)
F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
SparseMatrix reshape(const dim_vector &new_dims) const
octave_idx_type cols(void) const
octave_value resize(const dim_vector &dv, bool fill=false) const
SparseMatrix ipermute(const Array< octave_idx_type > &vec) const
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
SparseMatrix dinverse(MatrixType &mattype, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
#define SPARSE_ALL_OP(DIM)
bool any_element_is_negative(bool=false) const
octave_idx_type nnz(void) const
Actual number of nonzero terms.
SparseMatrix abs(void) const
SparseMatrix imag(const SparseComplexMatrix &a)
octave_idx_type nnz(void) const
Count nonzero elements.
#define SPARSE_ANY_OP(DIM)
bool all_integers(double &max_val, double &min_val) const
MSparse< T > squeeze(void) const
Matrix fsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
int type(bool quiet=true)
OCTAVE_EXPORT octave_value_list isdir nd deftypefn *std::string nm
Matrix trans_mul(const SparseMatrix &m, const Matrix &a)
SparseMatrix cumprod(int dim=-1) const
F77_RET_T const F77_INT F77_CMPLX * A
#define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN)
bool all_elements_are_zero(void) const
double & elem(octave_idx_type n)
#define FULL_SPARSE_MUL_TRANS(RET_TYPE, EL_TYPE, ZERO, CONJ_OP)
#define CHOLMOD_NAME(name)
SparseMatrix Pc(void) const
SparseMatrix inverse(void) const
#define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE)
int first_non_singleton(int def=0) const
SparseMatrix squeeze(void) const
void resize(const dim_vector &dv, const T &rfv)
Resizing (with fill).
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
void mark_as_unsymmetric(void)
SparseMatrix diag(octave_idx_type k=0) const
octave_idx_type * ridx(void)
Matrix utsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
SparseMatrix real(const SparseComplexMatrix &a)
ColumnVector column(octave_idx_type i) const
SparseMatrix operator-(const DiagMatrix &d, const SparseMatrix &a)
ColumnVector column(octave_idx_type i) const
SparseMatrix min(double d, const SparseMatrix &m)
Array< T > array_value(void) const
SparseMatrix operator+(const DiagMatrix &d, const SparseMatrix &a)
void mark_as_rectangular(void)
#define EMPTY_RETURN_CHECK(T)
SparseMatrix min(int dim=-1) const
With real return the complex result
#define SPARSE_SMSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2)
MSparse< T > permute(const Array< octave_idx_type > &vec, bool inv=false) const
octave_idx_type cols(void) const
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
int SparseCholPrint(const char *fmt,...)
bool all_elements_are_int_or_inf_or_nan(void) const
bool is_dense(void) const
#define END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
Matrix matrix_value(void) const
octave_f77_int_type F77_INT
#define SPARSE_SSM_CMP_OPS(S, SZ, SC, M, MZ, MC)
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol NaN(Not a Number). NaN is the result of operations which do not produce a well defined 0 result. Common operations which produce a NaN are arithmetic with infinity ex($\infty - \infty$)
bool issymmetric(void) const
bool operator==(const SparseMatrix &a) const
Matrix dsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
bool any_element_is_inf_or_nan(void) const
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
SparseMatrix operator*(const SparseMatrix &m, const SparseMatrix &a)
Sparse< T > maybe_compress(bool remove_zeros=false)
Matrix solve(MatrixType &typ, const Matrix &b) const
template ComplexMatrix dmsolve< ComplexMatrix, SparseMatrix, ComplexMatrix >(const SparseMatrix &, const ComplexMatrix &, octave_idx_type &)
bool xtoo_large_for_float(double x)
MSparse< T > ipermute(const Array< octave_idx_type > &vec) const
SparseBoolMatrix operator!(void) const
dim_vector dims(void) const
MatrixType transpose(void) const
OCTAVE_EXPORT octave_value_list error nd deftypefn *const octave_scalar_map err
#define SPARSE_SMS_BOOL_OPS(M, S, ZERO)
#define BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
octave_idx_type ndims(void) const
Number of dimensions.
SparseMatrix Q(void) const
std::complex< double > Complex
#define SPARSE_SSM_BOOL_OPS(S, M, ZERO)
write the output to stdout if nargout is
Matrix ltsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
Vector representing the dimensions (size) of an Array.
void warn_singular_matrix(double rcond)
SparseMatrix permute(const Array< octave_idx_type > &vec, bool inv=false) const
SparseBoolMatrix all(int dim=-1) const
octave_idx_type rows(void) const
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
bool operator!=(const SparseMatrix &a) const
SparseMatrix tinverse(MatrixType &mattype, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
std::ostream & operator<<(std::ostream &os, const SparseMatrix &a)