26 #if defined (HAVE_CONFIG_H)
61 #if ! defined (USE_QRSOLVE)
66 :
MSparse<double> (a.rows (), a.cols (), a.nnz ())
82 :
MSparse<double> (a.rows (), a.cols (), a.length ())
110 if (nr != nr_a || nc != nc_a || nz != nz_a)
127 return !(*
this == a);
136 if (nr == nc && nr > 0)
189 return max (dummy_idx, dim);
200 if (dim >= dv.
ndims ())
213 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
223 if (
ridx (i) != idx_j)
234 double tmp =
data (i);
254 result.
xcidx (0) = 0;
257 double tmp =
elem (idx_arg(j), j);
260 result.
xdata (ii) = tmp;
261 result.
xridx (ii++) = 0;
263 result.
xcidx (j+1) = ii;
271 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
281 if (found[
ridx (i)] == -j)
282 found[
ridx (i)] = -j - 1;
285 if (found[i] > -nc && found[i] < 0)
286 idx_arg.
elem (i) = -found[i];
294 double tmp =
data (i);
298 else if (ix == -1 || tmp >
elem (ir, ix))
299 idx_arg.
elem (ir) = j;
305 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
311 result.
xcidx (0) = 0;
312 result.
xcidx (1) = nel;
315 if (idx_arg(j) == -1)
319 result.
xridx (ii++) = j;
323 double tmp =
elem (j, idx_arg(j));
326 result.
xdata (ii) = tmp;
327 result.
xridx (ii++) = j;
340 return min (dummy_idx, dim);
351 if (dim >= dv.
ndims ())
364 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
374 if (
ridx (i) != idx_j)
385 double tmp =
data (i);
405 result.
xcidx (0) = 0;
408 double tmp =
elem (idx_arg(j), j);
411 result.
xdata (ii) = tmp;
412 result.
xridx (ii++) = 0;
414 result.
xcidx (j+1) = ii;
422 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
432 if (found[
ridx (i)] == -j)
433 found[
ridx (i)] = -j - 1;
436 if (found[i] > -nc && found[i] < 0)
437 idx_arg.
elem (i) = -found[i];
445 double tmp =
data (i);
449 else if (ix == -1 || tmp <
elem (ir, ix))
450 idx_arg.
elem (ir) = j;
456 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
462 result.
xcidx (0) = 0;
463 result.
xcidx (1) = nel;
466 if (idx_arg(j) == -1)
470 result.
xridx (ii++) = j;
474 double tmp =
elem (j, idx_arg(j));
477 result.
xdata (ii) = tmp;
478 result.
xridx (ii++) = j;
513 retval(j) =
data (k);
538 if (rb.
rows () > 0 && rb.
cols () > 0)
548 if (rb.
rows () > 0 && rb.
cols () > 0)
562 r.cidx (i) = a.
cidx (i);
567 r.ridx (i) = a.
ridx (i);
570 r.maybe_compress (
true);
583 r.cidx (i) = a.
cidx (i);
588 r.ridx (i) = a.
ridx (i);
591 r.maybe_compress (
true);
608 return inverse (mattype, info, rcond, 0, 0);
616 return inverse (mattype, info, rcond, 0, 0);
623 return inverse (mattype, info, rcond, 0, 0);
628 double& rcond,
const bool,
629 const bool calccond)
const
637 if (nr == 0 || nc == 0 || nr != nc)
638 (*current_liboctave_error_handler) (
"inverse requires square matrix");
641 int typ = mattype.
type ();
645 (*current_liboctave_error_handler) (
"incorrect matrix type");
653 double *v = retval.
data ();
661 double tmp = fabs (v[i]);
678 double& rcond,
const bool,
679 const bool calccond)
const
687 if (nr == 0 || nc == 0 || nr != nc)
688 (*current_liboctave_error_handler) (
"inverse requires square matrix");
691 int typ = mattype.
type ();
696 (*current_liboctave_error_handler) (
"incorrect matrix type");
699 double ainvnorm = 0.;
708 atmp += fabs (
data (i));
733 retval.
xcidx (i) = cx;
734 retval.
xridx (cx) = i;
735 retval.
xdata (cx) = 1.0;
748 (*current_liboctave_error_handler) (
"division by zero");
753 rpX = retval.
xridx (colXp);
762 v -= retval.
xdata (colXp) *
data (colUp);
767 while (rpX < j && rpU < j && colXp < cx && colUp < nz);
771 colUp =
cidx (j+1) - 1;
774 double pivot =
data (colUp);
775 if (pivot == 0. ||
ridx (colUp) != j)
776 (*current_liboctave_error_handler) (
"division by zero");
786 retval.
xridx (cx) = j;
787 retval.
xdata (cx) = v / pivot;
795 colUp =
cidx (i+1) - 1;
798 double pivot =
data (colUp);
799 if (pivot == 0. ||
ridx (colUp) != i)
800 (*current_liboctave_error_handler) (
"division by zero");
804 retval.
xdata (j) /= pivot;
806 retval.
xcidx (nr) = cx;
851 k <
cidx (jidx+1); k++)
864 (*current_liboctave_error_handler) (
"division by zero");
872 colUp =
cidx (perm[iidx]+1) - 1;
874 colUp =
cidx (perm[iidx]);
876 double pivot =
data (colUp);
878 (*current_liboctave_error_handler) (
"division by zero");
891 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
895 retval.
xcidx (i) = cx;
899 retval.
xridx (cx) = j;
900 retval.
xdata (cx++) = work[j];
904 retval.
xcidx (nr) = cx;
915 i < retval.
cidx (j+1); i++)
916 atmp += fabs (retval.
data (i));
921 rcond = 1. / ainvnorm / anorm;
929 double& rcond,
bool,
bool calc_cond)
const
933 (*current_liboctave_error_handler)
934 (
"inverse of the null matrix not defined");
937 int typ = mattype.
type (
false);
941 typ = mattype.
type (*
this);
944 ret =
dinverse (mattype, info, rcond,
true, calc_cond);
958 rcond = fact.rcond ();
964 info, rcond2,
true,
false);
965 ret =
Q * InvL.
transpose () * InvL *
Q.transpose ();
985 rcond = fact.rcond ();
1001 info, rcond2,
true,
false);
1004 ret = fact.Pc ().
transpose () * InvU * InvL * fact.Pr ();
1031 #if defined (HAVE_UMFPACK)
1036 if (nr == 0 || nc == 0 || nr != nc)
1045 Matrix Control (UMFPACK_CONTROL, 1);
1049 double tmp = octave::sparse_params::get_key (
"spumoni");
1051 Control (UMFPACK_PRL) = tmp;
1053 tmp = octave::sparse_params::get_key (
"piv_tol");
1056 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
1057 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
1061 tmp = octave::sparse_params::get_key (
"autoamd");
1063 Control (UMFPACK_FIXQ) = tmp;
1066 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1072 const double *Ax =
data ();
1080 Matrix Info (1, UMFPACK_INFO);
1085 Ax,
nullptr, &Symbolic, control, info);
1094 (*current_liboctave_error_handler)
1095 (
"SparseMatrix::determinant symbolic factorization failed");
1105 &Numeric, control, info);
1108 rcond = Info (UMFPACK_RCOND);
1116 (*current_liboctave_error_handler)
1117 (
"SparseMatrix::determinant numeric factorization failed");
1125 status =
UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric,
1133 (*current_liboctave_error_handler)
1134 (
"SparseMatrix::determinant error calculating determinant");
1137 retval =
DET (c10, e10, 10);
1146 octave_unused_parameter (err);
1147 octave_unused_parameter (rcond);
1149 (*current_liboctave_error_handler)
1150 (
"support for UMFPACK was unavailable or disabled "
1151 "when liboctave was built");
1161 double& rcond, solve_singularity_handler,
1162 bool calc_cond)
const
1171 if (nr != b.
rows ())
1173 (
"matrix dimension mismatch solution of linear equations");
1175 if (nr == 0 || nc == 0 || b.
cols () == 0)
1180 int typ = mattype.
type ();
1184 (*current_liboctave_error_handler) (
"incorrect matrix type");
1190 retval(i, j) = b(i, j) /
data (i);
1195 retval(k, j) = b(
ridx (i), j) /
data (i);
1203 double tmp = fabs (
data (i));
1209 rcond = dmin / dmax;
1221 solve_singularity_handler,
bool calc_cond)
const
1230 if (nr != b.
rows ())
1232 (
"matrix dimension mismatch solution of linear equations");
1234 if (nr == 0 || nc == 0 || b.
cols () == 0)
1239 int typ = mattype.
type ();
1243 (*current_liboctave_error_handler) (
"incorrect matrix type");
1249 retval.
xcidx (0) = 0;
1256 if (b.
ridx (i) >= nm)
1261 retval.
xcidx (j+1) = ii;
1271 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1279 retval.
xridx (ii) = l;
1283 retval.
xcidx (j+1) = ii;
1292 double tmp = fabs (
data (i));
1298 rcond = dmin / dmax;
1310 solve_singularity_handler,
bool calc_cond)
const
1319 if (nr != b.
rows ())
1321 (
"matrix dimension mismatch solution of linear equations");
1323 if (nr == 0 || nc == 0 || b.
cols () == 0)
1328 int typ = mattype.
type ();
1332 (*current_liboctave_error_handler) (
"incorrect matrix type");
1338 retval(i, j) = b(i, j) /
data (i);
1343 retval(k, j) = b(
ridx (i), j) /
data (i);
1351 double tmp = fabs (
data (i));
1357 rcond = dmin / dmax;
1369 solve_singularity_handler,
bool calc_cond)
const
1378 if (nr != b.
rows ())
1380 (
"matrix dimension mismatch solution of linear equations");
1382 if (nr == 0 || nc == 0 || b.
cols () == 0)
1387 int typ = mattype.
type ();
1391 (*current_liboctave_error_handler) (
"incorrect matrix type");
1397 retval.
xcidx (0) = 0;
1404 if (b.
ridx (i) >= nm)
1409 retval.
xcidx (j+1) = ii;
1419 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1427 retval.
xridx (ii) = l;
1431 retval.
xcidx (j+1) = ii;
1440 double tmp = fabs (
data (i));
1446 rcond = dmin / dmax;
1458 solve_singularity_handler sing_handler,
1459 bool calc_cond)
const
1468 if (nr != b.
rows ())
1470 (
"matrix dimension mismatch solution of linear equations");
1472 if (nr == 0 || nc == 0 || b.
cols () == 0)
1477 int typ = mattype.
type ();
1481 (*current_liboctave_error_handler) (
"incorrect matrix type");
1484 double ainvnorm = 0.;
1495 atmp += fabs (
data (i));
1503 retval.
resize (nc, b_nc);
1524 goto triangular_error;
1527 double tmp = work[k] /
data (
cidx (kidx+1)-1);
1530 i <
cidx (kidx+1)-1; i++)
1533 work[iidx] = work[iidx] - tmp *
data (i);
1539 retval.
xelem (perm[i], j) = work[i];
1558 double tmp = work[k] /
data (
cidx (iidx+1)-1);
1561 i <
cidx (iidx+1)-1; i++)
1564 work[idx2] = work[idx2] - tmp *
data (i);
1571 atmp += fabs (work[i]);
1574 if (atmp > ainvnorm)
1577 rcond = 1. / ainvnorm / anorm;
1583 retval.
resize (nc, b_nc);
1600 goto triangular_error;
1603 double tmp = work[k] /
data (
cidx (k+1)-1);
1608 work[iidx] = work[iidx] - tmp *
data (i);
1614 retval.
xelem (i, j) = work[i];
1631 double tmp = work[k] /
data (
cidx (k+1)-1);
1636 work[iidx] = work[iidx] - tmp *
data (i);
1643 atmp += fabs (work[i]);
1646 if (atmp > ainvnorm)
1649 rcond = 1. / ainvnorm / anorm;
1658 sing_handler (rcond);
1665 volatile double rcond_plus_one = rcond + 1.0;
1673 sing_handler (rcond);
1687 solve_singularity_handler sing_handler,
1688 bool calc_cond)
const
1697 if (nr != b.
rows ())
1699 (
"matrix dimension mismatch solution of linear equations");
1701 if (nr == 0 || nc == 0 || b.
cols () == 0)
1706 int typ = mattype.
type ();
1710 (*current_liboctave_error_handler) (
"incorrect matrix type");
1713 double ainvnorm = 0.;
1723 atmp += fabs (
data (i));
1732 retval.
xcidx (0) = 0;
1762 goto triangular_error;
1765 double tmp = work[k] /
data (
cidx (kidx+1)-1);
1768 i <
cidx (kidx+1)-1; i++)
1771 work[iidx] = work[iidx] - tmp *
data (i);
1783 if (ii + new_nnz > x_nz)
1792 if (work[rperm[i]] != 0.)
1794 retval.
xridx (ii) = i;
1795 retval.
xdata (ii++) = work[rperm[i]];
1797 retval.
xcidx (j+1) = ii;
1818 double tmp = work[k] /
data (
cidx (iidx+1)-1);
1821 i <
cidx (iidx+1)-1; i++)
1824 work[idx2] = work[idx2] - tmp *
data (i);
1831 atmp += fabs (work[i]);
1834 if (atmp > ainvnorm)
1837 rcond = 1. / ainvnorm / anorm;
1859 goto triangular_error;
1862 double tmp = work[k] /
data (
cidx (k+1)-1);
1867 work[iidx] = work[iidx] - tmp *
data (i);
1879 if (ii + new_nnz > x_nz)
1890 retval.
xridx (ii) = i;
1891 retval.
xdata (ii++) = work[i];
1893 retval.
xcidx (j+1) = ii;
1912 double tmp = work[k] /
data (
cidx (k+1)-1);
1915 i <
cidx (k+1)-1; i++)
1918 work[iidx] = work[iidx] - tmp *
data (i);
1925 atmp += fabs (work[i]);
1928 if (atmp > ainvnorm)
1931 rcond = 1. / ainvnorm / anorm;
1940 sing_handler (rcond);
1947 volatile double rcond_plus_one = rcond + 1.0;
1955 sing_handler (rcond);
1968 solve_singularity_handler sing_handler,
1969 bool calc_cond)
const
1978 if (nr != b.
rows ())
1980 (
"matrix dimension mismatch solution of linear equations");
1982 if (nr == 0 || nc == 0 || b.
cols () == 0)
1987 int typ = mattype.
type ();
1991 (*current_liboctave_error_handler) (
"incorrect matrix type");
1994 double ainvnorm = 0.;
2005 atmp += fabs (
data (i));
2013 retval.
resize (nc, b_nc);
2034 goto triangular_error;
2040 i <
cidx (kidx+1)-1; i++)
2043 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2049 retval.
xelem (perm[i], j) = cwork[i];
2069 double tmp = work[k] /
data (
cidx (iidx+1)-1);
2072 i <
cidx (iidx+1)-1; i++)
2075 work[idx2] = work[idx2] - tmp *
data (i);
2082 atmp += fabs (work[i]);
2085 if (atmp > ainvnorm)
2088 rcond = 1. / ainvnorm / anorm;
2094 retval.
resize (nc, b_nc);
2111 goto triangular_error;
2119 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2125 retval.
xelem (i, j) = cwork[i];
2143 double tmp = work[k] /
data (
cidx (k+1)-1);
2146 i <
cidx (k+1)-1; i++)
2149 work[iidx] = work[iidx] - tmp *
data (i);
2156 atmp += fabs (work[i]);
2159 if (atmp > ainvnorm)
2162 rcond = 1. / ainvnorm / anorm;
2171 sing_handler (rcond);
2178 volatile double rcond_plus_one = rcond + 1.0;
2186 sing_handler (rcond);
2200 solve_singularity_handler sing_handler,
2201 bool calc_cond)
const
2210 if (nr != b.
rows ())
2212 (
"matrix dimension mismatch solution of linear equations");
2214 if (nr == 0 || nc == 0 || b.
cols () == 0)
2219 int typ = mattype.
type ();
2223 (*current_liboctave_error_handler) (
"incorrect matrix type");
2226 double ainvnorm = 0.;
2236 atmp += fabs (
data (i));
2245 retval.
xcidx (0) = 0;
2275 goto triangular_error;
2281 i <
cidx (kidx+1)-1; i++)
2284 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2296 if (ii + new_nnz > x_nz)
2305 if (cwork[rperm[i]] != 0.)
2307 retval.
xridx (ii) = i;
2308 retval.
xdata (ii++) = cwork[rperm[i]];
2310 retval.
xcidx (j+1) = ii;
2332 double tmp = work[k] /
data (
cidx (iidx+1)-1);
2335 i <
cidx (iidx+1)-1; i++)
2338 work[idx2] = work[idx2] - tmp *
data (i);
2345 atmp += fabs (work[i]);
2348 if (atmp > ainvnorm)
2351 rcond = 1. / ainvnorm / anorm;
2373 goto triangular_error;
2381 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2393 if (ii + new_nnz > x_nz)
2404 retval.
xridx (ii) = i;
2405 retval.
xdata (ii++) = cwork[i];
2407 retval.
xcidx (j+1) = ii;
2427 double tmp = work[k] /
data (
cidx (k+1)-1);
2430 i <
cidx (k+1)-1; i++)
2433 work[iidx] = work[iidx] - tmp *
data (i);
2440 atmp += fabs (work[i]);
2443 if (atmp > ainvnorm)
2446 rcond = 1. / ainvnorm / anorm;
2455 sing_handler (rcond);
2462 volatile double rcond_plus_one = rcond + 1.0;
2470 sing_handler (rcond);
2484 solve_singularity_handler sing_handler,
2485 bool calc_cond)
const
2494 if (nr != b.
rows ())
2496 (
"matrix dimension mismatch solution of linear equations");
2498 if (nr == 0 || nc == 0 || b.
cols () == 0)
2503 int typ = mattype.
type ();
2507 (*current_liboctave_error_handler) (
"incorrect matrix type");
2510 double ainvnorm = 0.;
2521 atmp += fabs (
data (i));
2529 retval.
resize (nc, b_nc);
2539 work[perm[i]] = b(i, j);
2549 if (perm[
ridx (i)] < minr)
2551 minr = perm[
ridx (i)];
2555 if (minr != k ||
data (mini) == 0)
2558 goto triangular_error;
2561 double tmp = work[k] /
data (mini);
2569 work[iidx] = work[iidx] - tmp *
data (i);
2575 retval(i, j) = work[i];
2596 i <
cidx (k+1); i++)
2597 if (perm[
ridx (i)] < minr)
2599 minr = perm[
ridx (i)];
2603 double tmp = work[k] /
data (mini);
2606 i <
cidx (k+1); i++)
2612 work[iidx] = work[iidx] - tmp *
data (i);
2620 atmp += fabs (work[i]);
2623 if (atmp > ainvnorm)
2626 rcond = 1. / ainvnorm / anorm;
2632 retval.
resize (nc, b_nc, 0.);
2647 goto triangular_error;
2650 double tmp = work[k] /
data (
cidx (k));
2653 i <
cidx (k+1); i++)
2656 work[iidx] = work[iidx] - tmp *
data (i);
2662 retval.
xelem (i, j) = work[i];
2680 double tmp = work[k] /
data (
cidx (k));
2683 i <
cidx (k+1); i++)
2686 work[iidx] = work[iidx] - tmp *
data (i);
2693 atmp += fabs (work[i]);
2696 if (atmp > ainvnorm)
2699 rcond = 1. / ainvnorm / anorm;
2708 sing_handler (rcond);
2715 volatile double rcond_plus_one = rcond + 1.0;
2723 sing_handler (rcond);
2737 solve_singularity_handler sing_handler,
2738 bool calc_cond)
const
2747 if (nr != b.
rows ())
2749 (
"matrix dimension mismatch solution of linear equations");
2751 if (nr == 0 || nc == 0 || b.
cols () == 0)
2756 int typ = mattype.
type ();
2760 (*current_liboctave_error_handler) (
"incorrect matrix type");
2763 double ainvnorm = 0.;
2773 atmp += fabs (
data (i));
2782 retval.
xcidx (0) = 0;
2796 work[perm[b.
ridx (i)]] = b.
data (i);
2806 if (perm[
ridx (i)] < minr)
2808 minr = perm[
ridx (i)];
2812 if (minr != k ||
data (mini) == 0)
2815 goto triangular_error;
2818 double tmp = work[k] /
data (mini);
2826 work[iidx] = work[iidx] - tmp *
data (i);
2838 if (ii + new_nnz > x_nz)
2849 retval.
xridx (ii) = i;
2850 retval.
xdata (ii++) = work[i];
2852 retval.
xcidx (j+1) = ii;
2875 i <
cidx (k+1); i++)
2876 if (perm[
ridx (i)] < minr)
2878 minr = perm[
ridx (i)];
2882 double tmp = work[k] /
data (mini);
2885 i <
cidx (k+1); i++)
2891 work[iidx] = work[iidx] - tmp *
data (i);
2899 atmp += fabs (work[i]);
2902 if (atmp > ainvnorm)
2905 rcond = 1. / ainvnorm / anorm;
2926 goto triangular_error;
2929 double tmp = work[k] /
data (
cidx (k));
2934 work[iidx] = work[iidx] - tmp *
data (i);
2946 if (ii + new_nnz > x_nz)
2957 retval.
xridx (ii) = i;
2958 retval.
xdata (ii++) = work[i];
2960 retval.
xcidx (j+1) = ii;
2980 double tmp = work[k] /
data (
cidx (k));
2983 i <
cidx (k+1); i++)
2986 work[iidx] = work[iidx] - tmp *
data (i);
2993 atmp += fabs (work[i]);
2996 if (atmp > ainvnorm)
2999 rcond = 1. / ainvnorm / anorm;
3008 sing_handler (rcond);
3015 volatile double rcond_plus_one = rcond + 1.0;
3023 sing_handler (rcond);
3037 solve_singularity_handler sing_handler,
3038 bool calc_cond)
const
3047 if (nr != b.
rows ())
3049 (
"matrix dimension mismatch solution of linear equations");
3051 if (nr == 0 || nc == 0 || b.
cols () == 0)
3056 int typ = mattype.
type ();
3060 (*current_liboctave_error_handler) (
"incorrect matrix type");
3063 double ainvnorm = 0.;
3074 atmp += fabs (
data (i));
3082 retval.
resize (nc, b_nc);
3091 cwork[perm[i]] = b(i, j);
3101 if (perm[
ridx (i)] < minr)
3103 minr = perm[
ridx (i)];
3107 if (minr != k ||
data (mini) == 0)
3110 goto triangular_error;
3121 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3127 retval(i, j) = cwork[i];
3149 i <
cidx (k+1); i++)
3150 if (perm[
ridx (i)] < minr)
3152 minr = perm[
ridx (i)];
3156 double tmp = work[k] /
data (mini);
3159 i <
cidx (k+1); i++)
3165 work[iidx] = work[iidx] - tmp *
data (i);
3173 atmp += fabs (work[i]);
3176 if (atmp > ainvnorm)
3179 rcond = 1. / ainvnorm / anorm;
3185 retval.
resize (nc, b_nc, 0.);
3201 goto triangular_error;
3209 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3215 retval.
xelem (i, j) = cwork[i];
3234 double tmp = work[k] /
data (
cidx (k));
3237 i <
cidx (k+1); i++)
3240 work[iidx] = work[iidx] - tmp *
data (i);
3247 atmp += fabs (work[i]);
3250 if (atmp > ainvnorm)
3253 rcond = 1. / ainvnorm / anorm;
3262 sing_handler (rcond);
3269 volatile double rcond_plus_one = rcond + 1.0;
3277 sing_handler (rcond);
3291 solve_singularity_handler sing_handler,
3292 bool calc_cond)
const
3301 if (nr != b.
rows ())
3303 (
"matrix dimension mismatch solution of linear equations");
3305 if (nr == 0 || nc == 0 || b.
cols () == 0)
3310 int typ = mattype.
type ();
3314 (*current_liboctave_error_handler) (
"incorrect matrix type");
3317 double ainvnorm = 0.;
3327 atmp += fabs (
data (i));
3336 retval.
xcidx (0) = 0;
3350 cwork[perm[b.
ridx (i)]] = b.
data (i);
3360 if (perm[
ridx (i)] < minr)
3362 minr = perm[
ridx (i)];
3366 if (minr != k ||
data (mini) == 0)
3369 goto triangular_error;
3380 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3392 if (ii + new_nnz > x_nz)
3403 retval.
xridx (ii) = i;
3404 retval.
xdata (ii++) = cwork[i];
3406 retval.
xcidx (j+1) = ii;
3430 i <
cidx (k+1); i++)
3431 if (perm[
ridx (i)] < minr)
3433 minr = perm[
ridx (i)];
3437 double tmp = work[k] /
data (mini);
3440 i <
cidx (k+1); i++)
3446 work[iidx] = work[iidx] - tmp *
data (i);
3454 atmp += fabs (work[i]);
3457 if (atmp > ainvnorm)
3460 rcond = 1. / ainvnorm / anorm;
3481 goto triangular_error;
3489 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3501 if (ii + new_nnz > x_nz)
3512 retval.
xridx (ii) = i;
3513 retval.
xdata (ii++) = cwork[i];
3515 retval.
xcidx (j+1) = ii;
3536 double tmp = work[k] /
data (
cidx (k));
3539 i <
cidx (k+1); i++)
3542 work[iidx] = work[iidx] - tmp *
data (i);
3549 atmp += fabs (work[i]);
3552 if (atmp > ainvnorm)
3555 rcond = 1. / ainvnorm / anorm;
3564 sing_handler (rcond);
3571 volatile double rcond_plus_one = rcond + 1.0;
3579 sing_handler (rcond);
3593 solve_singularity_handler sing_handler,
3594 bool calc_cond)
const
3602 if (nr != nc || nr != b.
rows ())
3603 (*current_liboctave_error_handler)
3604 (
"matrix dimension mismatch solution of linear equations");
3606 if (nr == 0 || b.
cols () == 0)
3609 (*current_liboctave_error_handler)
3610 (
"calculation of condition number not implemented");
3614 volatile int typ = mattype.
type ();
3632 D[nc-1] =
data (ii);
3648 else if (
ridx (i) == j + 1)
3653 F77_INT tmp_nr = octave::to_f77_int (nr);
3663 F77_XFCN (dptsv, DPTSV, (tmp_nr, b_nc, D, DL, result,
3691 DL[j] =
data (ii++);
3692 DU[j] =
data (ii++);
3694 D[nc-1] =
data (ii);
3711 else if (
ridx (i) == j + 1)
3713 else if (
ridx (i) == j - 1)
3718 F77_INT tmp_nr = octave::to_f77_int (nr);
3728 F77_XFCN (dgtsv, DGTSV, (tmp_nr, b_nc, DL, D, DU, result,
3740 sing_handler (rcond);
3750 (*current_liboctave_error_handler) (
"incorrect matrix type");
3759 solve_singularity_handler sing_handler,
3760 bool calc_cond)
const
3768 if (nr != nc || nr != b.
rows ())
3769 (*current_liboctave_error_handler)
3770 (
"matrix dimension mismatch solution of linear equations");
3772 if (nr == 0 || b.
cols () == 0)
3775 (*current_liboctave_error_handler)
3776 (
"calculation of condition number not implemented");
3780 int typ = mattype.
type ();
3801 DL[j] =
data (ii++);
3802 DU[j] =
data (ii++);
3804 D[nc-1] =
data (ii);
3821 else if (
ridx (i) == j + 1)
3823 else if (
ridx (i) == j - 1)
3828 F77_INT tmp_nr = octave::to_f77_int (nr);
3832 F77_XFCN (dgttrf, DGTTRF, (tmp_nr, DL, D, DU, DU2, pipvt, tmp_err));
3841 sing_handler (rcond);
3854 retval.
xcidx (0) = 0;
3869 (F77_CONST_CHAR_ARG2 (&job, 1),
3870 tmp_nr, 1, DL, D, DU, DU2, pipvt,
3872 F77_CHAR_ARG_LEN (1)));
3883 if (ii + new_nnz > x_nz)
3894 retval.
xridx (ii) = i;
3895 retval.
xdata (ii++) = work[i];
3897 retval.
xcidx (j+1) = ii;
3904 (*current_liboctave_error_handler) (
"incorrect matrix type");
3913 solve_singularity_handler sing_handler,
3914 bool calc_cond)
const
3922 if (nr != nc || nr != b.
rows ())
3923 (*current_liboctave_error_handler)
3924 (
"matrix dimension mismatch solution of linear equations");
3926 if (nr == 0 || b.
cols () == 0)
3929 (*current_liboctave_error_handler)
3930 (
"calculation of condition number not implemented");
3934 volatile int typ = mattype.
type ();
3952 D[nc-1] =
data (ii);
3968 else if (
ridx (i) == j + 1)
3973 F77_INT tmp_nr = octave::to_f77_int (nr);
4012 DL[j] =
data (ii++);
4013 DU[j] =
data (ii++);
4015 D[nc-1] =
data (ii);
4032 else if (
ridx (i) == j + 1)
4034 else if (
ridx (i) == j - 1)
4039 F77_INT tmp_nr = octave::to_f77_int (nr);
4066 sing_handler (rcond);
4074 (*current_liboctave_error_handler) (
"incorrect matrix type");
4083 solve_singularity_handler sing_handler,
4084 bool calc_cond)
const
4092 if (nr != nc || nr != b.
rows ())
4093 (*current_liboctave_error_handler)
4094 (
"matrix dimension mismatch solution of linear equations");
4096 if (nr == 0 || b.
cols () == 0)
4099 (*current_liboctave_error_handler)
4100 (
"calculation of condition number not implemented");
4104 int typ = mattype.
type ();
4125 DL[j] =
data (ii++);
4126 DU[j] =
data (ii++);
4128 D[nc-1] =
data (ii);
4145 else if (
ridx (i) == j + 1)
4147 else if (
ridx (i) == j - 1)
4152 F77_INT tmp_nr = octave::to_f77_int (nr);
4156 F77_XFCN (dgttrf, DGTTRF, (tmp_nr, DL, D, DU, DU2, pipvt, tmp_err));
4167 sing_handler (rcond);
4188 retval.
xcidx (0) = 0;
4192 for (
F77_INT i = 0; i < b_nr; i++)
4200 (F77_CONST_CHAR_ARG2 (&job, 1),
4201 tmp_nr, 1, DL, D, DU, DU2, pipvt,
4203 F77_CHAR_ARG_LEN (1)));
4210 (*current_liboctave_error_handler)
4211 (
"SparseMatrix::solve solve failed");
4218 (F77_CONST_CHAR_ARG2 (&job, 1),
4219 tmp_nr, 1, DL, D, DU, DU2, pipvt,
4221 F77_CHAR_ARG_LEN (1)));
4228 (*current_liboctave_error_handler)
4229 (
"SparseMatrix::solve solve failed");
4239 if (Bx[i] != 0. || Bz[i] != 0.)
4242 if (ii + new_nnz > x_nz)
4251 if (Bx[i] != 0. || Bz[i] != 0.)
4253 retval.
xridx (ii) = i;
4257 retval.
xcidx (j+1) = ii;
4264 (*current_liboctave_error_handler) (
"incorrect matrix type");
4273 solve_singularity_handler sing_handler,
4274 bool calc_cond)
const
4282 if (nr != nc || nr != b.
rows ())
4283 (*current_liboctave_error_handler)
4284 (
"matrix dimension mismatch solution of linear equations");
4286 if (nr == 0 || b.
cols () == 0)
4291 volatile int typ = mattype.
type ();
4307 tmp_data[ii++] = 0.;
4315 m_band(ri - j, j) =
data (i);
4323 F77_INT tmp_nr = octave::to_f77_int (nr);
4328 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4329 tmp_nr, n_lower, tmp_data, ldm, tmp_err
4330 F77_CHAR_ARG_LEN (1)));
4353 (F77_CONST_CHAR_ARG2 (&job, 1),
4354 tmp_nr, n_lower, tmp_data, ldm,
4355 anorm, rcond, pz, piz, tmp_err
4356 F77_CHAR_ARG_LEN (1)));
4363 volatile double rcond_plus_one = rcond + 1.0;
4371 sing_handler (rcond);
4390 (F77_CONST_CHAR_ARG2 (&job, 1),
4391 tmp_nr, n_lower, b_nc, tmp_data,
4392 ldm, result, b_nr, tmp_err
4393 F77_CHAR_ARG_LEN (1)));
4400 (*current_liboctave_error_handler)
4401 (
"SparseMatrix::solve solve failed");
4413 F77_INT ldm = n_upper + 2 * n_lower + 1;
4422 for (
F77_INT j = 0; j < ldm; j++)
4424 tmp_data[ii++] = 0.;
4429 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
4439 atmp += fabs (
data (i));
4445 F77_INT tmp_nr = octave::to_f77_int (nr);
4452 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4453 tmp_data, ldm, pipvt, tmp_err));
4466 sing_handler (rcond);
4482 F77_INT tmp_nc = octave::to_f77_int (nc);
4485 (F77_CONST_CHAR_ARG2 (&job, 1),
4486 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4487 anorm, rcond, pz, piz, tmp_err
4488 F77_CHAR_ARG_LEN (1)));
4495 volatile double rcond_plus_one = rcond + 1.0;
4503 sing_handler (rcond);
4523 (F77_CONST_CHAR_ARG2 (&job, 1),
4524 tmp_nr, n_lower, n_upper, b_nc, tmp_data,
4525 ldm, pipvt, result, b_nr, tmp_err
4526 F77_CHAR_ARG_LEN (1)));
4533 (*current_liboctave_error_handler) (
"incorrect matrix type");
4542 solve_singularity_handler sing_handler,
4543 bool calc_cond)
const
4551 if (nr != nc || nr != b.
rows ())
4552 (*current_liboctave_error_handler)
4553 (
"matrix dimension mismatch solution of linear equations");
4555 if (nr == 0 || b.
cols () == 0)
4560 volatile int typ = mattype.
type ();
4566 F77_INT ldm = octave::to_f77_int (n_lower + 1);
4575 for (
F77_INT j = 0; j < ldm; j++)
4577 tmp_data[ii++] = 0.;
4585 m_band(ri - j, j) =
data (i);
4593 F77_INT tmp_nr = octave::to_f77_int (nr);
4598 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4599 tmp_nr, n_lower, tmp_data, ldm, tmp_err
4600 F77_CHAR_ARG_LEN (1)));
4621 (F77_CONST_CHAR_ARG2 (&job, 1),
4622 tmp_nr, n_lower, tmp_data, ldm,
4623 anorm, rcond, pz, piz, tmp_err
4624 F77_CHAR_ARG_LEN (1)));
4631 volatile double rcond_plus_one = rcond + 1.0;
4639 sing_handler (rcond);
4661 retval.
xcidx (0) = 0;
4664 for (
F77_INT i = 0; i < b_nr; i++)
4665 Bx[i] = b.
elem (i, j);
4668 (F77_CONST_CHAR_ARG2 (&job, 1),
4669 tmp_nr, n_lower, 1, tmp_data,
4670 ldm, Bx, b_nr, tmp_err
4671 F77_CHAR_ARG_LEN (1)));
4678 (*current_liboctave_error_handler)
4679 (
"SparseMatrix::solve solve failed");
4684 for (
F77_INT i = 0; i < b_nr; i++)
4693 sz = (
static_cast<double> (b_nc) - j) / b_nc
4695 sz = x_nz + (sz > 100 ? sz : 100);
4699 retval.
xdata (ii) = tmp;
4700 retval.
xridx (ii++) = i;
4703 retval.
xcidx (j+1) = ii;
4716 F77_INT ldm = octave::to_f77_int (n_upper + 2 * n_lower + 1);
4727 tmp_data[ii++] = 0.;
4732 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
4743 atmp += fabs (
data (i));
4749 F77_INT tmp_nr = octave::to_f77_int (nr);
4756 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4757 tmp_data, ldm, pipvt, tmp_err));
4768 sing_handler (rcond);
4784 F77_INT tmp_nc = octave::to_f77_int (nc);
4787 (F77_CONST_CHAR_ARG2 (&job, 1),
4788 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4789 anorm, rcond, pz, piz, tmp_err
4790 F77_CHAR_ARG_LEN (1)));
4797 volatile double rcond_plus_one = rcond + 1.0;
4805 sing_handler (rcond);
4821 retval.
xcidx (0) = 0;
4831 i < b.
cidx (j+1); i++)
4837 (F77_CONST_CHAR_ARG2 (&job, 1),
4838 tmp_nr, n_lower, n_upper, 1, tmp_data,
4839 ldm, pipvt, work, b_nr, tmp_err
4840 F77_CHAR_ARG_LEN (1)));
4851 if (ii + new_nnz > x_nz)
4862 retval.
xridx (ii) = i;
4863 retval.
xdata (ii++) = work[i];
4865 retval.
xcidx (j+1) = ii;
4873 (*current_liboctave_error_handler) (
"incorrect matrix type");
4882 solve_singularity_handler sing_handler,
4883 bool calc_cond)
const
4891 if (nr != nc || nr != b.
rows ())
4892 (*current_liboctave_error_handler)
4893 (
"matrix dimension mismatch solution of linear equations");
4895 if (nr == 0 || b.
cols () == 0)
4900 volatile int typ = mattype.
type ();
4915 for (
F77_INT j = 0; j < ldm; j++)
4917 tmp_data[ii++] = 0.;
4925 m_band(ri - j, j) =
data (i);
4933 F77_INT tmp_nr = octave::to_f77_int (nr);
4938 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4939 tmp_nr, n_lower, tmp_data, ldm, tmp_err
4940 F77_CHAR_ARG_LEN (1)));
4963 (F77_CONST_CHAR_ARG2 (&job, 1),
4964 tmp_nr, n_lower, tmp_data, ldm,
4965 anorm, rcond, pz, piz, tmp_err
4966 F77_CHAR_ARG_LEN (1)));
4973 volatile double rcond_plus_one = rcond + 1.0;
4981 sing_handler (rcond);
4999 retval.
resize (b_nr, b_nc);
5003 for (
F77_INT i = 0; i < b_nr; i++)
5011 (F77_CONST_CHAR_ARG2 (&job, 1),
5012 tmp_nr, n_lower, 1, tmp_data,
5013 ldm, Bx, b_nr, tmp_err
5014 F77_CHAR_ARG_LEN (1)));
5021 (*current_liboctave_error_handler)
5022 (
"SparseMatrix::solve solve failed");
5028 (F77_CONST_CHAR_ARG2 (&job, 1),
5029 tmp_nr, n_lower, 1, tmp_data,
5030 ldm, Bz, b_nr, tmp_err
5031 F77_CHAR_ARG_LEN (1)));
5038 (*current_liboctave_error_handler)
5039 (
"SparseMatrix::solve solve failed");
5045 retval(i, j) =
Complex (Bx[i], Bz[i]);
5056 F77_INT ldm = n_upper + 2 * n_lower + 1;
5065 for (
F77_INT j = 0; j < ldm; j++)
5067 tmp_data[ii++] = 0.;
5072 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
5083 atmp += fabs (
data (i));
5089 F77_INT tmp_nr = octave::to_f77_int (nr);
5096 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5097 tmp_data, ldm, pipvt, tmp_err));
5108 sing_handler (rcond);
5124 F77_INT tmp_nc = octave::to_f77_int (nc);
5127 (F77_CONST_CHAR_ARG2 (&job, 1),
5128 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5129 anorm, rcond, pz, piz, tmp_err
5130 F77_CHAR_ARG_LEN (1)));
5137 volatile double rcond_plus_one = rcond + 1.0;
5145 sing_handler (rcond);
5160 retval.
resize (nr, b_nc);
5175 (F77_CONST_CHAR_ARG2 (&job, 1),
5176 tmp_nr, n_lower, n_upper, 1, tmp_data,
5177 ldm, pipvt, Bx, b_nr, tmp_err
5178 F77_CHAR_ARG_LEN (1)));
5183 (F77_CONST_CHAR_ARG2 (&job, 1),
5184 tmp_nr, n_lower, n_upper, 1, tmp_data,
5185 ldm, pipvt, Bz, b_nr, tmp_err
5186 F77_CHAR_ARG_LEN (1)));
5191 retval(i, j) =
Complex (Bx[i], Bz[i]);
5197 (*current_liboctave_error_handler) (
"incorrect matrix type");
5206 solve_singularity_handler sing_handler,
5207 bool calc_cond)
const
5215 if (nr != nc || nr != b.
rows ())
5216 (*current_liboctave_error_handler)
5217 (
"matrix dimension mismatch solution of linear equations");
5219 if (nr == 0 || b.
cols () == 0)
5224 volatile int typ = mattype.
type ();
5239 for (
F77_INT j = 0; j < ldm; j++)
5241 tmp_data[ii++] = 0.;
5249 m_band(ri - j, j) =
data (i);
5257 F77_INT tmp_nr = octave::to_f77_int (nr);
5262 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5263 tmp_nr, n_lower, tmp_data, ldm, tmp_err
5264 F77_CHAR_ARG_LEN (1)));
5288 (F77_CONST_CHAR_ARG2 (&job, 1),
5289 tmp_nr, n_lower, tmp_data, ldm,
5290 anorm, rcond, pz, piz, tmp_err
5291 F77_CHAR_ARG_LEN (1)));
5298 volatile double rcond_plus_one = rcond + 1.0;
5306 sing_handler (rcond);
5329 retval.
xcidx (0) = 0;
5333 for (
F77_INT i = 0; i < b_nr; i++)
5341 (F77_CONST_CHAR_ARG2 (&job, 1),
5342 tmp_nr, n_lower, 1, tmp_data,
5343 ldm, Bx, b_nr, tmp_err
5344 F77_CHAR_ARG_LEN (1)));
5351 (*current_liboctave_error_handler)
5352 (
"SparseMatrix::solve solve failed");
5358 (F77_CONST_CHAR_ARG2 (&job, 1),
5359 tmp_nr, n_lower, 1, tmp_data,
5360 ldm, Bz, b_nr, tmp_err
5361 F77_CHAR_ARG_LEN (1)));
5368 (*current_liboctave_error_handler)
5369 (
"SparseMatrix::solve solve failed");
5379 if (Bx[i] != 0. || Bz[i] != 0.)
5382 if (ii + new_nnz > x_nz)
5391 if (Bx[i] != 0. || Bz[i] != 0.)
5393 retval.
xridx (ii) = i;
5397 retval.
xcidx (j+1) = ii;
5410 F77_INT ldm = n_upper + 2 * n_lower + 1;
5419 for (
F77_INT j = 0; j < ldm; j++)
5421 tmp_data[ii++] = 0.;
5426 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
5437 atmp += fabs (
data (i));
5443 F77_INT tmp_nr = octave::to_f77_int (nr);
5450 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5451 tmp_data, ldm, pipvt, tmp_err));
5462 sing_handler (rcond);
5478 F77_INT tmp_nc = octave::to_f77_int (nc);
5481 (F77_CONST_CHAR_ARG2 (&job, 1),
5482 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5483 anorm, rcond, pz, piz, tmp_err
5484 F77_CHAR_ARG_LEN (1)));
5491 volatile double rcond_plus_one = rcond + 1.0;
5499 sing_handler (rcond);
5516 retval.
xcidx (0) = 0;
5530 i < b.
cidx (j+1); i++)
5533 Bx[b.
ridx (i)] = c.real ();
5534 Bz[b.
ridx (i)] = c.imag ();
5538 (F77_CONST_CHAR_ARG2 (&job, 1),
5539 tmp_nr, n_lower, n_upper, 1, tmp_data,
5540 ldm, pipvt, Bx, b_nr, tmp_err
5541 F77_CHAR_ARG_LEN (1)));
5546 (F77_CONST_CHAR_ARG2 (&job, 1),
5547 tmp_nr, n_lower, n_upper, 1, tmp_data,
5548 ldm, pipvt, Bz, b_nr, tmp_err
5549 F77_CHAR_ARG_LEN (1)));
5557 if (Bx[i] != 0. || Bz[i] != 0.)
5560 if (ii + new_nnz > x_nz)
5569 if (Bx[i] != 0. || Bz[i] != 0.)
5571 retval.
xridx (ii) = i;
5574 retval.
xcidx (j+1) = ii;
5582 (*current_liboctave_error_handler) (
"incorrect matrix type");
5590 Matrix& Info, solve_singularity_handler sing_handler,
5591 bool calc_cond)
const
5594 void *Numeric =
nullptr;
5597 #if defined (HAVE_UMFPACK)
5600 Control =
Matrix (UMFPACK_CONTROL, 1);
5604 double tmp = octave::sparse_params::get_key (
"spumoni");
5606 Control (UMFPACK_PRL) = tmp;
5607 tmp = octave::sparse_params::get_key (
"piv_tol");
5610 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
5611 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
5615 tmp = octave::sparse_params::get_key (
"autoamd");
5617 Control (UMFPACK_FIXQ) = tmp;
5623 const double *Ax =
data ();
5633 Info =
Matrix (1, UMFPACK_INFO);
5638 Ax,
nullptr, &Symbolic, control, info);
5648 (*current_liboctave_error_handler)
5649 (
"SparseMatrix::solve symbolic factorization failed");
5658 Ax, Symbolic, &Numeric, control, info);
5662 rcond = Info (UMFPACK_RCOND);
5665 volatile double rcond_plus_one = rcond + 1.0;
5667 if (status == UMFPACK_WARNING_singular_matrix
5675 sing_handler (rcond);
5679 else if (status < 0)
5685 (*current_liboctave_error_handler)
5686 (
"SparseMatrix::solve numeric factorization failed");
5701 octave_unused_parameter (rcond);
5702 octave_unused_parameter (Control);
5703 octave_unused_parameter (Info);
5704 octave_unused_parameter (sing_handler);
5705 octave_unused_parameter (calc_cond);
5707 (*current_liboctave_error_handler)
5708 (
"support for UMFPACK was unavailable or disabled "
5709 "when liboctave was built");
5719 solve_singularity_handler sing_handler,
5720 bool calc_cond)
const
5728 if (nr != nc || nr != b.
rows ())
5729 (*current_liboctave_error_handler)
5730 (
"matrix dimension mismatch solution of linear equations");
5732 if (nr == 0 || b.
cols () == 0)
5737 volatile int typ = mattype.
type ();
5742 #if defined (HAVE_CHOLMOD)
5743 cholmod_common Common;
5744 cholmod_common *cm = &Common;
5748 cm->prefer_zomplex =
false;
5750 double spu = octave::sparse_params::get_key (
"spumoni");
5754 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
5758 cm->print =
static_cast<int> (spu) + 2;
5759 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5763 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5764 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5766 cm->final_ll =
true;
5768 cholmod_sparse Astore;
5769 cholmod_sparse *
A = &Astore;
5779 #if defined (OCTAVE_ENABLE_64)
5780 A->itype = CHOLMOD_LONG;
5782 A->itype = CHOLMOD_INT;
5784 A->dtype = CHOLMOD_DOUBLE;
5786 A->xtype = CHOLMOD_REAL;
5790 cholmod_dense Bstore;
5791 cholmod_dense *
B = &Bstore;
5792 B->nrow = b.
rows ();
5793 B->ncol = b.
cols ();
5795 B->nzmax =
B->nrow *
B->ncol;
5796 B->dtype = CHOLMOD_DOUBLE;
5797 B->xtype = CHOLMOD_REAL;
5799 B->x =
const_cast<double *
> (b.
data ());
5816 volatile double rcond_plus_one = rcond + 1.0;
5824 sing_handler (rcond);
5840 retval.
xelem (i, j) =
static_cast<double *
>(X->x)[jr + i];
5846 static char blank_name[] =
" ";
5850 (*current_liboctave_warning_with_id_handler)
5851 (
"Octave:missing-dependency",
5852 "support for CHOLMOD was unavailable or disabled "
5853 "when liboctave was built");
5862 #if defined (HAVE_UMFPACK)
5865 =
factorize (err, rcond, Control, Info, sing_handler, calc_cond);
5870 Control (UMFPACK_IRSTEP) = 1;
5871 const double *Bx = b.
data ();
5881 const double *Ax =
data ();
5888 Ax, &result[iidx], &Bx[iidx],
5889 Numeric, control, info);
5895 (*current_liboctave_error_handler)
5896 (
"SparseMatrix::solve solve failed");
5911 octave_unused_parameter (rcond);
5912 octave_unused_parameter (sing_handler);
5913 octave_unused_parameter (calc_cond);
5915 (*current_liboctave_error_handler)
5916 (
"support for UMFPACK was unavailable or disabled "
5917 "when liboctave was built");
5921 (*current_liboctave_error_handler) (
"incorrect matrix type");
5930 solve_singularity_handler sing_handler,
5931 bool calc_cond)
const
5939 if (nr != nc || nr != b.
rows ())
5940 (*current_liboctave_error_handler)
5941 (
"matrix dimension mismatch solution of linear equations");
5943 if (nr == 0 || b.
cols () == 0)
5948 volatile int typ = mattype.
type ();
5953 #if defined (HAVE_CHOLMOD)
5954 cholmod_common Common;
5955 cholmod_common *cm = &Common;
5959 cm->prefer_zomplex =
false;
5961 double spu = octave::sparse_params::get_key (
"spumoni");
5965 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
5969 cm->print =
static_cast<int> (spu) + 2;
5970 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5974 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5975 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5977 cm->final_ll =
true;
5979 cholmod_sparse Astore;
5980 cholmod_sparse *
A = &Astore;
5990 #if defined (OCTAVE_ENABLE_64)
5991 A->itype = CHOLMOD_LONG;
5993 A->itype = CHOLMOD_INT;
5995 A->dtype = CHOLMOD_DOUBLE;
5997 A->xtype = CHOLMOD_REAL;
6001 cholmod_sparse Bstore;
6002 cholmod_sparse *
B = &Bstore;
6003 B->nrow = b.
rows ();
6004 B->ncol = b.
cols ();
6007 B->nzmax = b.
nnz ();
6011 #if defined (OCTAVE_ENABLE_64)
6012 B->itype = CHOLMOD_LONG;
6014 B->itype = CHOLMOD_INT;
6016 B->dtype = CHOLMOD_DOUBLE;
6018 B->xtype = CHOLMOD_REAL;
6037 volatile double rcond_plus_one = rcond + 1.0;
6045 sing_handler (rcond);
6054 cholmod_sparse *X =
CHOLMOD_NAME(spsolve) (CHOLMOD_A, L,
B, cm);
6060 j <= static_cast<octave_idx_type> (X->ncol); j++)
6063 j < static_cast<octave_idx_type> (X->nzmax); j++)
6066 retval.
xdata (j) =
static_cast<double *
>(X->x)[j];
6072 static char blank_name[] =
" ";
6076 (*current_liboctave_warning_with_id_handler)
6077 (
"Octave:missing-dependency",
6078 "support for CHOLMOD was unavailable or disabled "
6079 "when liboctave was built");
6088 #if defined (HAVE_UMFPACK)
6090 void *Numeric =
factorize (err, rcond, Control, Info,
6091 sing_handler, calc_cond);
6096 Control (UMFPACK_IRSTEP) = 1;
6104 const double *Ax =
data ();
6115 retval.
xcidx (0) = 0;
6120 Bx[i] = b.
elem (i, j);
6125 Ax, Xx, Bx, Numeric,
6132 (*current_liboctave_error_handler)
6133 (
"SparseMatrix::solve solve failed");
6148 sz = (
static_cast<double> (b_nc) - j) / b_nc
6150 sz = x_nz + (sz > 100 ? sz : 100);
6154 retval.
xdata (ii) = tmp;
6155 retval.
xridx (ii++) = i;
6158 retval.
xcidx (j+1) = ii;
6171 octave_unused_parameter (rcond);
6172 octave_unused_parameter (sing_handler);
6173 octave_unused_parameter (calc_cond);
6175 (*current_liboctave_error_handler)
6176 (
"support for UMFPACK was unavailable or disabled "
6177 "when liboctave was built");
6181 (*current_liboctave_error_handler) (
"incorrect matrix type");
6190 solve_singularity_handler sing_handler,
6191 bool calc_cond)
const
6199 if (nr != nc || nr != b.
rows ())
6200 (*current_liboctave_error_handler)
6201 (
"matrix dimension mismatch solution of linear equations");
6203 if (nr == 0 || b.
cols () == 0)
6208 volatile int typ = mattype.
type ();
6213 #if defined (HAVE_CHOLMOD)
6214 cholmod_common Common;
6215 cholmod_common *cm = &Common;
6219 cm->prefer_zomplex =
false;
6221 double spu = octave::sparse_params::get_key (
"spumoni");
6225 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6229 cm->print =
static_cast<int> (spu) + 2;
6230 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6234 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6235 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6237 cm->final_ll =
true;
6239 cholmod_sparse Astore;
6240 cholmod_sparse *
A = &Astore;
6250 #if defined (OCTAVE_ENABLE_64)
6251 A->itype = CHOLMOD_LONG;
6253 A->itype = CHOLMOD_INT;
6255 A->dtype = CHOLMOD_DOUBLE;
6257 A->xtype = CHOLMOD_REAL;
6261 cholmod_dense Bstore;
6262 cholmod_dense *
B = &Bstore;
6263 B->nrow = b.
rows ();
6264 B->ncol = b.
cols ();
6266 B->nzmax =
B->nrow *
B->ncol;
6267 B->dtype = CHOLMOD_DOUBLE;
6268 B->xtype = CHOLMOD_COMPLEX;
6287 volatile double rcond_plus_one = rcond + 1.0;
6295 sing_handler (rcond);
6311 retval.
xelem (i, j) =
static_cast<Complex *
>(X->x)[jr + i];
6317 static char blank_name[] =
" ";
6321 (*current_liboctave_warning_with_id_handler)
6322 (
"Octave:missing-dependency",
6323 "support for CHOLMOD was unavailable or disabled "
6324 "when liboctave was built");
6333 #if defined (HAVE_UMFPACK)
6335 void *Numeric =
factorize (err, rcond, Control, Info,
6336 sing_handler, calc_cond);
6341 Control (UMFPACK_IRSTEP) = 1;
6349 const double *Ax =
data ();
6354 retval.
resize (b_nr, b_nc);
6371 Ax, Xx, Bx, Numeric,
6377 Numeric, control, info);
6379 if (status < 0 || status2 < 0)
6384 (*current_liboctave_error_handler)
6385 (
"SparseMatrix::solve solve failed");
6392 retval(i, j) =
Complex (Xx[i], Xz[i]);
6403 octave_unused_parameter (rcond);
6404 octave_unused_parameter (sing_handler);
6405 octave_unused_parameter (calc_cond);
6407 (*current_liboctave_error_handler)
6408 (
"support for UMFPACK was unavailable or disabled "
6409 "when liboctave was built");
6413 (*current_liboctave_error_handler) (
"incorrect matrix type");
6422 solve_singularity_handler sing_handler,
6423 bool calc_cond)
const
6431 if (nr != nc || nr != b.
rows ())
6432 (*current_liboctave_error_handler)
6433 (
"matrix dimension mismatch solution of linear equations");
6435 if (nr == 0 || b.
cols () == 0)
6440 volatile int typ = mattype.
type ();
6445 #if defined (HAVE_CHOLMOD)
6446 cholmod_common Common;
6447 cholmod_common *cm = &Common;
6451 cm->prefer_zomplex =
false;
6453 double spu = octave::sparse_params::get_key (
"spumoni");
6457 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6461 cm->print =
static_cast<int> (spu) + 2;
6462 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6466 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6467 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6469 cm->final_ll =
true;
6471 cholmod_sparse Astore;
6472 cholmod_sparse *
A = &Astore;
6482 #if defined (OCTAVE_ENABLE_64)
6483 A->itype = CHOLMOD_LONG;
6485 A->itype = CHOLMOD_INT;
6487 A->dtype = CHOLMOD_DOUBLE;
6489 A->xtype = CHOLMOD_REAL;
6493 cholmod_sparse Bstore;
6494 cholmod_sparse *
B = &Bstore;
6495 B->nrow = b.
rows ();
6496 B->ncol = b.
cols ();
6499 B->nzmax = b.
nnz ();
6503 #if defined (OCTAVE_ENABLE_64)
6504 B->itype = CHOLMOD_LONG;
6506 B->itype = CHOLMOD_INT;
6508 B->dtype = CHOLMOD_DOUBLE;
6510 B->xtype = CHOLMOD_COMPLEX;
6529 volatile double rcond_plus_one = rcond + 1.0;
6537 sing_handler (rcond);
6546 cholmod_sparse *X =
CHOLMOD_NAME(spsolve) (CHOLMOD_A, L,
B, cm);
6553 j <= static_cast<octave_idx_type> (X->ncol); j++)
6556 j < static_cast<octave_idx_type> (X->nzmax); j++)
6565 static char blank_name[] =
" ";
6569 (*current_liboctave_warning_with_id_handler)
6570 (
"Octave:missing-dependency",
6571 "support for CHOLMOD was unavailable or disabled "
6572 "when liboctave was built");
6581 #if defined (HAVE_UMFPACK)
6583 void *Numeric =
factorize (err, rcond, Control, Info,
6584 sing_handler, calc_cond);
6589 Control (UMFPACK_IRSTEP) = 1;
6597 const double *Ax =
data ();
6611 retval.
xcidx (0) = 0;
6624 Ax, Xx, Bx, Numeric,
6630 Numeric, control, info);
6632 if (status < 0 || status2 < 0)
6637 (*current_liboctave_error_handler)
6638 (
"SparseMatrix::solve solve failed");
6653 sz = (
static_cast<double> (b_nc) - j) / b_nc
6655 sz = x_nz + (sz > 100 ? sz : 100);
6659 retval.
xdata (ii) = tmp;
6660 retval.
xridx (ii++) = i;
6663 retval.
xcidx (j+1) = ii;
6675 octave_unused_parameter (rcond);
6676 octave_unused_parameter (sing_handler);
6677 octave_unused_parameter (calc_cond);
6679 (*current_liboctave_error_handler)
6680 (
"support for UMFPACK was unavailable or disabled "
6681 "when liboctave was built");
6685 (*current_liboctave_error_handler) (
"incorrect matrix type");
6696 return solve (mattype, b, info, rcond,
nullptr);
6704 return solve (mattype, b, info, rcond,
nullptr);
6711 return solve (mattype, b, info, rcond,
nullptr);
6716 double& rcond, solve_singularity_handler sing_handler,
6717 bool singular_fallback)
const
6720 int typ = mattype.
type (
false);
6723 typ = mattype.
type (*
this);
6727 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6729 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6731 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6733 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6736 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6738 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6740 (*current_liboctave_error_handler) (
"unknown matrix type");
6746 #if defined (USE_QRSOLVE)
6747 retval =
qrsolve (*
this, b, err);
6761 return solve (mattype, b, info, rcond,
nullptr);
6769 return solve (mattype, b, info, rcond,
nullptr);
6776 return solve (mattype, b, info, rcond,
nullptr);
6782 solve_singularity_handler sing_handler,
6783 bool singular_fallback)
const
6786 int typ = mattype.
type (
false);
6789 typ = mattype.
type (*
this);
6792 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6794 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6796 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6798 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6801 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6803 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6805 (*current_liboctave_error_handler) (
"unknown matrix type");
6810 #if defined (USE_QRSOLVE)
6811 retval =
qrsolve (*
this, b, err);
6826 return solve (mattype, b, info, rcond,
nullptr);
6834 return solve (mattype, b, info, rcond,
nullptr);
6841 return solve (mattype, b, info, rcond,
nullptr);
6847 solve_singularity_handler sing_handler,
6848 bool singular_fallback)
const
6851 int typ = mattype.
type (
false);
6854 typ = mattype.
type (*
this);
6857 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6859 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6861 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6863 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6866 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6868 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6870 (*current_liboctave_error_handler) (
"unknown matrix type");
6875 #if defined (USE_QRSOLVE)
6876 retval =
qrsolve (*
this, b, err);
6891 return solve (mattype, b, info, rcond,
nullptr);
6899 return solve (mattype, b, info, rcond,
nullptr);
6906 return solve (mattype, b, info, rcond,
nullptr);
6912 solve_singularity_handler sing_handler,
6913 bool singular_fallback)
const
6916 int typ = mattype.
type (
false);
6919 typ = mattype.
type (*
this);
6922 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6924 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6926 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6928 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6931 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6933 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6935 (*current_liboctave_error_handler) (
"unknown matrix type");
6940 #if defined (USE_QRSOLVE)
6941 retval =
qrsolve (*
this, b, err);
6955 return solve (mattype, b, info, rcond);
6963 return solve (mattype, b, info, rcond);
6970 return solve (mattype, b, info, rcond,
nullptr);
6976 solve_singularity_handler sing_handler)
const
6979 return solve (mattype, tmp, info, rcond,
6988 return solve (mattype, b, info, rcond,
nullptr);
6996 return solve (mattype, b, info, rcond,
nullptr);
7002 double& rcond)
const
7004 return solve (mattype, b, info, rcond,
nullptr);
7010 solve_singularity_handler sing_handler)
const
7013 return solve (mattype, tmp, info, rcond,
7022 return solve (b, info, rcond,
nullptr);
7029 return solve (b, info, rcond,
nullptr);
7034 double& rcond)
const
7036 return solve (b, info, rcond,
nullptr);
7041 solve_singularity_handler sing_handler)
const
7044 return solve (mattype, b, err, rcond, sing_handler);
7052 return solve (b, info, rcond,
nullptr);
7060 return solve (b, info, rcond,
nullptr);
7067 return solve (b, info, rcond,
nullptr);
7072 solve_singularity_handler sing_handler)
const
7075 return solve (mattype, b, err, rcond, sing_handler);
7082 return solve (b, info, rcond,
nullptr);
7087 double& rcond)
const
7089 return solve (b, info, rcond,
nullptr);
7095 solve_singularity_handler sing_handler)
const
7098 return solve (mattype, b, err, rcond, sing_handler);
7106 return solve (b, info, rcond,
nullptr);
7113 return solve (b, info, rcond,
nullptr);
7118 double& rcond)
const
7120 return solve (b, info, rcond,
nullptr);
7126 solve_singularity_handler sing_handler)
const
7129 return solve (mattype, b, err, rcond, sing_handler);
7136 return solve (b, info, rcond);
7143 return solve (b, info, rcond);
7148 double& rcond)
const
7150 return solve (b, info, rcond,
nullptr);
7156 solve_singularity_handler sing_handler)
const
7159 return solve (tmp, info, rcond,
7168 return solve (b, info, rcond,
nullptr);
7175 return solve (b, info, rcond,
nullptr);
7180 double& rcond)
const
7182 return solve (b, info, rcond,
nullptr);
7188 solve_singularity_handler sing_handler)
const
7191 return solve (tmp, info, rcond,
7225 double val =
data (i);
7240 double val =
data (i);
7255 double val =
data (i);
7256 if (val != 0.0 && val != 1.0)
7282 double val =
data (i);
7308 double val =
data (i);
7349 if (jj <
cidx (i+1) &&
ridx (jj) == j)
7392 if ((
rows () == 1 && dim == -1) || dim == 1)
7397 (
cidx (j+1) -
cidx (j) < nr ? 0.0 : 1.0), 1.0);
7411 double d = data (i); \
7412 tmp[ridx (i)] += d * d
7415 double d = data (i); \
7433 retval.
data (i) = fabs (retval.
data (i));
7462 os << a.
ridx (i) + 1 <<
' ' << j + 1 <<
' ';
7463 octave::write_value<double> (os, a.
data (i));
7476 return read_sparse_matrix<elt_type> (is, a, octave::read_value<double>);
7540 return do_mul_dm_sm<SparseMatrix> (
d, a);
7546 return do_mul_sm_dm<SparseMatrix> (a,
d);
7552 return do_add_dm_sm<SparseMatrix> (
d, a);
7558 return do_sub_dm_sm<SparseMatrix> (
d, a);
7564 return do_add_sm_dm<SparseMatrix> (a,
d);
7570 return do_sub_sm_dm<SparseMatrix> (a,
d);
7589 #define EMPTY_RETURN_CHECK(T) \
7590 if (nr == 0 || nc == 0) \
7614 result.
xdata (idx) = tmp;
7615 result.
xridx (idx) =
m.ridx (i);
7630 result.
xcidx (0) = 0;
7639 result.
xdata (ii) = tmp;
7640 result.
xridx (ii++) =
m.ridx (i);
7643 result.
xcidx (j+1) = ii;
7666 if (a_nr == b_nr && a_nc == b_nc)
7676 bool ja_lt_max = ja < ja_max;
7680 bool jb_lt_max = jb < jb_max;
7682 while (ja_lt_max || jb_lt_max)
7685 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
7690 r.ridx (jx) = a.
ridx (ja);
7695 ja_lt_max= ja < ja_max;
7697 else if ((! ja_lt_max)
7698 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
7703 r.ridx (jx) = b.
ridx (jb);
7708 jb_lt_max= jb < jb_max;
7716 r.ridx (jx) = a.
ridx (ja);
7720 ja_lt_max= ja < ja_max;
7722 jb_lt_max= jb < jb_max;
7728 r.maybe_compress ();
7732 if (a_nr == 0 || a_nc == 0)
7733 r.resize (a_nr, a_nc);
7734 else if (b_nr == 0 || b_nc == 0)
7735 r.resize (b_nr, b_nc);
7765 result.
xdata (idx) = tmp;
7766 result.
xridx (idx) =
m.ridx (i);
7781 result.
xcidx (0) = 0;
7789 result.
xdata (ii) = tmp;
7790 result.
xridx (ii++) =
m.ridx (i);
7793 result.
xcidx (j+1) = ii;
7816 if (a_nr == b_nr && a_nc == b_nc)
7826 bool ja_lt_max = ja < ja_max;
7830 bool jb_lt_max = jb < jb_max;
7832 while (ja_lt_max || jb_lt_max)
7835 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
7840 r.ridx (jx) = a.
ridx (ja);
7845 ja_lt_max= ja < ja_max;
7847 else if ((! ja_lt_max)
7848 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
7853 r.ridx (jx) = b.
ridx (jb);
7858 jb_lt_max= jb < jb_max;
7866 r.ridx (jx) = a.
ridx (ja);
7870 ja_lt_max= ja < ja_max;
7872 jb_lt_max= jb < jb_max;
7878 r.maybe_compress ();
7882 if (a_nr == 0 || a_nc == 0)
7883 r.resize (a_nr, a_nc);
7884 else if (b_nr == 0 || b_nc == 0)
7885 r.resize (b_nr, b_nc);
#define SPARSE_SSM_BOOL_OPS(S, M)
#define SPARSE_ANY_OP(DIM)
#define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, INIT_VAL, MT_RESULT)
#define SPARSE_FULL_MUL(RET_TYPE, EL_TYPE)
#define FULL_SPARSE_MUL_TRANS(RET_TYPE, EL_TYPE, CONJ_OP)
#define SPARSE_SSM_CMP_OPS(S, M)
#define SPARSE_SMSM_CMP_OPS(M1, M2)
#define SPARSE_REDUCTION_OP(RET_TYPE, EL_TYPE, OP, INIT_VAL, MT_RESULT)
#define SPARSE_SMS_CMP_OPS(M, S)
#define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE)
#define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN)
#define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN)
#define SPARSE_FULL_TRANS_MUL(RET_TYPE, EL_TYPE, CONJ_OP)
#define FULL_SPARSE_MUL(RET_TYPE, EL_TYPE)
#define SPARSE_SMSM_BOOL_OPS(M1, M2)
#define SPARSE_SMS_BOOL_OPS(M, S)
#define SPARSE_ALL_OP(DIM)
SM octinternal_do_mul_pm_sm(const PermMatrix &p, const SM &a)
SM octinternal_do_mul_sm_pm(const SM &a, const PermMatrix &p)
OCTARRAY_OVERRIDABLE_FUNC_API const T * data(void) const
Size of the specified dimension.
OCTARRAY_API void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type rows(void) const
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
OCTARRAY_OVERRIDABLE_FUNC_API T & elem(octave_idx_type n)
Size of the specified dimension.
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type cols(void) const
OCTARRAY_OVERRIDABLE_FUNC_API T & xelem(octave_idx_type n)
Size of the specified dimension.
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
octave_idx_type length(void) const
octave_idx_type cols(void) const
MSparse< T > permute(const Array< octave_idx_type > &vec, bool inv=false) const
MSparse< T > squeeze(void) const
MSparse< T > ipermute(const Array< octave_idx_type > &vec) const
MSparse< T > reshape(const dim_vector &new_dims) const
MSparse< T > & insert(const Sparse< T > &a, octave_idx_type r, octave_idx_type c)
MSparse< T > diag(octave_idx_type k=0) const
bool ishermitian(void) const
OCTAVE_API MatrixType transpose(void) const
void mark_as_rectangular(void)
bool is_dense(void) const
OCTAVE_API void mark_as_unsymmetric(void)
octave_idx_type * triangular_perm(void) const
OCTAVE_API void info(void) const
OCTAVE_API int type(bool quiet=true)
OCTAVE_API Matrix sum(int dim=-1) const
OCTAVE_API RowVector row(octave_idx_type i) const
OCTAVE_API Matrix abs(void) const
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
OCTAVE_API ColumnVector column(octave_idx_type i) const
OCTAVE_API double max(void) const
OCTAVE_API SparseComplexMatrix & insert(const SparseComplexMatrix &a, octave_idx_type r, octave_idx_type c)
OCTAVE_API bool issymmetric(void) const
OCTAVE_API SparseMatrix dinverse(MatrixType &mattype, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
OCTAVE_API SparseMatrix & insert(const SparseMatrix &a, octave_idx_type r, octave_idx_type c)
OCTAVE_API bool any_element_is_inf_or_nan(void) const
OCTAVE_API Matrix ltsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
OCTAVE_API SparseMatrix reshape(const dim_vector &new_dims) const
OCTAVE_API SparseMatrix sum(int dim=-1) const
OCTAVE_API SparseBoolMatrix any(int dim=-1) const
OCTAVE_API Matrix utsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
OCTAVE_API Matrix fsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
OCTAVE_API bool all_integers(double &max_val, double &min_val) const
OCTAVE_API bool operator==(const SparseMatrix &a) const
OCTAVE_API SparseMatrix abs(void) const
OCTAVE_API bool operator!=(const SparseMatrix &a) const
OCTAVE_API RowVector row(octave_idx_type i) const
OCTAVE_API SparseMatrix cumprod(int dim=-1) const
OCTAVE_API SparseMatrix diag(octave_idx_type k=0) const
OCTAVE_API SparseMatrix prod(int dim=-1) const
OCTAVE_API SparseBoolMatrix all(int dim=-1) const
OCTAVE_API Matrix matrix_value(void) const
OCTAVE_API SparseMatrix min(int dim=-1) const
OCTAVE_API bool any_element_is_nan(void) const
OCTAVE_API SparseMatrix permute(const Array< octave_idx_type > &vec, bool inv=false) const
OCTAVE_API SparseMatrix cumsum(int dim=-1) const
OCTAVE_API Matrix trisolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
OCTAVE_API ColumnVector column(octave_idx_type i) const
SparseMatrix transpose(void) const
OCTAVE_API SparseMatrix inverse(void) const
OCTAVE_API SparseMatrix concat(const SparseMatrix &rb, const Array< octave_idx_type > &ra_idx)
OCTAVE_API bool too_large_for_float(void) const
OCTAVE_API SparseMatrix ipermute(const Array< octave_idx_type > &vec) const
OCTAVE_API bool any_element_not_one_or_zero(void) const
OCTAVE_API DET determinant(void) const
OCTAVE_API bool all_elements_are_zero(void) const
OCTAVE_API SparseMatrix tinverse(MatrixType &mattype, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
OCTAVE_API SparseMatrix sumsq(int dim=-1) const
OCTAVE_API SparseMatrix squeeze(void) const
OCTAVE_API Matrix solve(MatrixType &typ, const Matrix &b) const
OCTAVE_API Matrix bsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
OCTAVE_API bool all_elements_are_int_or_inf_or_nan(void) const
OCTAVE_API Matrix dsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
OCTAVE_API SparseBoolMatrix operator!(void) const
OCTAVE_API SparseMatrix max(int dim=-1) const
OCTAVE_API bool any_element_is_negative(bool=false) const
OCTAVE_API void * factorize(octave_idx_type &err, double &rcond, Matrix &Control, Matrix &Info, solve_singularity_handler sing_handler, bool calc_cond=false) const
octave_idx_type rows(void) const
octave_idx_type * cidx(void)
octave_idx_type columns(void) const
octave_idx_type * ridx(void)
octave_idx_type nnz(void) const
Actual number of nonzero terms.
OCTAVE_API Array< T > array_value(void) const
dim_vector dims(void) const
bool test_any(F fcn) const
Sparse< T, Alloc > maybe_compress(bool remove_zeros=false)
T & elem(octave_idx_type n)
octave_idx_type cols(void) const
void change_capacity(octave_idx_type nz)
octave_idx_type * xcidx(void)
octave_idx_type * xridx(void)
Vector representing the dimensions (size) of an Array.
int first_non_singleton(int def=0) const
octave_idx_type ndims(void) const
Number of dimensions.
Matrix trans_mul(const SparseMatrix &m, const Matrix &a)
SparseMatrix real(const SparseComplexMatrix &a)
#define EMPTY_RETURN_CHECK(T)
SparseMatrix min(double d, const SparseMatrix &m)
Matrix mul_trans(const Matrix &m, const SparseMatrix &a)
SparseMatrix operator-(const DiagMatrix &d, const SparseMatrix &a)
SparseMatrix max(double d, const SparseMatrix &m)
std::istream & operator>>(std::istream &is, SparseMatrix &a)
SparseMatrix operator*(const SparseMatrix &m, const SparseMatrix &a)
SparseMatrix imag(const SparseComplexMatrix &a)
SparseMatrix operator+(const DiagMatrix &d, const SparseMatrix &a)
std::ostream & operator<<(std::ostream &os, const SparseMatrix &a)
#define F77_DBLE_CMPLX_ARG(x)
#define F77_XFCN(f, F, args)
octave_f77_int_type F77_INT
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
void err_nan_to_logical_conversion(void)
void warn_singular_matrix(double rcond)
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
#define lo_ieee_signbit(x)
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
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
F77_RET_T const F77_INT F77_CMPLX * A
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
bool too_large_for_float(double x)
class OCTAVE_API SparseMatrix
class OCTAVE_API SparseComplexMatrix
std::complex< double > Complex
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
#define CHOLMOD_NAME(name)
#define UMFPACK_DNAME(name)
const octave_base_value const Array< octave_idx_type > & ra_idx
template class OCTAVE_API sparse_chol< SparseMatrix >
template OCTAVE_API SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseMatrix, SparseComplexMatrix >(const SparseMatrix &, const SparseComplexMatrix &, octave_idx_type &)
template OCTAVE_API ComplexMatrix dmsolve< ComplexMatrix, SparseMatrix, ComplexMatrix >(const SparseMatrix &, const ComplexMatrix &, octave_idx_type &)
template OCTAVE_API Matrix dmsolve< Matrix, SparseMatrix, Matrix >(const SparseMatrix &, const Matrix &, octave_idx_type &)
template OCTAVE_API SparseMatrix dmsolve< SparseMatrix, SparseMatrix, SparseMatrix >(const SparseMatrix &, const SparseMatrix &, octave_idx_type &)
template class OCTAVE_API sparse_lu< SparseMatrix >
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
void SparseCholError(int status, char *file, int line, char *message)
int SparseCholPrint(const char *fmt,...)