73 const octave_idx_type&,
const octave_idx_type&,
74 double*,
const octave_idx_type&,
75 octave_idx_type*, octave_idx_type&);
79 const octave_idx_type&,
const octave_idx_type&,
80 const octave_idx_type&,
const octave_idx_type&,
81 const double*,
const octave_idx_type&,
82 const octave_idx_type*,
double*,
83 const octave_idx_type&, octave_idx_type&
88 const octave_idx_type&,
const octave_idx_type&,
89 const octave_idx_type&,
double*,
90 const octave_idx_type&,
const octave_idx_type*,
91 const double&,
double&,
double*,
92 octave_idx_type*, octave_idx_type&
97 const octave_idx_type&,
const octave_idx_type&,
98 double*,
const octave_idx_type&, octave_idx_type&
103 const octave_idx_type&,
const octave_idx_type&,
104 const octave_idx_type&,
double*,
105 const octave_idx_type&,
double*,
106 const octave_idx_type&, octave_idx_type&
111 const octave_idx_type&,
const octave_idx_type&,
112 double*,
const octave_idx_type&,
113 const double&,
double&,
double*,
114 octave_idx_type*, octave_idx_type&
117 F77_FUNC (dptsv, DPTSV) (
const octave_idx_type&,
const octave_idx_type&,
118 double*,
double*,
double*,
const octave_idx_type&,
122 F77_FUNC (dgtsv, DGTSV) (
const octave_idx_type&,
const octave_idx_type&,
123 double*,
double*,
double*,
double*,
124 const octave_idx_type&, octave_idx_type&);
127 F77_FUNC (dgttrf, DGTTRF) (
const octave_idx_type&,
double*,
double*,
128 double*,
double*, octave_idx_type*,
133 const octave_idx_type&,
const octave_idx_type&,
134 const double*,
const double*,
const double*,
135 const double*,
const octave_idx_type*,
136 double *,
const octave_idx_type&, octave_idx_type&
140 F77_FUNC (zptsv, ZPTSV) (
const octave_idx_type&,
const octave_idx_type&,
141 double*,
Complex*, Complex*,
const octave_idx_type&,
145 F77_FUNC (zgtsv, ZGTSV) (
const octave_idx_type&,
const octave_idx_type&,
146 Complex*, Complex*, Complex*, Complex*,
147 const octave_idx_type&, octave_idx_type&);
196 if (nr != nr_a || nc != nc_a || nz != nz_a)
213 return !(*
this == a);
222 if (nr == nc && nr > 0)
275 return max (dummy_idx, dim);
299 if (nr == 0 || nc == 0 || dim >= dv.
length ())
309 if (
ridx (i) != idx_j)
320 double tmp =
data (i);
324 else if (
xisnan (tmp_max) || tmp > tmp_max)
332 idx_arg.
elem (j) =
xisnan (tmp_max) ? 0 : idx_j;
340 result.
xcidx (0) = 0;
343 double tmp =
elem (idx_arg(j), j);
346 result.
xdata (ii) = tmp;
347 result.
xridx (ii++) = 0;
349 result.
xcidx (j+1) = ii;
357 if (nr == 0 || nc == 0 || dim >= dv.
length ())
367 if (found[
ridx (i)] == -j)
368 found[
ridx (i)] = -j - 1;
371 if (found[i] > -nc && found[i] < 0)
372 idx_arg.
elem (i) = -found[i];
380 double tmp =
data (i);
384 else if (ix == -1 || tmp >
elem (ir, ix))
385 idx_arg.
elem (ir) = j;
391 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
397 result.
xcidx (0) = 0;
398 result.
xcidx (1) = nel;
401 if (idx_arg(j) == -1)
405 result.
xridx (ii++) = j;
409 double tmp =
elem (j, idx_arg(j));
412 result.
xdata (ii) = tmp;
413 result.
xridx (ii++) = j;
426 return min (dummy_idx, dim);
450 if (nr == 0 || nc == 0 || dim >= dv.
length ())
460 if (
ridx (i) != idx_j)
471 double tmp =
data (i);
475 else if (
xisnan (tmp_min) || tmp < tmp_min)
483 idx_arg.
elem (j) =
xisnan (tmp_min) ? 0 : idx_j;
491 result.
xcidx (0) = 0;
494 double tmp =
elem (idx_arg(j), j);
497 result.
xdata (ii) = tmp;
498 result.
xridx (ii++) = 0;
500 result.
xcidx (j+1) = ii;
508 if (nr == 0 || nc == 0 || dim >= dv.
length ())
518 if (found[
ridx (i)] == -j)
519 found[
ridx (i)] = -j - 1;
522 if (found[i] > -nc && found[i] < 0)
523 idx_arg.
elem (i) = -found[i];
531 double tmp =
data (i);
535 else if (ix == -1 || tmp <
elem (ir, ix))
536 idx_arg.
elem (ir) = j;
542 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
548 result.
xcidx (0) = 0;
549 result.
xcidx (1) = nel;
552 if (idx_arg(j) == -1)
556 result.
xridx (ii++) = j;
560 double tmp =
elem (j, idx_arg(j));
563 result.
xdata (ii) = tmp;
564 result.
xridx (ii++) = j;
599 retval(j) =
data (k);
624 if (rb.
rows () > 0 && rb.
cols () > 0)
634 if (rb.
rows () > 0 && rb.
cols () > 0)
720 retval.
xcidx (0) = 0;
728 retval.
xdata (ii) = tmp;
732 retval.
xcidx (i+1) = ii;
764 if (x_nr != y_nr || x_nc != y_nc)
776 bool ja_lt_max= ja < ja_max;
780 bool jb_lt_max = jb < jb_max;
782 while (ja_lt_max || jb_lt_max)
786 || (ja_lt_max && (x.
ridx (ja) < y.
ridx (jb))))
792 ja_lt_max= ja < ja_max;
794 else if ((! ja_lt_max)
795 || (jb_lt_max && (y.
ridx (jb) < x.
ridx (ja))))
798 jb_lt_max= jb < jb_max;
810 ja_lt_max= ja < ja_max;
812 jb_lt_max= jb < jb_max;
833 return inverse (mattype, info, rcond, 0, 0);
841 return inverse (mattype, info, rcond, 0, 0);
848 return inverse (mattype, info, rcond, 0, 0);
853 double& rcond,
const bool,
854 const bool calccond)
const
862 if (nr == 0 || nc == 0 || nr != nc)
863 (*current_liboctave_error_handler) (
"inverse requires square matrix");
867 int typ = mattyp.
type ();
878 double *v = retval.
data ();
886 double tmp = fabs (v[i]);
907 double& rcond,
const bool,
908 const bool calccond)
const
916 if (nr == 0 || nc == 0 || nr != nc)
917 (*current_liboctave_error_handler) (
"inverse requires square matrix");
921 int typ = mattyp.
type ();
928 double ainvnorm = 0.;
937 atmp += fabs (
data (i));
962 retval.
xcidx (i) = cx;
963 retval.
xridx (cx) = i;
964 retval.
xdata (cx) = 1.0;
978 (*current_liboctave_error_handler)
979 (
"division by zero");
980 goto inverse_singular;
986 rpX = retval.
xridx (colXp);
995 v -= retval.
xdata (colXp) *
data (colUp);
1000 while (rpX < j && rpU < j && colXp < cx && colUp < nz);
1004 colUp =
cidx (j+1) - 1;
1007 double pivot =
data (colUp);
1008 if (pivot == 0. ||
ridx (colUp) != j)
1010 (*current_liboctave_error_handler)
1011 (
"division by zero");
1012 goto inverse_singular;
1023 retval.
xridx (cx) = j;
1024 retval.
xdata (cx) = v / pivot;
1032 colUp =
cidx (i+1) - 1;
1035 double pivot =
data (colUp);
1036 if (pivot == 0. ||
ridx (colUp) != i)
1038 (*current_liboctave_error_handler) (
"division by zero");
1039 goto inverse_singular;
1044 retval.
xdata (j) /= pivot;
1046 retval.
xcidx (nr) = cx;
1091 k <
cidx (jidx+1); k++)
1100 pivot =
data (cidx (jidx+1) - 1);
1102 pivot =
data (cidx (jidx));
1105 (*current_liboctave_error_handler)
1106 (
"division by zero");
1107 goto inverse_singular;
1110 work[j] = v / pivot;
1116 colUp =
cidx (perm[iidx]+1) - 1;
1118 colUp =
cidx (perm[iidx]);
1120 double pivot =
data (colUp);
1123 (*current_liboctave_error_handler)
1124 (
"division by zero");
1125 goto inverse_singular;
1139 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
1143 retval.
xcidx (i) = cx;
1147 retval.
xridx (cx) = j;
1148 retval.
xdata (cx++) = work[j];
1152 retval.
xcidx (nr) = cx;
1163 i < retval.
cidx (j+1); i++)
1164 atmp += fabs (retval.
data (i));
1165 if (atmp > ainvnorm)
1169 rcond = 1. / ainvnorm / anorm;
1184 double& rcond,
int,
int calc_cond)
const
1186 int typ = mattype.
type (
false);
1190 typ = mattype.
type (*
this);
1193 ret =
dinverse (mattype, info, rcond,
true, calc_cond);
1207 rcond = fact.
rcond ();
1213 info, rcond2,
true,
false);
1232 rcond = fact.
rcond ();
1235 info, rcond2,
true,
false);
1269 if (nr == 0 || nc == 0 || nr != nc)
1278 Matrix Control (UMFPACK_CONTROL, 1);
1284 Control (UMFPACK_PRL) = tmp;
1289 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
1290 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
1296 Control (UMFPACK_FIXQ) = tmp;
1299 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1305 const double *Ax =
data ();
1307 UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1, control);
1310 Matrix Info (1, UMFPACK_INFO);
1313 Ax, 0, &Symbolic, control, info);
1317 (*current_liboctave_error_handler)
1318 (
"SparseMatrix::determinant symbolic factorization failed");
1331 &Numeric, control, info);
1334 rcond = Info (UMFPACK_RCOND);
1338 (*current_liboctave_error_handler)
1339 (
"SparseMatrix::determinant numeric factorization failed");
1352 status =
UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric,
1357 (*current_liboctave_error_handler)
1358 (
"SparseMatrix::determinant error calculating determinant");
1364 retval =
DET (c10, e10, 10);
1371 (*current_liboctave_error_handler) (
"UMFPACK not installed");
1380 double& rcond, solve_singularity_handler,
1381 bool calc_cond)
const
1390 if (nr != b.
rows ())
1392 (
"matrix dimension mismatch solution of linear equations");
1393 else if (nr == 0 || nc == 0 || b.
cols () == 0)
1398 int typ = mattype.
type ();
1407 retval(i,j) = b(i,j) /
data (i);
1412 retval(k,j) = b(
ridx (i),j) /
data (i);
1420 double tmp = fabs (
data (i));
1426 rcond = dmin / dmax;
1441 solve_singularity_handler,
bool calc_cond)
const
1450 if (nr != b.
rows ())
1452 (
"matrix dimension mismatch solution of linear equations");
1453 else if (nr == 0 || nc == 0 || b.
cols () == 0)
1458 int typ = mattype.
type ();
1467 retval.
xcidx (0) = 0;
1474 if (b.
ridx (i) >= nm)
1479 retval.
xcidx (j+1) = ii;
1489 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1497 retval.
xridx (ii) = l;
1501 retval.
xcidx (j+1) = ii;
1510 double tmp = fabs (
data (i));
1516 rcond = dmin / dmax;
1531 solve_singularity_handler,
bool calc_cond)
const
1540 if (nr != b.
rows ())
1542 (
"matrix dimension mismatch solution of linear equations");
1543 else if (nr == 0 || nc == 0 || b.
cols () == 0)
1548 int typ = mattype.
type ();
1557 retval(i,j) = b(i,j) /
data (i);
1562 retval(k,j) = b(
ridx (i),j) /
data (i);
1570 double tmp = fabs (
data (i));
1576 rcond = dmin / dmax;
1591 solve_singularity_handler,
bool calc_cond)
const
1600 if (nr != b.
rows ())
1602 (
"matrix dimension mismatch solution of linear equations");
1603 else if (nr == 0 || nc == 0 || b.
cols () == 0)
1608 int typ = mattype.
type ();
1617 retval.
xcidx (0) = 0;
1624 if (b.
ridx (i) >= nm)
1629 retval.
xcidx (j+1) = ii;
1639 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1647 retval.
xridx (ii) = l;
1651 retval.
xcidx (j+1) = ii;
1660 double tmp = fabs (
data (i));
1666 rcond = dmin / dmax;
1681 solve_singularity_handler sing_handler,
1682 bool calc_cond)
const
1691 if (nr != b.
rows ())
1693 (
"matrix dimension mismatch solution of linear equations");
1694 else if (nr == 0 || nc == 0 || b.
cols () == 0)
1699 int typ = mattype.
type ();
1705 double ainvnorm = 0.;
1716 atmp += fabs (
data (i));
1724 retval.
resize (nc, b_nc);
1745 goto triangular_error;
1748 double tmp = work[k] /
data (
cidx (kidx+1)-1);
1751 i <
cidx (kidx+1)-1; i++)
1754 work[iidx] = work[iidx] - tmp *
data (i);
1760 retval.
xelem (perm[i], j) = work[i];
1779 double tmp = work[k] /
data (
cidx (iidx+1)-1);
1782 i <
cidx (iidx+1)-1; i++)
1785 work[idx2] = work[idx2] - tmp *
data (i);
1792 atmp += fabs (work[i]);
1795 if (atmp > ainvnorm)
1798 rcond = 1. / ainvnorm / anorm;
1804 retval.
resize (nc, b_nc);
1821 goto triangular_error;
1824 double tmp = work[k] /
data (
cidx (k+1)-1);
1829 work[iidx] = work[iidx] - tmp *
data (i);
1835 retval.
xelem (i, j) = work[i];
1852 double tmp = work[k] /
data (
cidx (k+1)-1);
1857 work[iidx] = work[iidx] - tmp *
data (i);
1864 atmp += fabs (work[i]);
1867 if (atmp > ainvnorm)
1870 rcond = 1. / ainvnorm / anorm;
1879 sing_handler (rcond);
1886 volatile double rcond_plus_one = rcond + 1.0;
1888 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
1894 sing_handler (rcond);
1911 solve_singularity_handler sing_handler,
1912 bool calc_cond)
const
1921 if (nr != b.
rows ())
1923 (
"matrix dimension mismatch solution of linear equations");
1924 else if (nr == 0 || nc == 0 || b.
cols () == 0)
1929 int typ = mattype.
type ();
1935 double ainvnorm = 0.;
1945 atmp += fabs (
data (i));
1954 retval.
xcidx (0) = 0;
1984 goto triangular_error;
1987 double tmp = work[k] /
data (
cidx (kidx+1)-1);
1990 i <
cidx (kidx+1)-1; i++)
1993 work[iidx] = work[iidx] - tmp *
data (i);
2005 if (ii + new_nnz > x_nz)
2014 if (work[rperm[i]] != 0.)
2016 retval.
xridx (ii) = i;
2017 retval.
xdata (ii++) = work[rperm[i]];
2019 retval.
xcidx (j+1) = ii;
2040 double tmp = work[k] /
data (
cidx (iidx+1)-1);
2043 i <
cidx (iidx+1)-1; i++)
2046 work[idx2] = work[idx2] - tmp *
data (i);
2053 atmp += fabs (work[i]);
2056 if (atmp > ainvnorm)
2059 rcond = 1. / ainvnorm / anorm;
2081 goto triangular_error;
2084 double tmp = work[k] /
data (
cidx (k+1)-1);
2089 work[iidx] = work[iidx] - tmp *
data (i);
2101 if (ii + new_nnz > x_nz)
2112 retval.
xridx (ii) = i;
2113 retval.
xdata (ii++) = work[i];
2115 retval.
xcidx (j+1) = ii;
2134 double tmp = work[k] /
data (
cidx (k+1)-1);
2137 i <
cidx (k+1)-1; i++)
2140 work[iidx] = work[iidx] - tmp *
data (i);
2147 atmp += fabs (work[i]);
2150 if (atmp > ainvnorm)
2153 rcond = 1. / ainvnorm / anorm;
2162 sing_handler (rcond);
2169 volatile double rcond_plus_one = rcond + 1.0;
2171 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
2177 sing_handler (rcond);
2193 solve_singularity_handler sing_handler,
2194 bool calc_cond)
const
2203 if (nr != b.
rows ())
2205 (
"matrix dimension mismatch solution of linear equations");
2206 else if (nr == 0 || nc == 0 || b.
cols () == 0)
2211 int typ = mattype.
type ();
2217 double ainvnorm = 0.;
2228 atmp += fabs (
data (i));
2236 retval.
resize (nc, b_nc);
2257 goto triangular_error;
2263 i <
cidx (kidx+1)-1; i++)
2266 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2272 retval.
xelem (perm[i], j) = cwork[i];
2292 double tmp = work[k] /
data (
cidx (iidx+1)-1);
2295 i <
cidx (iidx+1)-1; i++)
2298 work[idx2] = work[idx2] - tmp *
data (i);
2305 atmp += fabs (work[i]);
2308 if (atmp > ainvnorm)
2311 rcond = 1. / ainvnorm / anorm;
2317 retval.
resize (nc, b_nc);
2334 goto triangular_error;
2342 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2348 retval.
xelem (i, j) = cwork[i];
2366 double tmp = work[k] /
data (
cidx (k+1)-1);
2369 i <
cidx (k+1)-1; i++)
2372 work[iidx] = work[iidx] - tmp *
data (i);
2379 atmp += fabs (work[i]);
2382 if (atmp > ainvnorm)
2385 rcond = 1. / ainvnorm / anorm;
2394 sing_handler (rcond);
2401 volatile double rcond_plus_one = rcond + 1.0;
2403 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
2409 sing_handler (rcond);
2426 solve_singularity_handler sing_handler,
2427 bool calc_cond)
const
2436 if (nr != b.
rows ())
2438 (
"matrix dimension mismatch solution of linear equations");
2439 else if (nr == 0 || nc == 0 || b.
cols () == 0)
2444 int typ = mattype.
type ();
2450 double ainvnorm = 0.;
2460 atmp += fabs (
data (i));
2469 retval.
xcidx (0) = 0;
2499 goto triangular_error;
2505 i <
cidx (kidx+1)-1; i++)
2508 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2520 if (ii + new_nnz > x_nz)
2529 if (cwork[rperm[i]] != 0.)
2531 retval.
xridx (ii) = i;
2532 retval.
xdata (ii++) = cwork[rperm[i]];
2534 retval.
xcidx (j+1) = ii;
2556 double tmp = work[k] /
data (
cidx (iidx+1)-1);
2559 i <
cidx (iidx+1)-1; i++)
2562 work[idx2] = work[idx2] - tmp *
data (i);
2569 atmp += fabs (work[i]);
2572 if (atmp > ainvnorm)
2575 rcond = 1. / ainvnorm / anorm;
2597 goto triangular_error;
2605 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2617 if (ii + new_nnz > x_nz)
2628 retval.
xridx (ii) = i;
2629 retval.
xdata (ii++) = cwork[i];
2631 retval.
xcidx (j+1) = ii;
2651 double tmp = work[k] /
data (
cidx (k+1)-1);
2654 i <
cidx (k+1)-1; i++)
2657 work[iidx] = work[iidx] - tmp *
data (i);
2664 atmp += fabs (work[i]);
2667 if (atmp > ainvnorm)
2670 rcond = 1. / ainvnorm / anorm;
2679 sing_handler (rcond);
2686 volatile double rcond_plus_one = rcond + 1.0;
2688 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
2694 sing_handler (rcond);
2711 solve_singularity_handler sing_handler,
2712 bool calc_cond)
const
2721 if (nr != b.
rows ())
2723 (
"matrix dimension mismatch solution of linear equations");
2724 else if (nr == 0 || nc == 0 || b.
cols () == 0)
2729 int typ = mattype.
type ();
2735 double ainvnorm = 0.;
2746 atmp += fabs (
data (i));
2754 retval.
resize (nc, b_nc);
2764 work[perm[i]] = b(i,j);
2774 if (perm[
ridx (i)] < minr)
2776 minr = perm[
ridx (i)];
2780 if (minr != k ||
data (mini) == 0)
2783 goto triangular_error;
2786 double tmp = work[k] /
data (mini);
2794 work[iidx] = work[iidx] - tmp *
data (i);
2800 retval(i, j) = work[i];
2821 i <
cidx (k+1); i++)
2822 if (perm[
ridx (i)] < minr)
2824 minr = perm[
ridx (i)];
2828 double tmp = work[k] /
data (mini);
2831 i <
cidx (k+1); i++)
2837 work[iidx] = work[iidx] - tmp *
data (i);
2845 atmp += fabs (work[i]);
2848 if (atmp > ainvnorm)
2851 rcond = 1. / ainvnorm / anorm;
2857 retval.
resize (nc, b_nc, 0.);
2872 goto triangular_error;
2875 double tmp = work[k] /
data (
cidx (k));
2878 i <
cidx (k+1); i++)
2881 work[iidx] = work[iidx] - tmp *
data (i);
2887 retval.
xelem (i, j) = work[i];
2905 double tmp = work[k] /
data (
cidx (k));
2908 i <
cidx (k+1); i++)
2911 work[iidx] = work[iidx] - tmp *
data (i);
2918 atmp += fabs (work[i]);
2921 if (atmp > ainvnorm)
2924 rcond = 1. / ainvnorm / anorm;
2933 sing_handler (rcond);
2940 volatile double rcond_plus_one = rcond + 1.0;
2942 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
2948 sing_handler (rcond);
2965 solve_singularity_handler sing_handler,
2966 bool calc_cond)
const
2975 if (nr != b.
rows ())
2977 (
"matrix dimension mismatch solution of linear equations");
2978 else if (nr == 0 || nc == 0 || b.
cols () == 0)
2983 int typ = mattype.
type ();
2989 double ainvnorm = 0.;
2999 atmp += fabs (
data (i));
3008 retval.
xcidx (0) = 0;
3022 work[perm[b.
ridx (i)]] = b.
data (i);
3032 if (perm[
ridx (i)] < minr)
3034 minr = perm[
ridx (i)];
3038 if (minr != k ||
data (mini) == 0)
3041 goto triangular_error;
3044 double tmp = work[k] /
data (mini);
3052 work[iidx] = work[iidx] - tmp *
data (i);
3064 if (ii + new_nnz > x_nz)
3075 retval.
xridx (ii) = i;
3076 retval.
xdata (ii++) = work[i];
3078 retval.
xcidx (j+1) = ii;
3101 i <
cidx (k+1); i++)
3102 if (perm[
ridx (i)] < minr)
3104 minr = perm[
ridx (i)];
3108 double tmp = work[k] /
data (mini);
3111 i <
cidx (k+1); i++)
3117 work[iidx] = work[iidx] - tmp *
data (i);
3125 atmp += fabs (work[i]);
3128 if (atmp > ainvnorm)
3131 rcond = 1. / ainvnorm / anorm;
3152 goto triangular_error;
3155 double tmp = work[k] /
data (
cidx (k));
3160 work[iidx] = work[iidx] - tmp *
data (i);
3172 if (ii + new_nnz > x_nz)
3183 retval.
xridx (ii) = i;
3184 retval.
xdata (ii++) = work[i];
3186 retval.
xcidx (j+1) = ii;
3206 double tmp = work[k] /
data (
cidx (k));
3209 i <
cidx (k+1); i++)
3212 work[iidx] = work[iidx] - tmp *
data (i);
3219 atmp += fabs (work[i]);
3222 if (atmp > ainvnorm)
3225 rcond = 1. / ainvnorm / anorm;
3234 sing_handler (rcond);
3241 volatile double rcond_plus_one = rcond + 1.0;
3243 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
3249 sing_handler (rcond);
3266 solve_singularity_handler sing_handler,
3267 bool calc_cond)
const
3276 if (nr != b.
rows ())
3278 (
"matrix dimension mismatch solution of linear equations");
3279 else if (nr == 0 || nc == 0 || b.
cols () == 0)
3284 int typ = mattype.
type ();
3290 double ainvnorm = 0.;
3301 atmp += fabs (
data (i));
3309 retval.
resize (nc, b_nc);
3318 cwork[perm[i]] = b(i,j);
3328 if (perm[
ridx (i)] < minr)
3330 minr = perm[
ridx (i)];
3334 if (minr != k ||
data (mini) == 0)
3337 goto triangular_error;
3348 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3354 retval(i, j) = cwork[i];
3376 i <
cidx (k+1); i++)
3377 if (perm[
ridx (i)] < minr)
3379 minr = perm[
ridx (i)];
3383 double tmp = work[k] /
data (mini);
3386 i <
cidx (k+1); i++)
3392 work[iidx] = work[iidx] - tmp *
data (i);
3400 atmp += fabs (work[i]);
3403 if (atmp > ainvnorm)
3406 rcond = 1. / ainvnorm / anorm;
3412 retval.
resize (nc, b_nc, 0.);
3428 goto triangular_error;
3436 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3442 retval.
xelem (i, j) = cwork[i];
3461 double tmp = work[k] /
data (
cidx (k));
3464 i <
cidx (k+1); i++)
3467 work[iidx] = work[iidx] - tmp *
data (i);
3474 atmp += fabs (work[i]);
3477 if (atmp > ainvnorm)
3480 rcond = 1. / ainvnorm / anorm;
3489 sing_handler (rcond);
3496 volatile double rcond_plus_one = rcond + 1.0;
3498 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
3504 sing_handler (rcond);
3521 solve_singularity_handler sing_handler,
3522 bool calc_cond)
const
3531 if (nr != b.
rows ())
3533 (
"matrix dimension mismatch solution of linear equations");
3534 else if (nr == 0 || nc == 0 || b.
cols () == 0)
3539 int typ = mattype.
type ();
3545 double ainvnorm = 0.;
3555 atmp += fabs (
data (i));
3564 retval.
xcidx (0) = 0;
3578 cwork[perm[b.
ridx (i)]] = b.
data (i);
3588 if (perm[
ridx (i)] < minr)
3590 minr = perm[
ridx (i)];
3594 if (minr != k ||
data (mini) == 0)
3597 goto triangular_error;
3608 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3620 if (ii + new_nnz > x_nz)
3631 retval.
xridx (ii) = i;
3632 retval.
xdata (ii++) = cwork[i];
3634 retval.
xcidx (j+1) = ii;
3658 i <
cidx (k+1); i++)
3659 if (perm[
ridx (i)] < minr)
3661 minr = perm[
ridx (i)];
3665 double tmp = work[k] /
data (mini);
3668 i <
cidx (k+1); i++)
3674 work[iidx] = work[iidx] - tmp *
data (i);
3682 atmp += fabs (work[i]);
3685 if (atmp > ainvnorm)
3688 rcond = 1. / ainvnorm / anorm;
3709 goto triangular_error;
3717 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3729 if (ii + new_nnz > x_nz)
3740 retval.
xridx (ii) = i;
3741 retval.
xdata (ii++) = cwork[i];
3743 retval.
xcidx (j+1) = ii;
3764 double tmp = work[k] /
data (
cidx (k));
3767 i <
cidx (k+1); i++)
3770 work[iidx] = work[iidx] - tmp *
data (i);
3777 atmp += fabs (work[i]);
3780 if (atmp > ainvnorm)
3783 rcond = 1. / ainvnorm / anorm;
3792 sing_handler (rcond);
3799 volatile double rcond_plus_one = rcond + 1.0;
3801 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
3807 sing_handler (rcond);
3824 solve_singularity_handler sing_handler,
3825 bool calc_cond)
const
3833 if (nr != nc || nr != b.
rows ())
3835 (
"matrix dimension mismatch solution of linear equations");
3836 else if (nr == 0 || b.
cols () == 0)
3839 (*current_liboctave_error_handler)
3840 (
"calculation of condition number not implemented");
3844 volatile int typ = mattype.
type ();
3862 D[nc-1] =
data (ii);
3878 else if (
ridx (i) == j + 1)
3887 F77_XFCN (dptsv, DPTSV, (nr, b_nc, D, DL, result,
3913 DL[j] =
data (ii++);
3914 DU[j] =
data (ii++);
3916 D[nc-1] =
data (ii);
3933 else if (
ridx (i) == j + 1)
3935 else if (
ridx (i) == j - 1)
3944 F77_XFCN (dgtsv, DGTSV, (nr, b_nc, DL, D, DU, result,
3954 sing_handler (rcond);
3964 (*current_liboctave_error_handler) (
"incorrect matrix type");
3973 solve_singularity_handler sing_handler,
3974 bool calc_cond)
const
3982 if (nr != nc || nr != b.
rows ())
3984 (
"matrix dimension mismatch solution of linear equations");
3985 else if (nr == 0 || b.
cols () == 0)
3988 (*current_liboctave_error_handler)
3989 (
"calculation of condition number not implemented");
3993 int typ = mattype.
type ();
4014 DL[j] =
data (ii++);
4015 DU[j] =
data (ii++);
4017 D[nc-1] =
data (ii);
4034 else if (
ridx (i) == j + 1)
4036 else if (
ridx (i) == j - 1)
4041 F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
4050 sing_handler (rcond);
4063 retval.
xcidx (0) = 0;
4077 nr, 1, DL, D, DU, DU2, pipvt,
4078 work, b.
rows (), err
4088 if (ii + new_nnz > x_nz)
4099 retval.
xridx (ii) = i;
4100 retval.
xdata (ii++) = work[i];
4102 retval.
xcidx (j+1) = ii;
4109 (*current_liboctave_error_handler) (
"incorrect matrix type");
4118 solve_singularity_handler sing_handler,
4119 bool calc_cond)
const
4127 if (nr != nc || nr != b.
rows ())
4129 (
"matrix dimension mismatch solution of linear equations");
4130 else if (nr == 0 || b.
cols () == 0)
4133 (*current_liboctave_error_handler)
4134 (
"calculation of condition number not implemented");
4138 volatile int typ = mattype.
type ();
4156 D[nc-1] =
data (ii);
4172 else if (
ridx (i) == j + 1)
4184 F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result,
4208 DL[j] =
data (ii++);
4209 DU[j] =
data (ii++);
4211 D[nc-1] =
data (ii);
4228 else if (
ridx (i) == j + 1)
4230 else if (
ridx (i) == j - 1)
4242 F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result,
4252 sing_handler (rcond);
4260 (*current_liboctave_error_handler) (
"incorrect matrix type");
4269 solve_singularity_handler sing_handler,
4270 bool calc_cond)
const
4278 if (nr != nc || nr != b.
rows ())
4280 (
"matrix dimension mismatch solution of linear equations");
4281 else if (nr == 0 || b.
cols () == 0)
4284 (*current_liboctave_error_handler)
4285 (
"calculation of condition number not implemented");
4289 int typ = mattype.
type ();
4310 DL[j] =
data (ii++);
4311 DU[j] =
data (ii++);
4313 D[nc-1] =
data (ii);
4330 else if (
ridx (i) == j + 1)
4332 else if (
ridx (i) == j - 1)
4337 F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
4346 sing_handler (rcond);
4367 retval.
xcidx (0) = 0;
4380 nr, 1, DL, D, DU, DU2, pipvt,
4386 (*current_liboctave_error_handler)
4387 (
"SparseMatrix::solve solve failed");
4395 nr, 1, DL, D, DU, DU2, pipvt,
4401 (*current_liboctave_error_handler)
4402 (
"SparseMatrix::solve solve failed");
4412 if (Bx[i] != 0. || Bz[i] != 0.)
4415 if (ii + new_nnz > x_nz)
4424 if (Bx[i] != 0. || Bz[i] != 0.)
4426 retval.
xridx (ii) = i;
4427 retval.
xdata (ii++) =
4431 retval.
xcidx (j+1) = ii;
4438 (*current_liboctave_error_handler) (
"incorrect matrix type");
4447 solve_singularity_handler sing_handler,
4448 bool calc_cond)
const
4456 if (nr != nc || nr != b.
rows ())
4458 (
"matrix dimension mismatch solution of linear equations");
4459 else if (nr == 0 || b.
cols () == 0)
4464 volatile int typ = mattype.
type ();
4480 tmp_data[ii++] = 0.;
4488 m_band(ri - j, j) =
data (i);
4498 nr, n_lower, tmp_data, ldm, err
4521 nr, n_lower, tmp_data, ldm,
4522 anorm, rcond, pz, piz, err
4528 volatile double rcond_plus_one = rcond + 1.0;
4530 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
4536 sing_handler (rcond);
4555 nr, n_lower, b_nc, tmp_data,
4556 ldm, result, b.
rows (), err
4561 (*current_liboctave_error_handler)
4562 (
"SparseMatrix::solve solve failed");
4585 tmp_data[ii++] = 0.;
4590 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
4600 atmp += fabs (
data (i));
4609 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
4621 sing_handler (rcond);
4639 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4640 anorm, rcond, pz, piz, err
4646 volatile double rcond_plus_one = rcond + 1.0;
4648 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
4654 sing_handler (rcond);
4674 nr, n_lower, n_upper, b_nc, tmp_data,
4675 ldm, pipvt, result, b.
rows (), err
4681 (*current_liboctave_error_handler) (
"incorrect matrix type");
4690 solve_singularity_handler sing_handler,
4691 bool calc_cond)
const
4699 if (nr != nc || nr != b.
rows ())
4701 (
"matrix dimension mismatch solution of linear equations");
4702 else if (nr == 0 || b.
cols () == 0)
4707 volatile int typ = mattype.
type ();
4724 tmp_data[ii++] = 0.;
4732 m_band(ri - j, j) =
data (i);
4742 nr, n_lower, tmp_data, ldm, err
4763 nr, n_lower, tmp_data, ldm,
4764 anorm, rcond, pz, piz, err
4770 volatile double rcond_plus_one = rcond + 1.0;
4772 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
4778 sing_handler (rcond);
4800 retval.
xcidx (0) = 0;
4804 Bx[i] = b.
elem (i, j);
4808 nr, n_lower, 1, tmp_data,
4814 (*current_liboctave_error_handler)
4815 (
"SparseMatrix::solve solve failed");
4830 sz = (sz > 10 ? sz : 10) + x_nz;
4834 retval.
xdata (ii) = tmp;
4835 retval.
xridx (ii++) = i;
4838 retval.
xcidx (j+1) = ii;
4862 tmp_data[ii++] = 0.;
4867 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
4877 atmp += fabs (
data (i));
4886 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
4896 sing_handler (rcond);
4914 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4915 anorm, rcond, pz, piz, err
4921 volatile double rcond_plus_one = rcond + 1.0;
4923 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
4929 sing_handler (rcond);
4945 retval.
xcidx (0) = 0;
4955 i < b.
cidx (j+1); i++)
4960 nr, n_lower, n_upper, 1, tmp_data,
4961 ldm, pipvt, work, b.
rows (), err
4971 if (ii + new_nnz > x_nz)
4982 retval.
xridx (ii) = i;
4983 retval.
xdata (ii++) = work[i];
4985 retval.
xcidx (j+1) = ii;
4993 (*current_liboctave_error_handler) (
"incorrect matrix type");
5002 solve_singularity_handler sing_handler,
5003 bool calc_cond)
const
5011 if (nr != nc || nr != b.
rows ())
5013 (
"matrix dimension mismatch solution of linear equations");
5014 else if (nr == 0 || b.
cols () == 0)
5019 volatile int typ = mattype.
type ();
5036 tmp_data[ii++] = 0.;
5044 m_band(ri - j, j) =
data (i);
5054 nr, n_lower, tmp_data, ldm, err
5077 nr, n_lower, tmp_data, ldm,
5078 anorm, rcond, pz, piz, err
5084 volatile double rcond_plus_one = rcond + 1.0;
5086 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
5092 sing_handler (rcond);
5110 retval.
resize (b_nr, b_nc);
5123 nr, n_lower, 1, tmp_data,
5129 (*current_liboctave_error_handler)
5130 (
"SparseMatrix::solve solve failed");
5137 nr, n_lower, 1, tmp_data,
5138 ldm, Bz, b.
rows (), err
5143 (*current_liboctave_error_handler)
5144 (
"SparseMatrix::solve solve failed");
5150 retval(i, j) =
Complex (Bx[i], Bz[i]);
5172 tmp_data[ii++] = 0.;
5177 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
5187 atmp += fabs (
data (i));
5196 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
5206 sing_handler (rcond);
5224 nr, n_lower, tmp_data, ldm,
5225 anorm, rcond, pz, piz, err
5231 volatile double rcond_plus_one = rcond + 1.0;
5233 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
5239 sing_handler (rcond);
5269 nr, n_lower, n_upper, 1, tmp_data,
5270 ldm, pipvt, Bx, b.
rows (), err
5275 nr, n_lower, n_upper, 1, tmp_data,
5276 ldm, pipvt, Bz, b.
rows (), err
5280 retval(i, j) =
Complex (Bx[i], Bz[i]);
5286 (*current_liboctave_error_handler) (
"incorrect matrix type");
5295 solve_singularity_handler sing_handler,
5296 bool calc_cond)
const
5304 if (nr != nc || nr != b.
rows ())
5306 (
"matrix dimension mismatch solution of linear equations");
5307 else if (nr == 0 || b.
cols () == 0)
5312 volatile int typ = mattype.
type ();
5329 tmp_data[ii++] = 0.;
5337 m_band(ri - j, j) =
data (i);
5347 nr, n_lower, tmp_data, ldm, err
5371 nr, n_lower, tmp_data, ldm,
5372 anorm, rcond, pz, piz, err
5378 volatile double rcond_plus_one = rcond + 1.0;
5380 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
5386 sing_handler (rcond);
5409 retval.
xcidx (0) = 0;
5422 nr, n_lower, 1, tmp_data,
5428 (*current_liboctave_error_handler)
5429 (
"SparseMatrix::solve solve failed");
5436 nr, n_lower, 1, tmp_data,
5442 (*current_liboctave_error_handler)
5443 (
"SparseMatrix::solve solve failed");
5453 if (Bx[i] != 0. || Bz[i] != 0.)
5456 if (ii + new_nnz > x_nz)
5465 if (Bx[i] != 0. || Bz[i] != 0.)
5467 retval.
xridx (ii) = i;
5468 retval.
xdata (ii++) =
5472 retval.
xcidx (j+1) = ii;
5496 tmp_data[ii++] = 0.;
5501 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
5511 atmp += fabs (
data (i));
5520 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
5530 sing_handler (rcond);
5548 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5549 anorm, rcond, pz, piz, err
5555 volatile double rcond_plus_one = rcond + 1.0;
5557 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
5563 sing_handler (rcond);
5579 retval.
xcidx (0) = 0;
5593 i < b.
cidx (j+1); i++)
5602 nr, n_lower, n_upper, 1, tmp_data,
5603 ldm, pipvt, Bx, b.
rows (), err
5608 nr, n_lower, n_upper, 1, tmp_data,
5609 ldm, pipvt, Bz, b.
rows (), err
5616 if (Bx[i] != 0. || Bz[i] != 0.)
5619 if (ii + new_nnz > x_nz)
5628 if (Bx[i] != 0. || Bz[i] != 0.)
5630 retval.
xridx (ii) = i;
5631 retval.
xdata (ii++) =
5634 retval.
xcidx (j+1) = ii;
5642 (*current_liboctave_error_handler) (
"incorrect matrix type");
5650 Matrix &Info, solve_singularity_handler sing_handler,
5651 bool calc_cond)
const
5659 Control =
Matrix (UMFPACK_CONTROL, 1);
5665 Control (UMFPACK_PRL) = tmp;
5669 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
5670 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
5676 Control (UMFPACK_FIXQ) = tmp;
5682 const double *Ax =
data ();
5686 UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1, control);
5689 Info =
Matrix (1, UMFPACK_INFO);
5691 int status =
UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, 0,
5692 &Symbolic, control, info);
5696 (*current_liboctave_error_handler)
5697 (
"SparseMatrix::solve symbolic factorization failed");
5710 &Numeric, control, info);
5714 rcond = Info (UMFPACK_RCOND);
5717 volatile double rcond_plus_one = rcond + 1.0;
5719 if (status == UMFPACK_WARNING_singular_matrix
5720 || rcond_plus_one == 1.0 ||
xisnan (rcond))
5727 sing_handler (rcond);
5731 else if (status < 0)
5733 (*current_liboctave_error_handler)
5734 (
"SparseMatrix::solve numeric factorization failed");
5751 (*current_liboctave_error_handler) (
"UMFPACK not installed");
5760 solve_singularity_handler sing_handler,
5761 bool calc_cond)
const
5769 if (nr != nc || nr != b.
rows ())
5771 (
"matrix dimension mismatch solution of linear equations");
5772 else if (nr == 0 || b.
cols () == 0)
5777 volatile int typ = mattype.
type ();
5783 cholmod_common Common;
5784 cholmod_common *cm = &Common;
5787 CHOLMOD_NAME(start) (cm);
5788 cm->prefer_zomplex =
false;
5794 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
5798 cm->print =
static_cast<int> (spu) + 2;
5799 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5803 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5804 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5806 cm->final_ll =
true;
5808 cholmod_sparse Astore;
5809 cholmod_sparse *
A = &Astore;
5820 #ifdef USE_64_BIT_IDX_T
5821 A->itype = CHOLMOD_LONG;
5823 A->itype = CHOLMOD_INT;
5825 A->dtype = CHOLMOD_DOUBLE;
5827 A->xtype = CHOLMOD_REAL;
5834 cholmod_dense Bstore;
5835 cholmod_dense *
B = &Bstore;
5836 B->nrow = b.
rows ();
5837 B->ncol = b.
cols ();
5839 B->nzmax = B->nrow * B->ncol;
5840 B->dtype = CHOLMOD_DOUBLE;
5841 B->xtype = CHOLMOD_REAL;
5842 if (nc < 1 || b.
cols () < 1)
5850 L = CHOLMOD_NAME(analyze) (
A, cm);
5853 rcond = CHOLMOD_NAME(rcond)(L, cm);
5867 volatile double rcond_plus_one = rcond + 1.0;
5869 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
5875 sing_handler (rcond);
5886 X = CHOLMOD_NAME(
solve) (CHOLMOD_A, L,
B, cm);
5894 retval.
xelem (i,j) =
static_cast<double *
>(X->x)[jr + i];
5898 CHOLMOD_NAME(free_dense) (&X, cm);
5899 CHOLMOD_NAME(free_factor) (&L, cm);
5900 CHOLMOD_NAME(finish) (cm);
5901 static char tmp[] =
" ";
5902 CHOLMOD_NAME(print_common) (tmp, cm);
5906 (*current_liboctave_warning_with_id_handler)
5907 (
"Octave:missing-dependency",
"CHOLMOD not installed");
5919 factorize (err, rcond, Control, Info, sing_handler, calc_cond);
5933 const double *Ax =
data ();
5938 Ai, Ax, &result[iidx],
5939 &Bx[iidx], Numeric, control,
5943 (*current_liboctave_error_handler)
5944 (
"SparseMatrix::solve solve failed");
5962 (*current_liboctave_error_handler) (
"UMFPACK not installed");
5966 (*current_liboctave_error_handler) (
"incorrect matrix type");
5975 solve_singularity_handler sing_handler,
5976 bool calc_cond)
const
5984 if (nr != nc || nr != b.
rows ())
5986 (
"matrix dimension mismatch solution of linear equations");
5987 else if (nr == 0 || b.
cols () == 0)
5992 volatile int typ = mattype.
type ();
5998 cholmod_common Common;
5999 cholmod_common *cm = &Common;
6002 CHOLMOD_NAME(start) (cm);
6003 cm->prefer_zomplex =
false;
6009 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
6013 cm->print =
static_cast<int> (spu) + 2;
6014 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6018 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6019 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6021 cm->final_ll =
true;
6023 cholmod_sparse Astore;
6024 cholmod_sparse *
A = &Astore;
6035 #ifdef USE_64_BIT_IDX_T
6036 A->itype = CHOLMOD_LONG;
6038 A->itype = CHOLMOD_INT;
6040 A->dtype = CHOLMOD_DOUBLE;
6042 A->xtype = CHOLMOD_REAL;
6049 cholmod_sparse Bstore;
6050 cholmod_sparse *
B = &Bstore;
6051 B->nrow = b.
rows ();
6052 B->ncol = b.
cols ();
6055 B->nzmax = b.
nnz ();
6059 #ifdef USE_64_BIT_IDX_T
6060 B->itype = CHOLMOD_LONG;
6062 B->itype = CHOLMOD_INT;
6064 B->dtype = CHOLMOD_DOUBLE;
6066 B->xtype = CHOLMOD_REAL;
6068 if (b.
rows () < 1 || b.
cols () < 1)
6075 L = CHOLMOD_NAME(analyze) (
A, cm);
6078 rcond = CHOLMOD_NAME(rcond)(L, cm);
6091 volatile double rcond_plus_one = rcond + 1.0;
6093 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
6099 sing_handler (rcond);
6110 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L,
B, cm);
6113 retval =
SparseMatrix (static_cast<octave_idx_type>(X->nrow),
6114 static_cast<octave_idx_type>(X->ncol),
6115 static_cast<octave_idx_type>(X->nzmax));
6117 j <= static_cast<octave_idx_type>(X->ncol); j++)
6120 j < static_cast<octave_idx_type>(X->nzmax); j++)
6123 retval.
xdata (j) =
static_cast<double *
>(X->x)[j];
6127 CHOLMOD_NAME(free_sparse) (&X, cm);
6128 CHOLMOD_NAME(free_factor) (&L, cm);
6129 CHOLMOD_NAME(finish) (cm);
6130 static char tmp[] =
" ";
6131 CHOLMOD_NAME(print_common) (tmp, cm);
6135 (*current_liboctave_warning_with_id_handler)
6136 (
"Octave:missing-dependency",
"CHOLMOD not installed");
6147 void *Numeric =
factorize (err, rcond, Control, Info,
6148 sing_handler, calc_cond);
6159 const double *Ax =
data ();
6170 retval.
xcidx (0) = 0;
6175 Bx[i] = b.
elem (i, j);
6178 Ai, Ax, Xx, Bx, Numeric,
6182 (*current_liboctave_error_handler)
6183 (
"SparseMatrix::solve solve failed");
6201 sz = (sz > 10 ? sz : 10) + x_nz;
6205 retval.
xdata (ii) = tmp;
6206 retval.
xridx (ii++) = i;
6209 retval.
xcidx (j+1) = ii;
6222 (*current_liboctave_error_handler) (
"UMFPACK not installed");
6226 (*current_liboctave_error_handler) (
"incorrect matrix type");
6235 solve_singularity_handler sing_handler,
6236 bool calc_cond)
const
6244 if (nr != nc || nr != b.
rows ())
6246 (
"matrix dimension mismatch solution of linear equations");
6247 else if (nr == 0 || b.
cols () == 0)
6252 volatile int typ = mattype.
type ();
6258 cholmod_common Common;
6259 cholmod_common *cm = &Common;
6262 CHOLMOD_NAME(start) (cm);
6263 cm->prefer_zomplex =
false;
6269 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
6273 cm->print =
static_cast<int> (spu) + 2;
6274 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6278 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6279 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6281 cm->final_ll =
true;
6283 cholmod_sparse Astore;
6284 cholmod_sparse *
A = &Astore;
6295 #ifdef USE_64_BIT_IDX_T
6296 A->itype = CHOLMOD_LONG;
6298 A->itype = CHOLMOD_INT;
6300 A->dtype = CHOLMOD_DOUBLE;
6302 A->xtype = CHOLMOD_REAL;
6309 cholmod_dense Bstore;
6310 cholmod_dense *
B = &Bstore;
6311 B->nrow = b.
rows ();
6312 B->ncol = b.
cols ();
6314 B->nzmax = B->nrow * B->ncol;
6315 B->dtype = CHOLMOD_DOUBLE;
6316 B->xtype = CHOLMOD_COMPLEX;
6317 if (nc < 1 || b.
cols () < 1)
6325 L = CHOLMOD_NAME(analyze) (
A, cm);
6328 rcond = CHOLMOD_NAME(rcond)(L, cm);
6341 volatile double rcond_plus_one = rcond + 1.0;
6343 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
6349 sing_handler (rcond);
6360 X = CHOLMOD_NAME(
solve) (CHOLMOD_A, L,
B, cm);
6368 retval.
xelem (i,j) =
static_cast<Complex *
>(X->x)[jr + i];
6372 CHOLMOD_NAME(free_dense) (&X, cm);
6373 CHOLMOD_NAME(free_factor) (&L, cm);
6374 CHOLMOD_NAME(finish) (cm);
6375 static char tmp[] =
" ";
6376 CHOLMOD_NAME(print_common) (tmp, cm);
6380 (*current_liboctave_warning_with_id_handler)
6381 (
"Octave:missing-dependency",
"CHOLMOD not installed");
6392 void *Numeric =
factorize (err, rcond, Control, Info,
6393 sing_handler, calc_cond);
6404 const double *Ax =
data ();
6409 retval.
resize (b_nr, b_nc);
6424 Ai, Ax, Xx, Bx, Numeric,
6428 Numeric, control, info);
6430 if (status < 0 || status2 < 0)
6432 (*current_liboctave_error_handler)
6433 (
"SparseMatrix::solve solve failed");
6443 retval(i, j) =
Complex (Xx[i], Xz[i]);
6454 (*current_liboctave_error_handler) (
"UMFPACK not installed");
6458 (*current_liboctave_error_handler) (
"incorrect matrix type");
6467 solve_singularity_handler sing_handler,
6468 bool calc_cond)
const
6476 if (nr != nc || nr != b.
rows ())
6478 (
"matrix dimension mismatch solution of linear equations");
6479 else if (nr == 0 || b.
cols () == 0)
6484 volatile int typ = mattype.
type ();
6490 cholmod_common Common;
6491 cholmod_common *cm = &Common;
6494 CHOLMOD_NAME(start) (cm);
6495 cm->prefer_zomplex =
false;
6501 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
6505 cm->print =
static_cast<int> (spu) + 2;
6506 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6510 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6511 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6513 cm->final_ll =
true;
6515 cholmod_sparse Astore;
6516 cholmod_sparse *
A = &Astore;
6527 #ifdef USE_64_BIT_IDX_T
6528 A->itype = CHOLMOD_LONG;
6530 A->itype = CHOLMOD_INT;
6532 A->dtype = CHOLMOD_DOUBLE;
6534 A->xtype = CHOLMOD_REAL;
6541 cholmod_sparse Bstore;
6542 cholmod_sparse *
B = &Bstore;
6543 B->nrow = b.
rows ();
6544 B->ncol = b.
cols ();
6547 B->nzmax = b.
nnz ();
6551 #ifdef USE_64_BIT_IDX_T
6552 B->itype = CHOLMOD_LONG;
6554 B->itype = CHOLMOD_INT;
6556 B->dtype = CHOLMOD_DOUBLE;
6558 B->xtype = CHOLMOD_COMPLEX;
6560 if (b.
rows () < 1 || b.
cols () < 1)
6567 L = CHOLMOD_NAME(analyze) (
A, cm);
6570 rcond = CHOLMOD_NAME(rcond)(L, cm);
6583 volatile double rcond_plus_one = rcond + 1.0;
6585 if (rcond_plus_one == 1.0 ||
xisnan (rcond))
6591 sing_handler (rcond);
6602 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L,
B, cm);
6606 (static_cast<octave_idx_type>(X->nrow),
6607 static_cast<octave_idx_type>(X->ncol),
6608 static_cast<octave_idx_type>(X->nzmax));
6610 j <= static_cast<octave_idx_type>(X->ncol); j++)
6613 j < static_cast<octave_idx_type>(X->nzmax); j++)
6620 CHOLMOD_NAME(free_sparse) (&X, cm);
6621 CHOLMOD_NAME(free_factor) (&L, cm);
6622 CHOLMOD_NAME(finish) (cm);
6623 static char tmp[] =
" ";
6624 CHOLMOD_NAME(print_common) (tmp, cm);
6628 (*current_liboctave_warning_with_id_handler)
6629 (
"Octave:missing-dependency",
"CHOLMOD not installed");
6640 void *Numeric =
factorize (err, rcond, Control, Info,
6641 sing_handler, calc_cond);
6652 const double *Ax =
data ();
6666 retval.
xcidx (0) = 0;
6677 Ai, Ax, Xx, Bx, Numeric,
6681 Numeric, control, info);
6683 if (status < 0 || status2 < 0)
6685 (*current_liboctave_error_handler)
6686 (
"SparseMatrix::solve solve failed");
6704 sz = (sz > 10 ? sz : 10) + x_nz;
6708 retval.
xdata (ii) = tmp;
6709 retval.
xridx (ii++) = i;
6712 retval.
xcidx (j+1) = ii;
6724 (*current_liboctave_error_handler) (
"UMFPACK not installed");
6728 (*current_liboctave_error_handler) (
"incorrect matrix type");
6739 return solve (mattype, b, info, rcond, 0);
6747 return solve (mattype, b, info, rcond, 0);
6754 return solve (mattype, b, info, rcond, 0);
6759 double& rcond, solve_singularity_handler sing_handler,
6760 bool singular_fallback)
const
6763 int typ = mattype.
type (
false);
6766 typ = mattype.
type (*
this);
6770 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6772 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6774 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6776 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6779 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6781 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6784 (*current_liboctave_error_handler) (
"unknown matrix type");
6793 retval =
qrsolve (*
this, b, err);
6795 retval = dmsolve<Matrix, SparseMatrix, Matrix> (*
this, b, err);
6807 return solve (mattype, b, info, rcond, 0);
6815 return solve (mattype, b, info, rcond, 0);
6822 return solve (mattype, b, info, rcond, 0);
6828 solve_singularity_handler sing_handler,
6829 bool singular_fallback)
const
6832 int typ = mattype.
type (
false);
6835 typ = mattype.
type (*
this);
6838 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6840 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6842 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6844 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6847 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6849 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6852 (*current_liboctave_error_handler) (
"unknown matrix type");
6860 retval =
qrsolve (*
this, b, err);
6862 retval = dmsolve<SparseMatrix, SparseMatrix, SparseMatrix>
6875 return solve (mattype, b, info, rcond, 0);
6883 return solve (mattype, b, info, rcond, 0);
6890 return solve (mattype, b, info, rcond, 0);
6896 solve_singularity_handler sing_handler,
6897 bool singular_fallback)
const
6900 int typ = mattype.
type (
false);
6903 typ = mattype.
type (*
this);
6906 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6908 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6910 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6912 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6915 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6917 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6920 (*current_liboctave_error_handler) (
"unknown matrix type");
6928 retval =
qrsolve (*
this, b, err);
6930 retval = dmsolve<ComplexMatrix, SparseMatrix, ComplexMatrix>
6943 return solve (mattype, b, info, rcond, 0);
6951 return solve (mattype, b, info, rcond, 0);
6958 return solve (mattype, b, info, rcond, 0);
6964 solve_singularity_handler sing_handler,
6965 bool singular_fallback)
const
6968 int typ = mattype.
type (
false);
6971 typ = mattype.
type (*
this);
6974 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6976 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6978 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6980 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6983 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6985 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6988 (*current_liboctave_error_handler) (
"unknown matrix type");
6996 retval =
qrsolve (*
this, b, err);
6998 retval = dmsolve<SparseComplexMatrix, SparseMatrix, SparseComplexMatrix>
7010 return solve (mattype, b, info, rcond);
7018 return solve (mattype, b, info, rcond);
7025 return solve (mattype, b, info, rcond, 0);
7031 solve_singularity_handler sing_handler)
const
7034 return solve (mattype, tmp, info, rcond,
7035 sing_handler).
column (static_cast<octave_idx_type> (0));
7043 return solve (mattype, b, info, rcond, 0);
7051 return solve (mattype, b, info, rcond, 0);
7057 double& rcond)
const
7059 return solve (mattype, b, info, rcond, 0);
7065 solve_singularity_handler sing_handler)
const
7068 return solve (mattype, tmp, info, rcond,
7069 sing_handler).
column (static_cast<octave_idx_type> (0));
7077 return solve (b, info, rcond, 0);
7084 return solve (b, info, rcond, 0);
7089 double& rcond)
const
7091 return solve (b, info, rcond, 0);
7096 solve_singularity_handler sing_handler)
const
7099 return solve (mattype, b, err, rcond, sing_handler);
7107 return solve (b, info, rcond, 0);
7115 return solve (b, info, rcond, 0);
7122 return solve (b, info, rcond, 0);
7127 solve_singularity_handler sing_handler)
const
7130 return solve (mattype, b, err, rcond, sing_handler);
7137 return solve (b, info, rcond, 0);
7142 double& rcond)
const
7144 return solve (b, info, rcond, 0);
7150 solve_singularity_handler sing_handler)
const
7153 return solve (mattype, b, err, rcond, sing_handler);
7161 return solve (b, info, rcond, 0);
7168 return solve (b, info, rcond, 0);
7173 double& rcond)
const
7175 return solve (b, info, rcond, 0);
7181 solve_singularity_handler sing_handler)
const
7184 return solve (mattype, b, err, rcond, sing_handler);
7191 return solve (b, info, rcond);
7198 return solve (b, info, rcond);
7203 double& rcond)
const
7205 return solve (b, info, rcond, 0);
7211 solve_singularity_handler sing_handler)
const
7214 return solve (tmp, info, rcond,
7215 sing_handler).
column (static_cast<octave_idx_type> (0));
7223 return solve (b, info, rcond, 0);
7230 return solve (b, info, rcond, 0);
7235 double& rcond)
const
7237 return solve (b, info, rcond, 0);
7243 solve_singularity_handler sing_handler)
const
7246 return solve (tmp, info, rcond,
7247 sing_handler).
column (static_cast<octave_idx_type> (0));
7280 double val =
data (i);
7295 double val =
data (i);
7310 double val =
data (i);
7311 if (val != 0.0 && val != 1.0)
7337 double val =
data (i);
7363 double val =
data (i);
7404 if (jj <
cidx (i+1) &&
ridx (jj) == j)
7447 if ((
rows () == 1 && dim == -1) || dim == 1)
7452 (
cidx (j+1) -
cidx (j) < nr ? 0.0 : 1.0), 1.0);
7466 double d = data (i); \
7467 tmp[ridx (i)] += d * d
7470 double d = data (i); \
7488 retval.
data (i) = fabs (retval.
data (i));
7517 os << a.
ridx (i) + 1 <<
" " << j + 1 <<
" ";
7531 return read_sparse_matrix<elt_type> (is, a, octave_read_value<double>);
7595 return do_mul_dm_sm<SparseMatrix> (
d, a);
7601 return do_mul_sm_dm<SparseMatrix> (a,
d);
7607 return do_add_dm_sm<SparseMatrix> (
d, a);
7613 return do_sub_dm_sm<SparseMatrix> (
d, a);
7619 return do_add_sm_dm<SparseMatrix> (a,
d);
7625 return do_sub_sm_dm<SparseMatrix> (a,
d);
7644 #define EMPTY_RETURN_CHECK(T) \
7645 if (nr == 0 || nc == 0) \
7669 result.
xdata (idx) = tmp;
7685 result.
xcidx (0) = 0;
7694 result.
xdata (ii) = tmp;
7698 result.
xcidx (j+1) = ii;
7721 if (a_nr == b_nr && a_nc == b_nc)
7731 bool ja_lt_max= ja < ja_max;
7735 bool jb_lt_max = jb < jb_max;
7737 while (ja_lt_max || jb_lt_max)
7740 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
7742 double tmp =
xmin (a.
data (ja), 0.);
7750 ja_lt_max= ja < ja_max;
7752 else if ((! ja_lt_max)
7753 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
7755 double tmp =
xmin (0., b.
data (jb));
7763 jb_lt_max= jb < jb_max;
7775 ja_lt_max= ja < ja_max;
7777 jb_lt_max= jb < jb_max;
7787 if (a_nr == 0 || a_nc == 0)
7789 else if (b_nr == 0 || b_nc == 0)
7820 result.
xdata (idx) = tmp;
7836 result.
xcidx (0) = 0;
7844 result.
xdata (ii) = tmp;
7848 result.
xcidx (j+1) = ii;
7871 if (a_nr == b_nr && a_nc == b_nc)
7881 bool ja_lt_max= ja < ja_max;
7885 bool jb_lt_max = jb < jb_max;
7887 while (ja_lt_max || jb_lt_max)
7890 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
7892 double tmp =
xmax (a.
data (ja), 0.);
7900 ja_lt_max= ja < ja_max;
7902 else if ((! ja_lt_max)
7903 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
7905 double tmp =
xmax (0., b.
data (jb));
7913 jb_lt_max= jb < jb_max;
7925 ja_lt_max= ja < ja_max;
7927 jb_lt_max= jb < jb_max;
7937 if (a_nr == 0 || a_nc == 0)
7939 else if (b_nr == 0 || b_nc == 0)
void octave_write_double(std::ostream &os, double d)
octave_idx_type * xridx(void)
SparseMatrix min(int dim=-1) const
octave_idx_type cols(void) const
#define F77_CHAR_ARG_LEN(l)
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
SparseMatrix reshape(const dim_vector &new_dims) const
bool operator==(const SparseMatrix &a) const
SparseMatrix concat(const SparseMatrix &rb, const Array< octave_idx_type > &ra_idx)
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type octave_idx_type octave_idx_type &F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type const octave_idx_type const double const octave_idx_type const octave_idx_type double const octave_idx_type octave_idx_type & F77_CHAR_ARG_LEN_DECL
octave_idx_type rows(void) const
#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)
double & elem(octave_idx_type n)
SparseBoolMatrix all(int dim=-1) const
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
SparseMatrix transpose(void) const
dim_vector dims(void) const
#define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, INIT_VAL, MT_RESULT)
ColumnVector column(octave_idx_type i) const
bool all_elements_are_zero(void) const
SparseMatrix sum(int dim=-1) const
Matrix dsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
Matrix bsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
void SparseCholError(int status, char *file, int line, char *message)
SparseMatrix inverse(void) const
F77_RET_T const octave_idx_type Complex * A
SM octinternal_do_mul_pm_sm(const PermMatrix &p, const SM &a)
Complex xmax(const Complex &x, const Complex &y)
Matrix mul_trans(const Matrix &m, const SparseMatrix &a)
#define SPARSE_SMSM_BOOL_OPS(M1, M2, ZERO)
SparseMatrix ipermute(const Array< octave_idx_type > &vec) const
octave_idx_type * xcidx(void)
F77_RET_T F77_FUNC(dgbtrf, DGBTRF)(const octave_idx_type &
SparseBoolMatrix operator!(void) const
SparseMatrix & insert(const SparseMatrix &a, octave_idx_type r, octave_idx_type c)
RowVector row(octave_idx_type i) 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)
Complex xmin(const Complex &x, const Complex &y)
octave_idx_type * cidx(void)
T & elem(octave_idx_type n)
bool is_hermitian(void) const
#define SPARSE_REDUCTION_OP(RET_TYPE, EL_TYPE, OP, INIT_VAL, MT_RESULT)
octave_idx_type columns(void) const
SparseMatrix max(int dim=-1) const
std::istream & operator>>(std::istream &is, SparseMatrix &a)
bool is_dense(void) const
SparseMatrix abs(void) const
#define FULL_SPARSE_MUL(RET_TYPE, EL_TYPE, ZERO)
bool operator!=(const SparseMatrix &a) const
#define F77_XFCN(f, F, args)
RowVector row(octave_idx_type i) const
static double get_key(const std::string &key)
bool is_symmetric(void) const
Matrix solve(MatrixType &typ, const Matrix &b) const
SparseMatrix cumprod(int dim=-1) const
int first_non_singleton(int def=0) const
octave_idx_type rows(void) const
#define SPARSE_SMS_CMP_OPS(M, MZ, CM, S, SZ, CS)
F77_RET_T const double const double double * d
#define SPARSE_FULL_MUL(RET_TYPE, EL_TYPE, ZERO)
#define F77_CONST_CHAR_ARG2(x, l)
octave_idx_type * triangular_perm(void) const
octave_idx_type nnz(void) const
Matrix ltsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
#define SPARSE_ALL_OP(DIM)
MSparse< T > squeeze(void) const
SparseMatrix Q(void) const
SparseMatrix imag(const SparseComplexMatrix &a)
#define SPARSE_ANY_OP(DIM)
liboctave_error_handler current_liboctave_error_handler
int type(bool quiet=true)
Matrix trans_mul(const SparseMatrix &m, const Matrix &a)
#define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN)
DET determinant(void) const
bool all_integers(double &max_val, double &min_val) const
#define FULL_SPARSE_MUL_TRANS(RET_TYPE, EL_TYPE, ZERO, CONJ_OP)
Matrix matrix_value(void) const
void change_capacity(octave_idx_type nz)
SparseMatrix diag(octave_idx_type k=0) const
MatrixType transpose(void) const
#define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE)
Sparse< T > maybe_compress(bool remove_zeros=false)
void gripe_nan_to_logical_conversion(void)
void resize(octave_idx_type r, octave_idx_type c)
bool any_element_is_negative(bool=false) const
MSparse< T > & insert(const Sparse< T > &a, octave_idx_type r, octave_idx_type c)
SparseMatrix sumsq(int dim=-1) const
void resize(const dim_vector &dv, const T &rfv)
void mark_as_unsymmetric(void)
bool too_large_for_float(void) const
MSparse< T > permute(const Array< octave_idx_type > &vec, bool inv=false) const
SparseBoolMatrix any(int dim=-1) const
SparseMatrix real(const SparseComplexMatrix &a)
bool any_element_not_one_or_zero(void) const
void gripe_singular_matrix(double rcond)
SparseMatrix L(void) const
SparseMatrix operator-(const DiagMatrix &d, const SparseMatrix &a)
SparseMatrix cumsum(int dim=-1) const
MSparse< T > ipermute(const Array< octave_idx_type > &vec) const
SparseMatrix min(double d, const SparseMatrix &m)
SparseComplexMatrix & insert(const SparseComplexMatrix &a, octave_idx_type r, octave_idx_type c)
OCTAVE_API double D_NINT(double x)
SparseMatrix operator+(const DiagMatrix &d, const SparseMatrix &a)
void mark_as_rectangular(void)
#define EMPTY_RETURN_CHECK(T)
octave_idx_type nnz(void) const
Count nonzero elements.
#define SPARSE_SMSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2)
T & xelem(octave_idx_type n)
SparseMatrix squeeze(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
octave_idx_type cols(void) const
Matrix utsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
int SparseCholPrint(const char *fmt,...)
SparseMatrix atan2(const double &x, const SparseMatrix &y)
octave_idx_type * ridx(void)
F77_RET_T const octave_idx_type Complex const octave_idx_type Complex * B
#define END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
ComplexMatrix qrsolve(const SparseComplexMatrix &a, const Matrix &b, octave_idx_type &info)
#define F77_CONST_CHAR_ARG_DECL
#define SPARSE_SSM_CMP_OPS(S, SZ, SC, M, MZ, MC)
bool test_any(F fcn) const
bool any_element_is_inf_or_nan(void) const
bool all_elements_are_int_or_inf_or_nan(void) const
void * factorize(octave_idx_type &err, double &rcond, Matrix &Control, Matrix &Info, solve_singularity_handler sing_handler, bool calc_cond=false) const
ColumnVector column(octave_idx_type i) const
SparseMatrix permute(const Array< octave_idx_type > &vec, bool inv=false) const
SparseMatrix prod(int dim=-1) const
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
SparseMatrix operator*(const SparseMatrix &m, const SparseMatrix &a)
bool any_element_is_nan(void) const
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double * Q
bool xtoo_large_for_float(double x)
#define SPARSE_SMS_BOOL_OPS(M, S, ZERO)
#define BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
Matrix sum(int dim=-1) const
std::complex< double > Complex
const T * fortran_vec(void) const
MSparse< T > diag(octave_idx_type k=0) const
#define SPARSE_SSM_BOOL_OPS(S, M, ZERO)
octave_idx_type cols(void) const
SparseMatrix dinverse(MatrixType &mattyp, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
Matrix fsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
SparseMatrix tinverse(MatrixType &mattyp, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
octave_idx_type length(void) const
#define UMFPACK_DNAME(name)
F77_RET_T const double * x
MSparse< T > reshape(const dim_vector &new_dims) const
Array< T > array_value(void) const
std::ostream & operator<<(std::ostream &os, const SparseMatrix &a)