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);
946 ret = tinverse (mattype, info, rcond,
true, calc_cond).
transpose ();
950 ret =
transpose ().tinverse (newtype, info, rcond,
true, calc_cond);
957 octave::math::sparse_chol<SparseMatrix> fact (*
this, info,
false);
958 rcond = fact.rcond ();
964 info, rcond2,
true,
false);
965 ret =
Q * InvL.
transpose () * InvL *
Q.transpose ();
982 octave::math::sparse_lu<SparseMatrix> fact (*
this,
985 rcond = fact.rcond ();
1001 info, rcond2,
true,
false);
1002 SparseMatrix InvU = fact.U ().tinverse (tmp_typ, info, rcond2,
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 ();
3792 F77_INT *pipvt = ipvt.fortran_vec ();
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 ();
4116 F77_INT *pipvt = ipvt.fortran_vec ();
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 ();
4299 double *tmp_data = m_band.fortran_vec ();
4307 tmp_data[ii++] = 0.;
4315 m_band(ri - j, j) =
data (i);
4321 anorm = m_band.abs ().sum ().row (0).max ();
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)));
4348 double *pz = z.fortran_vec ();
4350 F77_INT *piz = iz.fortran_vec ();
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;
4416 double *tmp_data = m_band.fortran_vec ();
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);
4448 F77_INT *pipvt = ipvt.fortran_vec ();
4452 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4453 tmp_data, ldm, pipvt, tmp_err));
4466 sing_handler (rcond);
4478 double *pz = z.fortran_vec ();
4480 F77_INT *piz = iz.fortran_vec ();
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);
4569 double *tmp_data = m_band.fortran_vec ();
4575 for (
F77_INT j = 0; j < ldm; j++)
4577 tmp_data[ii++] = 0.;
4585 m_band(ri - j, j) =
data (i);
4591 anorm = m_band.abs ().sum ().row (0).max ();
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)));
4616 double *pz = z.fortran_vec ();
4618 F77_INT *piz = iz.fortran_vec ();
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);
4719 double *tmp_data = m_band.fortran_vec ();
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);
4752 F77_INT *pipvt = ipvt.fortran_vec ();
4756 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4757 tmp_data, ldm, pipvt, tmp_err));
4768 sing_handler (rcond);
4780 double *pz = z.fortran_vec ();
4782 F77_INT *piz = iz.fortran_vec ();
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 ();
4909 double *tmp_data = m_band.fortran_vec ();
4915 for (
F77_INT j = 0; j < ldm; j++)
4917 tmp_data[ii++] = 0.;
4925 m_band(ri - j, j) =
data (i);
4931 anorm = m_band.abs ().sum ().row (0).max ();
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)));
4958 double *pz = z.fortran_vec ();
4960 F77_INT *piz = iz.fortran_vec ();
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;
5059 double *tmp_data = m_band.fortran_vec ();
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);
5092 F77_INT *pipvt = ipvt.fortran_vec ();
5096 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5097 tmp_data, ldm, pipvt, tmp_err));
5108 sing_handler (rcond);
5120 double *pz = z.fortran_vec ();
5122 F77_INT *piz = iz.fortran_vec ();
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 ();
5233 double *tmp_data = m_band.fortran_vec ();
5239 for (
F77_INT j = 0; j < ldm; j++)
5241 tmp_data[ii++] = 0.;
5249 m_band(ri - j, j) =
data (i);
5255 anorm = m_band.abs ().sum ().row (0).max ();
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)));
5283 double *pz = z.fortran_vec ();
5285 F77_INT *piz = iz.fortran_vec ();
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;
5413 double *tmp_data = m_band.fortran_vec ();
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);
5446 F77_INT *pipvt = ipvt.fortran_vec ();
5450 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5451 tmp_data, ldm, pipvt, tmp_err));
5462 sing_handler (rcond);
5474 double *pz = z.fortran_vec ();
5476 F77_INT *piz = iz.fortran_vec ();
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);
6062 j <= static_cast<octave_idx_type> (X->ncol); j++)
6069 retval.
xdata (j) =
static_cast<double *
> (X->x)[j];
6075 static char blank_name[] =
" ";
6079 (*current_liboctave_warning_with_id_handler)
6080 (
"Octave:missing-dependency",
6081 "support for CHOLMOD was unavailable or disabled "
6082 "when liboctave was built");
6091 #if defined (HAVE_UMFPACK)
6093 void *Numeric = factorize (err, rcond, Control, Info,
6094 sing_handler, calc_cond);
6099 Control (UMFPACK_IRSTEP) = 1;
6107 const double *Ax =
data ();
6118 retval.
xcidx (0) = 0;
6123 Bx[i] = b.
elem (i, j);
6128 Ax, Xx, Bx, Numeric,
6135 (*current_liboctave_error_handler)
6136 (
"SparseMatrix::solve solve failed");
6151 sz = (
static_cast<double> (b_nc) - j) / b_nc
6153 sz = x_nz + (sz > 100 ? sz : 100);
6157 retval.
xdata (ii) = tmp;
6158 retval.
xridx (ii++) = i;
6161 retval.
xcidx (j+1) = ii;
6174 octave_unused_parameter (rcond);
6175 octave_unused_parameter (sing_handler);
6176 octave_unused_parameter (calc_cond);
6178 (*current_liboctave_error_handler)
6179 (
"support for UMFPACK was unavailable or disabled "
6180 "when liboctave was built");
6184 (*current_liboctave_error_handler) (
"incorrect matrix type");
6193 solve_singularity_handler sing_handler,
6194 bool calc_cond)
const
6202 if (nr != nc || nr != b.
rows ())
6203 (*current_liboctave_error_handler)
6204 (
"matrix dimension mismatch solution of linear equations");
6206 if (nr == 0 || b.
cols () == 0)
6211 volatile int typ = mattype.
type ();
6216 #if defined (HAVE_CHOLMOD)
6217 cholmod_common Common;
6218 cholmod_common *cm = &Common;
6222 cm->prefer_zomplex =
false;
6224 double spu = octave::sparse_params::get_key (
"spumoni");
6228 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6232 cm->print =
static_cast<int> (spu) + 2;
6233 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6237 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6238 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6240 cm->final_ll =
true;
6242 cholmod_sparse Astore;
6243 cholmod_sparse *
A = &Astore;
6253 #if defined (OCTAVE_ENABLE_64)
6254 A->itype = CHOLMOD_LONG;
6256 A->itype = CHOLMOD_INT;
6258 A->dtype = CHOLMOD_DOUBLE;
6260 A->xtype = CHOLMOD_REAL;
6264 cholmod_dense Bstore;
6265 cholmod_dense *
B = &Bstore;
6266 B->nrow = b.
rows ();
6267 B->ncol = b.
cols ();
6269 B->nzmax =
B->nrow *
B->ncol;
6270 B->dtype = CHOLMOD_DOUBLE;
6271 B->xtype = CHOLMOD_COMPLEX;
6290 volatile double rcond_plus_one = rcond + 1.0;
6298 sing_handler (rcond);
6314 retval.
xelem (i, j) =
static_cast<Complex *
> (X->x)[jr + i];
6320 static char blank_name[] =
" ";
6324 (*current_liboctave_warning_with_id_handler)
6325 (
"Octave:missing-dependency",
6326 "support for CHOLMOD was unavailable or disabled "
6327 "when liboctave was built");
6336 #if defined (HAVE_UMFPACK)
6338 void *Numeric = factorize (err, rcond, Control, Info,
6339 sing_handler, calc_cond);
6344 Control (UMFPACK_IRSTEP) = 1;
6352 const double *Ax =
data ();
6357 retval.
resize (b_nr, b_nc);
6374 Ax, Xx, Bx, Numeric,
6380 Numeric, control, info);
6382 if (status < 0 || status2 < 0)
6387 (*current_liboctave_error_handler)
6388 (
"SparseMatrix::solve solve failed");
6395 retval(i, j) =
Complex (Xx[i], Xz[i]);
6406 octave_unused_parameter (rcond);
6407 octave_unused_parameter (sing_handler);
6408 octave_unused_parameter (calc_cond);
6410 (*current_liboctave_error_handler)
6411 (
"support for UMFPACK was unavailable or disabled "
6412 "when liboctave was built");
6416 (*current_liboctave_error_handler) (
"incorrect matrix type");
6425 solve_singularity_handler sing_handler,
6426 bool calc_cond)
const
6434 if (nr != nc || nr != b.
rows ())
6435 (*current_liboctave_error_handler)
6436 (
"matrix dimension mismatch solution of linear equations");
6438 if (nr == 0 || b.
cols () == 0)
6443 volatile int typ = mattype.
type ();
6448 #if defined (HAVE_CHOLMOD)
6449 cholmod_common Common;
6450 cholmod_common *cm = &Common;
6454 cm->prefer_zomplex =
false;
6456 double spu = octave::sparse_params::get_key (
"spumoni");
6460 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6464 cm->print =
static_cast<int> (spu) + 2;
6465 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6469 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6470 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6472 cm->final_ll =
true;
6474 cholmod_sparse Astore;
6475 cholmod_sparse *
A = &Astore;
6485 #if defined (OCTAVE_ENABLE_64)
6486 A->itype = CHOLMOD_LONG;
6488 A->itype = CHOLMOD_INT;
6490 A->dtype = CHOLMOD_DOUBLE;
6492 A->xtype = CHOLMOD_REAL;
6496 cholmod_sparse Bstore;
6497 cholmod_sparse *
B = &Bstore;
6498 B->nrow = b.
rows ();
6499 B->ncol = b.
cols ();
6502 B->nzmax = b.
nnz ();
6506 #if defined (OCTAVE_ENABLE_64)
6507 B->itype = CHOLMOD_LONG;
6509 B->itype = CHOLMOD_INT;
6511 B->dtype = CHOLMOD_DOUBLE;
6513 B->xtype = CHOLMOD_COMPLEX;
6532 volatile double rcond_plus_one = rcond + 1.0;
6540 sing_handler (rcond);
6549 cholmod_sparse *X =
CHOLMOD_NAME(spsolve) (CHOLMOD_A, L,
B, cm);
6557 j <= static_cast<octave_idx_type> (X->ncol); j++)
6570 static char blank_name[] =
" ";
6574 (*current_liboctave_warning_with_id_handler)
6575 (
"Octave:missing-dependency",
6576 "support for CHOLMOD was unavailable or disabled "
6577 "when liboctave was built");
6586 #if defined (HAVE_UMFPACK)
6588 void *Numeric = factorize (err, rcond, Control, Info,
6589 sing_handler, calc_cond);
6594 Control (UMFPACK_IRSTEP) = 1;
6602 const double *Ax =
data ();
6616 retval.
xcidx (0) = 0;
6629 Ax, Xx, Bx, Numeric,
6635 Numeric, control, info);
6637 if (status < 0 || status2 < 0)
6642 (*current_liboctave_error_handler)
6643 (
"SparseMatrix::solve solve failed");
6658 sz = (
static_cast<double> (b_nc) - j) / b_nc
6660 sz = x_nz + (sz > 100 ? sz : 100);
6664 retval.
xdata (ii) = tmp;
6665 retval.
xridx (ii++) = i;
6668 retval.
xcidx (j+1) = ii;
6680 octave_unused_parameter (rcond);
6681 octave_unused_parameter (sing_handler);
6682 octave_unused_parameter (calc_cond);
6684 (*current_liboctave_error_handler)
6685 (
"support for UMFPACK was unavailable or disabled "
6686 "when liboctave was built");
6690 (*current_liboctave_error_handler) (
"incorrect matrix type");
6701 return solve (mattype, b, info, rcond,
nullptr);
6709 return solve (mattype, b, info, rcond,
nullptr);
6716 return solve (mattype, b, info, rcond,
nullptr);
6721 double& rcond, solve_singularity_handler sing_handler,
6722 bool singular_fallback)
const
6725 int typ = mattype.
type (
false);
6728 typ = mattype.
type (*
this);
6732 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
6734 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
6736 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
6738 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
6741 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
6743 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
6745 (*current_liboctave_error_handler) (
"unknown matrix type");
6751 #if defined (USE_QRSOLVE)
6752 retval =
qrsolve (*
this, b, err);
6766 return solve (mattype, b, info, rcond,
nullptr);
6774 return solve (mattype, b, info, rcond,
nullptr);
6781 return solve (mattype, b, info, rcond,
nullptr);
6787 solve_singularity_handler sing_handler,
6788 bool singular_fallback)
const
6791 int typ = mattype.
type (
false);
6794 typ = mattype.
type (*
this);
6797 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
6799 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
6801 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
6803 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
6806 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
6808 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
6810 (*current_liboctave_error_handler) (
"unknown matrix type");
6815 #if defined (USE_QRSOLVE)
6816 retval =
qrsolve (*
this, b, err);
6831 return solve (mattype, b, info, rcond,
nullptr);
6839 return solve (mattype, b, info, rcond,
nullptr);
6846 return solve (mattype, b, info, rcond,
nullptr);
6852 solve_singularity_handler sing_handler,
6853 bool singular_fallback)
const
6856 int typ = mattype.
type (
false);
6859 typ = mattype.
type (*
this);
6862 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
6864 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
6866 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
6868 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
6871 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
6873 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
6875 (*current_liboctave_error_handler) (
"unknown matrix type");
6880 #if defined (USE_QRSOLVE)
6881 retval =
qrsolve (*
this, b, err);
6896 return solve (mattype, b, info, rcond,
nullptr);
6904 return solve (mattype, b, info, rcond,
nullptr);
6911 return solve (mattype, b, info, rcond,
nullptr);
6917 solve_singularity_handler sing_handler,
6918 bool singular_fallback)
const
6921 int typ = mattype.
type (
false);
6924 typ = mattype.
type (*
this);
6927 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
6929 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
6931 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
6933 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
6936 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
6938 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
6940 (*current_liboctave_error_handler) (
"unknown matrix type");
6945 #if defined (USE_QRSOLVE)
6946 retval =
qrsolve (*
this, b, err);
6960 return solve (mattype, b, info, rcond);
6968 return solve (mattype, b, info, rcond);
6975 return solve (mattype, b, info, rcond,
nullptr);
6981 solve_singularity_handler sing_handler)
const
6984 return solve (mattype, tmp, info, rcond,
6993 return solve (mattype, b, info, rcond,
nullptr);
7001 return solve (mattype, b, info, rcond,
nullptr);
7007 double& rcond)
const
7009 return solve (mattype, b, info, rcond,
nullptr);
7015 solve_singularity_handler sing_handler)
const
7018 return solve (mattype, tmp, info, rcond,
7027 return solve (b, info, rcond,
nullptr);
7034 return solve (b, info, rcond,
nullptr);
7039 double& rcond)
const
7041 return solve (b, info, rcond,
nullptr);
7046 solve_singularity_handler sing_handler)
const
7049 return solve (mattype, b, err, rcond, sing_handler);
7057 return solve (b, info, rcond,
nullptr);
7065 return solve (b, info, rcond,
nullptr);
7072 return solve (b, info, rcond,
nullptr);
7077 solve_singularity_handler sing_handler)
const
7080 return solve (mattype, b, err, rcond, sing_handler);
7087 return solve (b, info, rcond,
nullptr);
7092 double& rcond)
const
7094 return solve (b, info, rcond,
nullptr);
7100 solve_singularity_handler sing_handler)
const
7103 return solve (mattype, b, err, rcond, sing_handler);
7111 return solve (b, info, rcond,
nullptr);
7118 return solve (b, info, rcond,
nullptr);
7123 double& rcond)
const
7125 return solve (b, info, rcond,
nullptr);
7131 solve_singularity_handler sing_handler)
const
7134 return solve (mattype, b, err, rcond, sing_handler);
7141 return solve (b, info, rcond);
7148 return solve (b, info, rcond);
7153 double& rcond)
const
7155 return solve (b, info, rcond,
nullptr);
7161 solve_singularity_handler sing_handler)
const
7164 return solve (tmp, info, rcond,
7173 return solve (b, info, rcond,
nullptr);
7180 return solve (b, info, rcond,
nullptr);
7185 double& rcond)
const
7187 return solve (b, info, rcond,
nullptr);
7193 solve_singularity_handler sing_handler)
const
7196 return solve (tmp, info, rcond,
7230 double val =
data (i);
7245 double val =
data (i);
7260 double val =
data (i);
7261 if (val != 0.0 && val != 1.0)
7287 double val =
data (i);
7313 double val =
data (i);
7354 if (jj <
cidx (i+1) &&
ridx (jj) == j)
7397 if ((
rows () == 1 && dim == -1) || dim == 1)
7402 (
cidx (j+1) -
cidx (j) < nr ? 0.0 : 1.0), 1.0);
7416 double d = data (i); \
7417 tmp[ridx (i)] += d * d
7420 double d = data (i); \
7438 retval.
data (i) = fabs (retval.
data (i));
7467 os << a.
ridx (i) + 1 <<
' ' << j + 1 <<
' ';
7468 octave::write_value<double> (os, a.
data (i));
7481 return read_sparse_matrix<elt_type> (is, a, octave::read_value<double>);
7545 return do_mul_dm_sm<SparseMatrix> (
d, a);
7551 return do_mul_sm_dm<SparseMatrix> (a,
d);
7557 return do_add_dm_sm<SparseMatrix> (
d, a);
7563 return do_sub_dm_sm<SparseMatrix> (
d, a);
7569 return do_add_sm_dm<SparseMatrix> (a,
d);
7575 return do_sub_sm_dm<SparseMatrix> (a,
d);
7594 #define EMPTY_RETURN_CHECK(T) \
7595 if (nr == 0 || nc == 0) \
7619 result.
xdata (idx) = tmp;
7620 result.
xridx (idx) =
m.ridx (i);
7635 result.
xcidx (0) = 0;
7644 result.
xdata (ii) = tmp;
7645 result.
xridx (ii++) =
m.ridx (i);
7648 result.
xcidx (j+1) = ii;
7671 if (a_nr == b_nr && a_nc == b_nc)
7681 bool ja_lt_max = ja < ja_max;
7685 bool jb_lt_max = jb < jb_max;
7687 while (ja_lt_max || jb_lt_max)
7690 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
7695 r.ridx (jx) = a.
ridx (ja);
7700 ja_lt_max= ja < ja_max;
7702 else if ((! ja_lt_max)
7703 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
7708 r.ridx (jx) = b.
ridx (jb);
7713 jb_lt_max= jb < jb_max;
7721 r.ridx (jx) = a.
ridx (ja);
7725 ja_lt_max= ja < ja_max;
7727 jb_lt_max= jb < jb_max;
7733 r.maybe_compress ();
7737 if (a_nr == 0 || a_nc == 0)
7738 r.resize (a_nr, a_nc);
7739 else if (b_nr == 0 || b_nc == 0)
7740 r.resize (b_nr, b_nc);
7770 result.
xdata (idx) = tmp;
7771 result.
xridx (idx) =
m.ridx (i);
7786 result.
xcidx (0) = 0;
7794 result.
xdata (ii) = tmp;
7795 result.
xridx (ii++) =
m.ridx (i);
7798 result.
xcidx (j+1) = ii;
7821 if (a_nr == b_nr && a_nc == b_nc)
7831 bool ja_lt_max = ja < ja_max;
7835 bool jb_lt_max = jb < jb_max;
7837 while (ja_lt_max || jb_lt_max)
7840 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
7845 r.ridx (jx) = a.
ridx (ja);
7850 ja_lt_max= ja < ja_max;
7852 else if ((! ja_lt_max)
7853 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
7858 r.ridx (jx) = b.
ridx (jb);
7863 jb_lt_max= jb < jb_max;
7871 r.ridx (jx) = a.
ridx (ja);
7875 ja_lt_max= ja < ja_max;
7877 jb_lt_max= jb < jb_max;
7883 r.maybe_compress ();
7887 if (a_nr == 0 || a_nc == 0)
7888 r.resize (a_nr, a_nc);
7889 else if (b_nr == 0 || b_nc == 0)
7890 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)
T & elem(octave_idx_type n)
Size of the specified dimension.
T * fortran_vec()
Size of the specified dimension.
octave_idx_type rows() const
void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
const T * data() const
Size of the specified dimension.
octave_idx_type cols() const
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() const
octave_idx_type cols() const
MSparse< T > permute(const Array< octave_idx_type > &vec, bool inv=false) const
MSparse< T > squeeze() 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
void mark_as_unsymmetric()
octave_idx_type * triangular_perm() const
MatrixType transpose() const
int type(bool quiet=true)
void mark_as_rectangular()
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
ColumnVector column(octave_idx_type i) const
SparseComplexMatrix & insert(const SparseComplexMatrix &a, octave_idx_type r, octave_idx_type c)
SparseMatrix & insert(const SparseMatrix &a, octave_idx_type r, octave_idx_type c)
SparseMatrix squeeze() const
bool any_element_not_one_or_zero() const
SparseMatrix reshape(const dim_vector &new_dims) const
SparseMatrix sum(int dim=-1) const
SparseBoolMatrix any(int dim=-1) const
bool any_element_is_inf_or_nan() const
bool all_integers(double &max_val, double &min_val) const
bool any_element_is_nan() const
bool too_large_for_float() const
bool operator==(const SparseMatrix &a) const
SparseMatrix transpose() const
SparseMatrix inverse() const
bool operator!=(const SparseMatrix &a) const
RowVector row(octave_idx_type i) const
SparseMatrix cumprod(int dim=-1) const
SparseMatrix diag(octave_idx_type k=0) const
SparseMatrix prod(int dim=-1) const
bool all_elements_are_zero() const
SparseBoolMatrix all(int dim=-1) const
SparseMatrix min(int dim=-1) const
SparseMatrix permute(const Array< octave_idx_type > &vec, bool inv=false) const
SparseMatrix cumsum(int dim=-1) const
ColumnVector column(octave_idx_type i) const
SparseMatrix concat(const SparseMatrix &rb, const Array< octave_idx_type > &ra_idx)
SparseMatrix ipermute(const Array< octave_idx_type > &vec) const
SparseBoolMatrix operator!() const
SparseMatrix sumsq(int dim=-1) const
bool all_elements_are_int_or_inf_or_nan() const
Matrix solve(MatrixType &typ, const Matrix &b) const
SparseMatrix max(int dim=-1) const
bool any_element_is_negative(bool=false) const
Matrix matrix_value() const
octave_idx_type cols() const
octave_idx_type * xcidx()
Array< T > array_value() const
octave_idx_type columns() 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 nnz() const
Actual number of nonzero terms.
octave_idx_type rows() const
void change_capacity(octave_idx_type nz)
octave_idx_type * xridx()
Vector representing the dimensions (size) of an Array.
octave_idx_type ndims() const
Number of dimensions.
int first_non_singleton(int def=0) const
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 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)
std::complex< double > Complex
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
#define CHOLMOD_NAME(name)
octave_idx_type from_suitesparse_long(SuiteSparse_long x)
#define UMFPACK_DNAME(name)
octave_idx_type from_size_t(std::size_t x)
const octave_base_value const Array< octave_idx_type > & ra_idx
template SparseMatrix dmsolve< SparseMatrix, SparseMatrix, SparseMatrix >(const SparseMatrix &, const SparseMatrix &, octave_idx_type &)
template Matrix dmsolve< Matrix, SparseMatrix, Matrix >(const SparseMatrix &, const Matrix &, octave_idx_type &)
template ComplexMatrix dmsolve< ComplexMatrix, SparseMatrix, ComplexMatrix >(const SparseMatrix &, const ComplexMatrix &, octave_idx_type &)
template SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseMatrix, SparseComplexMatrix >(const SparseMatrix &, const SparseComplexMatrix &, octave_idx_type &)
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,...)