26#if defined (HAVE_CONFIG_H)
119 if (nr != nr_a || nc != nc_a || nz != nz_a)
136 return !(*
this == a);
145 if (nr == nc && nr > 0)
180Complex_NaN_result (octave::numeric_limits<double>::NaN (),
181 octave::numeric_limits<double>::NaN ());
187 return max (dummy_idx, dim, nanflag, realabs);
192 bool nanflag,
bool realabs)
const
199 if (dim >= dv.
ndims ())
212 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
221 double real_max = octave::numeric_limits<double>::NaN ();
222 double imag_max = octave::numeric_limits<double>::NaN ();
227 if (
ridx (i) != idx_j)
243 if (octave::math::isnan (tmp) && ! nanflag)
249 else if (octave::math::isnan (tmp))
252 double real_tmp = tmp.real ();
253 double imag_tmp = tmp.imag ();
255 if (octave::math::isnan (real_max) || real_tmp > real_max)
262 else if (real_tmp == real_max && imag_tmp > imag_max)
271 if (octave::math::isnan (tmp_max) && ! nanflag)
272 idx_arg.
elem (j) = idx_j;
274 idx_arg.
elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
275 if (real_max != 0. || imag_max != 0.)
284 double abs_max = octave::numeric_limits<double>::NaN ();
285 double arg_max = octave::numeric_limits<double>::NaN ();
289 if (
ridx (i) != idx_j)
305 if (octave::math::isnan (tmp) && ! nanflag)
311 else if (octave::math::isnan (tmp))
314 double abs_tmp = std::abs (tmp);
315 double arg_tmp = std::arg (tmp);
317 if (octave::math::isnan (abs_max) || abs_tmp > abs_max)
324 else if (abs_tmp == abs_max && arg_tmp > arg_max)
333 if (octave::math::isnan (tmp_max) && ! nanflag)
334 idx_arg.
elem (j) = idx_j;
336 idx_arg.
elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
337 if (abs_max != 0. || arg_max != 0.)
345 result.
xcidx (0) = 0;
351 result.
xdata (ii) = tmp;
352 result.
xridx (ii++) = 0;
354 result.
xcidx (j+1) = ii;
361 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
371 if (found[
ridx (i)] == -j)
372 found[
ridx (i)] = -j - 1;
375 if (found[i] > -nc && found[i] < 0)
376 idx_arg.
elem (i) = -found[i];
388 if (octave::math::isnan (tmp) && ! nanflag)
390 idx_arg.
elem (ir) = j;
393 else if (octave::math::isnan (tmp))
395 else if (ix == -1 || octave::math::isnan (
elem (ir, ix)) ||
396 tmp.real () >
elem (ir, ix).
real ())
397 idx_arg.
elem (ir) = j;
398 else if (tmp.real () ==
elem (ir, ix).
real () &&
400 idx_arg.
elem (ir) = j;
414 if (octave::math::isnan (tmp) && ! nanflag)
416 idx_arg.
elem (ir) = j;
419 else if (octave::math::isnan (tmp))
421 else if (ix == -1 || octave::math::isnan (
elem (ir, ix)) ||
422 std::abs (tmp) > std::abs (
elem (ir, ix)))
423 idx_arg.
elem (ir) = j;
424 else if (std::abs (tmp) == std::abs (
elem (ir, ix)) &&
425 std::arg (tmp) > std::arg (
elem (ir, ix)))
426 idx_arg.
elem (ir) = j;
433 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
439 result.
xcidx (0) = 0;
440 result.
xcidx (1) = nel;
443 if (idx_arg(j) == -1)
446 result.
xdata (ii) = Complex_NaN_result;
447 result.
xridx (ii++) = j;
454 result.
xdata (ii) = tmp;
455 result.
xridx (ii++) = j;
468 return min (dummy_idx, dim, nanflag, realabs);
473 bool nanflag,
bool realabs)
const
480 if (dim >= dv.
ndims ())
493 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
502 double real_min = octave::numeric_limits<double>::NaN ();
503 double imag_min = octave::numeric_limits<double>::NaN ();
508 if (
ridx (i) != idx_j)
524 if (octave::math::isnan (tmp) && ! nanflag)
530 else if (octave::math::isnan (tmp))
533 double real_tmp = tmp.real ();
534 double imag_tmp = tmp.imag ();
536 if (octave::math::isnan (real_min) || real_tmp < real_min)
543 else if (real_tmp == real_min && imag_tmp < imag_min)
552 if (octave::math::isnan (tmp_min) && ! nanflag)
553 idx_arg.
elem (j) = idx_j;
555 idx_arg.
elem (j) = (octave::math::isnan (tmp_min) ? 0 : idx_j);
556 if (real_min != 0. || imag_min != 0.)
565 double abs_min = octave::numeric_limits<double>::NaN ();
566 double arg_min = octave::numeric_limits<double>::NaN ();
571 if (
ridx (i) != idx_j)
587 if (octave::math::isnan (tmp) && ! nanflag)
593 else if (octave::math::isnan (tmp))
596 double abs_tmp = std::abs (tmp);
597 double arg_tmp = std::arg (tmp);
599 if (octave::math::isnan (abs_min) || abs_tmp < abs_min)
606 else if (abs_tmp == abs_min && arg_tmp < arg_min)
615 if (octave::math::isnan (tmp_min) && ! nanflag)
616 idx_arg.
elem (j) = idx_j;
618 idx_arg.
elem (j) = (octave::math::isnan (tmp_min) ? 0 : idx_j);
619 if (abs_min != 0. || arg_min != 0.)
627 result.
xcidx (0) = 0;
633 result.
xdata (ii) = tmp;
634 result.
xridx (ii++) = 0;
636 result.
xcidx (j+1) = ii;
643 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
653 if (found[
ridx (i)] == -j)
654 found[
ridx (i)] = -j - 1;
657 if (found[i] > -nc && found[i] < 0)
658 idx_arg.
elem (i) = -found[i];
670 if (octave::math::isnan (tmp) && ! nanflag)
672 idx_arg.
elem (ir) = j;
675 else if (octave::math::isnan (tmp))
677 else if (ix == -1 || octave::math::isnan (
elem (ir, ix)) ||
678 tmp.real () <
elem (ir, ix).
real ())
679 idx_arg.
elem (ir) = j;
680 else if (tmp.real () ==
elem (ir, ix).
real () &&
682 idx_arg.
elem (ir) = j;
696 if (octave::math::isnan (tmp) && ! nanflag)
698 idx_arg.
elem (ir) = j;
701 else if (octave::math::isnan (tmp))
703 else if (ix == -1 || octave::math::isnan (
elem (ir, ix)) ||
704 std::abs (tmp) < std::abs (
elem (ir, ix)))
705 idx_arg.
elem (ir) = j;
706 else if (std::abs (tmp) == std::abs (
elem (ir, ix)) &&
707 std::arg (tmp) < std::arg (
elem (ir, ix)))
708 idx_arg.
elem (ir) = j;
715 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
721 result.
xcidx (0) = 0;
722 result.
xcidx (1) = nel;
725 if (idx_arg(j) == -1)
728 result.
xdata (ii) = Complex_NaN_result;
729 result.
xridx (ii++) = j;
736 result.
xdata (ii) = tmp;
737 result.
xridx (ii++) = j;
772 retval(j) =
data (k);
799 return insert (tmp , r, c);
815 return insert (tmp , indx);
831 if (rb.
rows () > 0 && rb.
cols () > 0)
841 if (rb.
rows () > 0 && rb.
cols () > 0)
867 retval.
xcidx (i) = nz;
876 retval.
xridx (q) = j;
909is_singular (
const double rcond)
911 return (std::abs (rcond) <= std::numeric_limits<double>::epsilon ());
920 return inverse (mattype, info, rcond, 0, 0);
928 return inverse (mattype, info, rcond, 0, 0);
935 return inverse (mattype, info, rcond, 0, 0);
940 double& rcond,
const bool,
941 const bool calccond)
const
949 if (nr == 0 || nc == 0 || nr != nc)
950 (*current_liboctave_error_handler) (
"inverse requires square matrix");
953 int typ = mattype.
type ();
957 (*current_liboctave_error_handler) (
"incorrect matrix type");
970 double dmin = octave::numeric_limits<double>::Inf ();
973 double tmp = std::abs (v[i]);
990 double& rcond,
const bool,
991 const bool calccond)
const
999 if (nr == 0 || nc == 0 || nr != nc)
1000 (*current_liboctave_error_handler) (
"inverse requires square matrix");
1003 int typ = mattype.
type ();
1008 (*current_liboctave_error_handler) (
"incorrect matrix type");
1011 double ainvnorm = 0.;
1020 atmp += std::abs (
data (i));
1045 retval.
xcidx (i) = cx;
1046 retval.
xridx (cx) = i;
1047 retval.
xdata (cx) = 1.0;
1060 (*current_liboctave_error_handler) (
"division by zero");
1065 rpX = retval.
xridx (colXp);
1074 v -= retval.
xdata (colXp) *
data (colUp);
1079 while (rpX < j && rpU < j && colXp < cx && colUp < nz);
1083 colUp =
cidx (j+1) - 1;
1087 if (pivot == 0. ||
ridx (colUp) != j)
1088 (*current_liboctave_error_handler) (
"division by zero");
1098 retval.
xridx (cx) = j;
1099 retval.
xdata (cx) = v / pivot;
1107 colUp =
cidx (i+1) - 1;
1111 if (pivot == 0. ||
ridx (colUp) != i)
1112 (*current_liboctave_error_handler) (
"division by zero");
1116 retval.
xdata (j) /= pivot;
1118 retval.
xcidx (nr) = cx;
1163 k <
cidx (jidx+1); k++)
1176 (*current_liboctave_error_handler) (
"division by zero");
1178 work[j] = v / pivot;
1184 colUp =
cidx (perm[iidx]+1) - 1;
1186 colUp =
cidx (perm[iidx]);
1190 (*current_liboctave_error_handler) (
"division by zero");
1203 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
1207 retval.
xcidx (i) = cx;
1211 retval.
xridx (cx) = j;
1212 retval.
xdata (cx++) = work[j];
1216 retval.
xcidx (nr) = cx;
1227 i < retval.
cidx (j+1); i++)
1228 atmp += std::abs (retval.
data (i));
1229 if (atmp > ainvnorm)
1233 rcond = 1. / ainvnorm / anorm;
1241 double& rcond,
bool,
bool calc_cond)
const
1245 (*current_liboctave_error_handler)
1246 (
"inverse of the null matrix not defined");
1249 int typ = mattype.
type (
false);
1253 typ = mattype.
type (*
this);
1256 ret = dinverse (mattype, info, rcond,
true, calc_cond);
1258 ret = tinverse (mattype, info, rcond,
true, calc_cond).
transpose ();
1262 ret =
transpose ().tinverse (newtype, info, rcond,
true, calc_cond);
1269 octave::math::sparse_chol<SparseComplexMatrix> fact (*
this, info,
false);
1270 rcond = fact.rcond ();
1276 tinverse (tmp_typ, info, rcond2,
1295 octave::math::sparse_lu<SparseComplexMatrix> fact (*
this,
1298 rcond = fact.rcond ();
1305 octave::numeric_limits<double>::Inf ());
1306 std::copy_n (
ridx (), nz, ret.
xridx ());
1313 tinverse (tmp_typ, info, rcond2,
1316 tinverse (tmp_typ, info, rcond2,
1318 ret = fact.Pc ().
transpose () * InvU * InvL * fact.Pr ();
1346#if defined (HAVE_UMFPACK)
1351 if (nr == 0 || nc == 0 || nr != nc)
1360 Matrix Control (UMFPACK_CONTROL, 1);
1361 double *control = Control.
rwdata ();
1364 double tmp = octave::sparse_params::get_key (
"spumoni");
1365 if (! octave::math::isnan (tmp))
1366 Control (UMFPACK_PRL) = tmp;
1368 tmp = octave::sparse_params::get_key (
"piv_tol");
1369 if (! octave::math::isnan (tmp))
1371 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
1372 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
1376 tmp = octave::sparse_params::get_key (
"autoamd");
1377 if (! octave::math::isnan (tmp))
1378 Control (UMFPACK_FIXQ) = tmp;
1381 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1390 octave::to_suitesparse_intptr (Ap),
1391 octave::to_suitesparse_intptr (Ai),
1392 reinterpret_cast<const double *
> (Ax),
1393 nullptr, 1, control);
1396 Matrix Info (1, UMFPACK_INFO);
1397 double *info = Info.
rwdata ();
1399 octave::to_suitesparse_intptr (Ap),
1400 octave::to_suitesparse_intptr (Ai),
1401 reinterpret_cast<const double *
> (Ax),
1402 nullptr,
nullptr, &Symbolic, control, info);
1411 (*current_liboctave_error_handler)
1412 (
"SparseComplexMatrix::determinant symbolic factorization failed");
1420 =
UMFPACK_ZNAME (numeric) (octave::to_suitesparse_intptr (Ap),
1421 octave::to_suitesparse_intptr (Ai),
1422 reinterpret_cast<const double *
> (Ax),
1423 nullptr, Symbolic, &Numeric, control, info);
1426 rcond = Info (UMFPACK_RCOND);
1435 (*current_liboctave_error_handler)
1436 (
"SparseComplexMatrix::determinant numeric factorization failed");
1444 status =
UMFPACK_ZNAME (get_determinant) (c10,
nullptr, &e10,
1452 (*current_liboctave_error_handler)
1453 (
"SparseComplexMatrix::determinant error calculating determinant");
1465 octave_unused_parameter (err);
1466 octave_unused_parameter (rcond);
1468 (*current_liboctave_error_handler)
1469 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
1479 solve_singularity_handler,
bool calc_cond)
const
1488 if (nr != b.
rows ())
1490 (
"matrix dimension mismatch solution of linear equations");
1492 if (nr == 0 || nc == 0 || b.
cols () == 0)
1497 int typ = mattype.
type ();
1501 (*current_liboctave_error_handler) (
"incorrect matrix type");
1507 retval(i, j) = b(i, j) /
data (i);
1512 retval(k, j) = b(
ridx (i), j) /
data (i);
1517 double dmin = octave::numeric_limits<double>::Inf ();
1520 double tmp = std::abs (
data (i));
1526 rcond = dmin / dmax;
1538 solve_singularity_handler,
1539 bool calc_cond)
const
1548 if (nr != b.
rows ())
1550 (
"matrix dimension mismatch solution of linear equations");
1552 if (nr == 0 || nc == 0 || b.
cols () == 0)
1557 int typ = mattype.
type ();
1561 (*current_liboctave_error_handler) (
"incorrect matrix type");
1567 retval.
xcidx (0) = 0;
1574 if (b.
ridx (i) >= nm)
1579 retval.
xcidx (j+1) = ii;
1589 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1597 retval.
xridx (ii) = l;
1601 retval.
xcidx (j+1) = ii;
1607 double dmin = octave::numeric_limits<double>::Inf ();
1610 double tmp = std::abs (
data (i));
1616 rcond = dmin / dmax;
1628 solve_singularity_handler,
1629 bool calc_cond)
const
1638 if (nr != b.
rows ())
1640 (
"matrix dimension mismatch solution of linear equations");
1642 if (nr == 0 || nc == 0 || b.
cols () == 0)
1647 int typ = mattype.
type ();
1651 (*current_liboctave_error_handler) (
"incorrect matrix type");
1657 retval(i, j) = b(i, j) /
data (i);
1662 retval(k, j) = b(
ridx (i), j) /
data (i);
1667 double dmin = octave::numeric_limits<double>::Inf ();
1670 double tmp = std::abs (
data (i));
1676 rcond = dmin / dmax;
1688 solve_singularity_handler,
1689 bool calc_cond)
const
1698 if (nr != b.
rows ())
1700 (
"matrix dimension mismatch solution of linear equations");
1702 if (nr == 0 || nc == 0 || b.
cols () == 0)
1707 int typ = mattype.
type ();
1711 (*current_liboctave_error_handler) (
"incorrect matrix type");
1717 retval.
xcidx (0) = 0;
1724 if (b.
ridx (i) >= nm)
1729 retval.
xcidx (j+1) = ii;
1739 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1747 retval.
xridx (ii) = l;
1751 retval.
xcidx (j+1) = ii;
1757 double dmin = octave::numeric_limits<double>::Inf ();
1760 double tmp = std::abs (
data (i));
1766 rcond = dmin / dmax;
1778 solve_singularity_handler sing_handler,
1779 bool calc_cond)
const
1788 if (nr != b.
rows ())
1789 (*current_liboctave_error_handler)
1790 (
"matrix dimension mismatch solution of linear equations");
1792 if (nr == 0 || nc == 0 || b.
cols () == 0)
1797 int typ = mattype.
type ();
1801 (*current_liboctave_error_handler) (
"incorrect matrix type");
1804 double ainvnorm = 0.;
1815 atmp += std::abs (
data (i));
1823 retval.
resize (nc, b_nc);
1844 goto triangular_error;
1850 i <
cidx (kidx+1)-1; i++)
1853 work[iidx] = work[iidx] - tmp *
data (i);
1859 retval(perm[i], j) = work[i];
1881 i <
cidx (iidx+1)-1; i++)
1884 work[idx2] = work[idx2] - tmp *
data (i);
1891 atmp += std::abs (work[i]);
1894 if (atmp > ainvnorm)
1897 rcond = 1. / ainvnorm / anorm;
1903 retval.
resize (nc, b_nc);
1920 goto triangular_error;
1928 work[iidx] = work[iidx] - tmp *
data (i);
1934 retval.
xelem (i, j) = work[i];
1954 i <
cidx (k+1)-1; i++)
1957 work[iidx] = work[iidx] - tmp *
data (i);
1964 atmp += std::abs (work[i]);
1967 if (atmp > ainvnorm)
1970 rcond = 1. / ainvnorm / anorm;
1979 sing_handler (rcond);
1983 octave::warn_singular_matrix (rcond);
1986 if (is_singular (rcond) || octave::math::isnan (rcond))
1992 sing_handler (rcond);
1996 octave::warn_singular_matrix (rcond);
2006 solve_singularity_handler sing_handler,
2007 bool calc_cond)
const
2016 if (nr != b.
rows ())
2017 (*current_liboctave_error_handler)
2018 (
"matrix dimension mismatch solution of linear equations");
2020 if (nr == 0 || nc == 0 || b.
cols () == 0)
2025 int typ = mattype.
type ();
2029 (*current_liboctave_error_handler) (
"incorrect matrix type");
2032 double ainvnorm = 0.;
2042 atmp += std::abs (
data (i));
2051 retval.
xcidx (0) = 0;
2081 goto triangular_error;
2087 i <
cidx (kidx+1)-1; i++)
2090 work[iidx] = work[iidx] - tmp *
data (i);
2102 if (ii + new_nnz > x_nz)
2111 if (work[rperm[i]] != 0.)
2113 retval.
xridx (ii) = i;
2114 retval.
xdata (ii++) = work[rperm[i]];
2116 retval.
xcidx (j+1) = ii;
2140 i <
cidx (iidx+1)-1; i++)
2143 work[idx2] = work[idx2] - tmp *
data (i);
2150 atmp += std::abs (work[i]);
2153 if (atmp > ainvnorm)
2156 rcond = 1. / ainvnorm / anorm;
2178 goto triangular_error;
2186 work[iidx] = work[iidx] - tmp *
data (i);
2198 if (ii + new_nnz > x_nz)
2209 retval.
xridx (ii) = i;
2210 retval.
xdata (ii++) = work[i];
2212 retval.
xcidx (j+1) = ii;
2234 i <
cidx (k+1)-1; i++)
2237 work[iidx] = work[iidx] - tmp *
data (i);
2244 atmp += std::abs (work[i]);
2247 if (atmp > ainvnorm)
2250 rcond = 1. / ainvnorm / anorm;
2259 sing_handler (rcond);
2263 octave::warn_singular_matrix (rcond);
2266 if (is_singular (rcond) || octave::math::isnan (rcond))
2272 sing_handler (rcond);
2276 octave::warn_singular_matrix (rcond);
2285 solve_singularity_handler sing_handler,
2286 bool calc_cond)
const
2295 if (nr != b.
rows ())
2296 (*current_liboctave_error_handler)
2297 (
"matrix dimension mismatch solution of linear equations");
2299 if (nr == 0 || nc == 0 || b.
cols () == 0)
2304 int typ = mattype.
type ();
2308 (*current_liboctave_error_handler) (
"incorrect matrix type");
2311 double ainvnorm = 0.;
2322 atmp += std::abs (
data (i));
2330 retval.
resize (nc, b_nc);
2351 goto triangular_error;
2357 i <
cidx (kidx+1)-1; i++)
2360 work[iidx] = work[iidx] - tmp *
data (i);
2366 retval(perm[i], j) = work[i];
2388 i <
cidx (iidx+1)-1; i++)
2391 work[idx2] = work[idx2] - tmp *
data (i);
2398 atmp += std::abs (work[i]);
2401 if (atmp > ainvnorm)
2404 rcond = 1. / ainvnorm / anorm;
2410 retval.
resize (nc, b_nc);
2427 goto triangular_error;
2435 work[iidx] = work[iidx] - tmp *
data (i);
2441 retval.
xelem (i, j) = work[i];
2461 i <
cidx (k+1)-1; i++)
2464 work[iidx] = work[iidx] - tmp *
data (i);
2471 atmp += std::abs (work[i]);
2474 if (atmp > ainvnorm)
2477 rcond = 1. / ainvnorm / anorm;
2486 sing_handler (rcond);
2490 octave::warn_singular_matrix (rcond);
2493 if (is_singular (rcond) || octave::math::isnan (rcond))
2499 sing_handler (rcond);
2503 octave::warn_singular_matrix (rcond);
2513 solve_singularity_handler sing_handler,
2514 bool calc_cond)
const
2523 if (nr != b.
rows ())
2524 (*current_liboctave_error_handler)
2525 (
"matrix dimension mismatch solution of linear equations");
2527 if (nr == 0 || nc == 0 || b.
cols () == 0)
2532 int typ = mattype.
type ();
2536 (*current_liboctave_error_handler) (
"incorrect matrix type");
2539 double ainvnorm = 0.;
2549 atmp += std::abs (
data (i));
2558 retval.
xcidx (0) = 0;
2588 goto triangular_error;
2594 i <
cidx (kidx+1)-1; i++)
2597 work[iidx] = work[iidx] - tmp *
data (i);
2609 if (ii + new_nnz > x_nz)
2618 if (work[rperm[i]] != 0.)
2620 retval.
xridx (ii) = i;
2621 retval.
xdata (ii++) = work[rperm[i]];
2623 retval.
xcidx (j+1) = ii;
2647 i <
cidx (iidx+1)-1; i++)
2650 work[idx2] = work[idx2] - tmp *
data (i);
2657 atmp += std::abs (work[i]);
2660 if (atmp > ainvnorm)
2663 rcond = 1. / ainvnorm / anorm;
2685 goto triangular_error;
2693 work[iidx] = work[iidx] - tmp *
data (i);
2705 if (ii + new_nnz > x_nz)
2716 retval.
xridx (ii) = i;
2717 retval.
xdata (ii++) = work[i];
2719 retval.
xcidx (j+1) = ii;
2741 i <
cidx (k+1)-1; i++)
2744 work[iidx] = work[iidx] - tmp *
data (i);
2751 atmp += std::abs (work[i]);
2754 if (atmp > ainvnorm)
2757 rcond = 1. / ainvnorm / anorm;
2766 sing_handler (rcond);
2770 octave::warn_singular_matrix (rcond);
2773 if (is_singular (rcond) || octave::math::isnan (rcond))
2779 sing_handler (rcond);
2783 octave::warn_singular_matrix (rcond);
2793 solve_singularity_handler sing_handler,
2794 bool calc_cond)
const
2803 if (nr != b.
rows ())
2804 (*current_liboctave_error_handler)
2805 (
"matrix dimension mismatch solution of linear equations");
2807 if (nr == 0 || nc == 0 || b.
cols () == 0)
2812 int typ = mattype.
type ();
2816 (*current_liboctave_error_handler) (
"incorrect matrix type");
2819 double ainvnorm = 0.;
2830 atmp += std::abs (
data (i));
2838 retval.
resize (nc, b_nc);
2844 inv_perm[perm[i]] = i;
2851 work[inv_perm[i]] = b(i, j);
2862 if (inv_perm[
ridx (i)] < minr)
2864 minr = inv_perm[
ridx (i)];
2868 if (minr != k ||
data (mini) == 0.)
2871 goto triangular_error;
2882 work[iidx] = work[iidx] - tmp *
data (i);
2888 retval(i, j) = work[i];
2909 i <
cidx (k+1); i++)
2910 if (inv_perm[
ridx (i)] < minr)
2912 minr = inv_perm[
ridx (i)];
2919 i <
cidx (k+1); i++)
2925 work[iidx] = work[iidx] - tmp *
data (i);
2933 atmp += std::abs (work[i]);
2936 if (atmp > ainvnorm)
2939 rcond = 1. / ainvnorm / anorm;
2945 retval.
resize (nc, b_nc, 0.);
2960 goto triangular_error;
2968 work[iidx] = work[iidx] - tmp *
data (i);
2973 retval.
xelem (i, j) = work[i];
2994 i <
cidx (k+1); i++)
2997 work[iidx] = work[iidx] - tmp *
data (i);
3004 atmp += std::abs (work[i]);
3007 if (atmp > ainvnorm)
3010 rcond = 1. / ainvnorm / anorm;
3018 sing_handler (rcond);
3022 octave::warn_singular_matrix (rcond);
3025 if (is_singular (rcond) || octave::math::isnan (rcond))
3031 sing_handler (rcond);
3035 octave::warn_singular_matrix (rcond);
3045 solve_singularity_handler sing_handler,
3046 bool calc_cond)
const
3056 if (nr != b.
rows ())
3057 (*current_liboctave_error_handler)
3058 (
"matrix dimension mismatch solution of linear equations");
3060 if (nr == 0 || nc == 0 || b.
cols () == 0)
3065 int typ = mattype.
type ();
3069 (*current_liboctave_error_handler) (
"incorrect matrix type");
3072 double ainvnorm = 0.;
3082 atmp += std::abs (
data (i));
3091 retval.
xcidx (0) = 0;
3102 inv_perm[perm[i]] = i;
3109 work[inv_perm[b.
ridx (i)]] = b.
data (i);
3120 if (inv_perm[
ridx (i)] < minr)
3122 minr = inv_perm[
ridx (i)];
3126 if (minr != k ||
data (mini) == 0.)
3129 goto triangular_error;
3140 work[iidx] = work[iidx] - tmp *
data (i);
3152 if (ii + new_nnz > x_nz)
3163 retval.
xridx (ii) = i;
3164 retval.
xdata (ii++) = work[i];
3166 retval.
xcidx (j+1) = ii;
3189 i <
cidx (k+1); i++)
3190 if (inv_perm[
ridx (i)] < minr)
3192 minr = inv_perm[
ridx (i)];
3199 i <
cidx (k+1); i++)
3205 work[iidx] = work[iidx] - tmp *
data (i);
3213 atmp += std::abs (work[i]);
3216 if (atmp > ainvnorm)
3219 rcond = 1. / ainvnorm / anorm;
3240 goto triangular_error;
3248 work[iidx] = work[iidx] - tmp *
data (i);
3260 if (ii + new_nnz > x_nz)
3271 retval.
xridx (ii) = i;
3272 retval.
xdata (ii++) = work[i];
3274 retval.
xcidx (j+1) = ii;
3297 i <
cidx (k+1); i++)
3300 work[iidx] = work[iidx] - tmp *
data (i);
3307 atmp += std::abs (work[i]);
3310 if (atmp > ainvnorm)
3313 rcond = 1. / ainvnorm / anorm;
3322 sing_handler (rcond);
3326 octave::warn_singular_matrix (rcond);
3329 if (is_singular (rcond) || octave::math::isnan (rcond))
3335 sing_handler (rcond);
3339 octave::warn_singular_matrix (rcond);
3349 solve_singularity_handler sing_handler,
3350 bool calc_cond)
const
3359 if (nr != b.
rows ())
3360 (*current_liboctave_error_handler)
3361 (
"matrix dimension mismatch solution of linear equations");
3363 if (nr == 0 || nc == 0 || b.
cols () == 0)
3368 int typ = mattype.
type ();
3372 (*current_liboctave_error_handler) (
"incorrect matrix type");
3375 double ainvnorm = 0.;
3386 atmp += std::abs (
data (i));
3394 retval.
resize (nc, b_nc);
3400 inv_perm[perm[i]] = i;
3407 work[inv_perm[i]] = b(i, j);
3417 if (inv_perm[
ridx (i)] < minr)
3419 minr = inv_perm[
ridx (i)];
3423 if (minr != k ||
data (mini) == 0.)
3426 goto triangular_error;
3437 work[iidx] = work[iidx] - tmp *
data (i);
3443 retval(i, j) = work[i];
3464 i <
cidx (k+1); i++)
3465 if (inv_perm[
ridx (i)] < minr)
3467 minr = inv_perm[
ridx (i)];
3474 i <
cidx (k+1); i++)
3480 work[iidx] = work[iidx] - tmp *
data (i);
3488 atmp += std::abs (work[i]);
3491 if (atmp > ainvnorm)
3494 rcond = 1. / ainvnorm / anorm;
3500 retval.
resize (nc, b_nc, 0.);
3516 goto triangular_error;
3524 work[iidx] = work[iidx] - tmp *
data (i);
3530 retval.
xelem (i, j) = work[i];
3551 i <
cidx (k+1); i++)
3554 work[iidx] = work[iidx] - tmp *
data (i);
3561 atmp += std::abs (work[i]);
3564 if (atmp > ainvnorm)
3567 rcond = 1. / ainvnorm / anorm;
3576 sing_handler (rcond);
3580 octave::warn_singular_matrix (rcond);
3583 if (is_singular (rcond) || octave::math::isnan (rcond))
3589 sing_handler (rcond);
3593 octave::warn_singular_matrix (rcond);
3603 solve_singularity_handler sing_handler,
3604 bool calc_cond)
const
3613 if (nr != b.
rows ())
3614 (*current_liboctave_error_handler)
3615 (
"matrix dimension mismatch solution of linear equations");
3617 if (nr == 0 || nc == 0 || b.
cols () == 0)
3622 int typ = mattype.
type ();
3626 (*current_liboctave_error_handler) (
"incorrect matrix type");
3629 double ainvnorm = 0.;
3639 atmp += std::abs (
data (i));
3648 retval.
xcidx (0) = 0;
3659 inv_perm[perm[i]] = i;
3666 work[inv_perm[b.
ridx (i)]] = b.
data (i);
3676 if (inv_perm[
ridx (i)] < minr)
3678 minr = inv_perm[
ridx (i)];
3682 if (minr != k ||
data (mini) == 0.)
3685 goto triangular_error;
3696 work[iidx] = work[iidx] - tmp *
data (i);
3708 if (ii + new_nnz > x_nz)
3719 retval.
xridx (ii) = i;
3720 retval.
xdata (ii++) = work[i];
3722 retval.
xcidx (j+1) = ii;
3745 i <
cidx (k+1); i++)
3746 if (inv_perm[
ridx (i)] < minr)
3748 minr = inv_perm[
ridx (i)];
3755 i <
cidx (k+1); i++)
3761 work[iidx] = work[iidx] - tmp *
data (i);
3769 atmp += std::abs (work[i]);
3772 if (atmp > ainvnorm)
3775 rcond = 1. / ainvnorm / anorm;
3796 goto triangular_error;
3804 work[iidx] = work[iidx] - tmp *
data (i);
3816 if (ii + new_nnz > x_nz)
3827 retval.
xridx (ii) = i;
3828 retval.
xdata (ii++) = work[i];
3830 retval.
xcidx (j+1) = ii;
3853 i <
cidx (k+1); i++)
3856 work[iidx] = work[iidx] - tmp *
data (i);
3863 atmp += std::abs (work[i]);
3866 if (atmp > ainvnorm)
3869 rcond = 1. / ainvnorm / anorm;
3878 sing_handler (rcond);
3882 octave::warn_singular_matrix (rcond);
3885 if (is_singular (rcond) || octave::math::isnan (rcond))
3891 sing_handler (rcond);
3895 octave::warn_singular_matrix (rcond);
3905 solve_singularity_handler sing_handler,
3906 bool calc_cond)
const
3914 if (nr != nc || nr != b.
rows ())
3915 (*current_liboctave_error_handler)
3916 (
"matrix dimension mismatch solution of linear equations");
3918 if (nr == 0 || b.
cols () == 0)
3921 (*current_liboctave_error_handler)
3922 (
"calculation of condition number not implemented");
3926 int typ = mattype.
type ();
3940 D[j] = std::real (
data (ii++));
3944 D[nc-1] = std::real (
data (ii));
3959 D[j] = std::real (
data (i));
3960 else if (
ridx (i) == j + 1)
3971 F77_INT tmp_nr = octave::to_f77_int (nr);
4004 DL[j] =
data (ii++);
4005 DU[j] =
data (ii++);
4007 D[nc-1] =
data (ii);
4024 else if (
ridx (i) == j + 1)
4026 else if (
ridx (i) == j - 1)
4037 F77_INT tmp_nr = octave::to_f77_int (nr);
4054 sing_handler (rcond);
4058 octave::warn_singular_matrix ();
4065 (*current_liboctave_error_handler) (
"incorrect matrix type");
4074 solve_singularity_handler sing_handler,
4075 bool calc_cond)
const
4083 if (nr != nc || nr != b.
rows ())
4084 (*current_liboctave_error_handler)
4085 (
"matrix dimension mismatch solution of linear equations");
4087 if (nr == 0 || b.
cols () == 0)
4090 (*current_liboctave_error_handler)
4091 (
"calculation of condition number not implemented");
4095 int typ = mattype.
type ();
4107 F77_INT *pipvt = ipvt.rwdata ();
4116 DL[j] =
data (ii++);
4117 DU[j] =
data (ii++);
4119 D[nc-1] =
data (ii);
4136 else if (
ridx (i) == j + 1)
4138 else if (
ridx (i) == j - 1)
4143 F77_INT tmp_nr = octave::to_f77_int (nr);
4162 sing_handler (rcond);
4166 octave::warn_singular_matrix ();
4175 retval.
xcidx (0) = 0;
4189 (F77_CONST_CHAR_ARG2 (&job, 1),
4195 F77_CHAR_ARG_LEN (1)));
4206 if (ii + new_nnz > x_nz)
4217 retval.
xridx (ii) = i;
4218 retval.
xdata (ii++) = work[i];
4220 retval.
xcidx (j+1) = ii;
4227 (*current_liboctave_error_handler) (
"incorrect matrix type");
4236 solve_singularity_handler sing_handler,
4237 bool calc_cond)
const
4245 if (nr != nc || nr != b.
rows ())
4246 (*current_liboctave_error_handler)
4247 (
"matrix dimension mismatch solution of linear equations");
4249 if (nr == 0 || b.
cols () == 0)
4252 (*current_liboctave_error_handler)
4253 (
"calculation of condition number not implemented");
4257 int typ = mattype.
type ();
4271 D[j] = std::real (
data (ii++));
4275 D[nc-1] = std::real (
data (ii));
4290 D[j] = std::real (
data (i));
4291 else if (
ridx (i) == j + 1)
4304 F77_INT tmp_nr = octave::to_f77_int (nr);
4335 DL[j] =
data (ii++);
4336 DU[j] =
data (ii++);
4338 D[nc-1] =
data (ii);
4355 else if (
ridx (i) == j + 1)
4357 else if (
ridx (i) == j - 1)
4370 F77_INT tmp_nr = octave::to_f77_int (nr);
4387 sing_handler (rcond);
4391 octave::warn_singular_matrix ();
4395 (*current_liboctave_error_handler) (
"incorrect matrix type");
4402SparseComplexMatrix::trisolve (
MatrixType& mattype,
4405 solve_singularity_handler sing_handler,
4406 bool calc_cond)
const
4414 if (nr != nc || nr != b.
rows ())
4415 (*current_liboctave_error_handler)
4416 (
"matrix dimension mismatch solution of linear equations");
4418 if (nr == 0 || b.
cols () == 0)
4421 (*current_liboctave_error_handler)
4422 (
"calculation of condition number not implemented");
4426 int typ = mattype.
type ();
4438 F77_INT *pipvt = ipvt.rwdata ();
4447 DL[j] =
data (ii++);
4448 DU[j] =
data (ii++);
4450 D[nc-1] =
data (ii);
4467 else if (
ridx (i) == j + 1)
4469 else if (
ridx (i) == j - 1)
4474 F77_INT tmp_nr = octave::to_f77_int (nr);
4493 sing_handler (rcond);
4497 octave::warn_singular_matrix ();
4513 retval.
xcidx (0) = 0;
4517 for (
F77_INT i = 0; i < b_nr; i++)
4521 (F77_CONST_CHAR_ARG2 (&job, 1),
4527 F77_CHAR_ARG_LEN (1)));
4535 (*current_liboctave_error_handler)
4536 (
"SparseComplexMatrix::solve solve failed");
4549 if (ii + new_nnz > x_nz)
4560 retval.
xridx (ii) = i;
4561 retval.
xdata (ii++) = Bx[i];
4564 retval.
xcidx (j+1) = ii;
4571 (*current_liboctave_error_handler) (
"incorrect matrix type");
4580 solve_singularity_handler sing_handler,
4581 bool calc_cond)
const
4589 if (nr != nc || nr != b.
rows ())
4590 (*current_liboctave_error_handler)
4591 (
"matrix dimension mismatch solution of linear equations");
4593 if (nr == 0 || b.
cols () == 0)
4598 int typ = mattype.
type ();
4606 Complex *tmp_data = m_band.rwdata ();
4612 for (
F77_INT j = 0; j < ldm; j++)
4614 tmp_data[ii++] = 0.;
4622 m_band(ri - j, j) =
data (i);
4628 anorm = m_band.abs ().sum ().row (0).max ();
4630 F77_INT tmp_nr = octave::to_f77_int (nr);
4635 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4637 F77_CHAR_ARG_LEN (1)));
4657 double *piz = iz.rwdata ();
4660 (F77_CONST_CHAR_ARG2 (&job, 1),
4663 F77_CHAR_ARG_LEN (1)));
4670 if (is_singular (rcond) || octave::math::isnan (rcond))
4676 sing_handler (rcond);
4680 octave::warn_singular_matrix (rcond);
4695 (F77_CONST_CHAR_ARG2 (&job, 1),
4698 F77_CHAR_ARG_LEN (1)));
4705 (*current_liboctave_error_handler)
4706 (
"SparseMatrix::solve solve failed");
4718 F77_INT ldm = n_upper + 2 * n_lower + 1;
4721 Complex *tmp_data = m_band.rwdata ();
4727 for (
F77_INT j = 0; j < ldm; j++)
4729 tmp_data[ii++] = 0.;
4734 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
4744 atmp += std::abs (
data (i));
4751 F77_INT *pipvt = ipvt.rwdata ();
4753 F77_INT tmp_nr = octave::to_f77_int (nr);
4757 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4759 ldm, pipvt, tmp_err));
4772 sing_handler (rcond);
4776 octave::warn_singular_matrix ();
4786 double *piz = iz.rwdata ();
4788 F77_INT tmp_nc = octave::to_f77_int (nc);
4791 (F77_CONST_CHAR_ARG2 (&job, 1),
4794 F77_CHAR_ARG_LEN (1)));
4801 if (is_singular (rcond) || octave::math::isnan (rcond))
4807 sing_handler (rcond);
4811 octave::warn_singular_matrix (rcond);
4827 (F77_CONST_CHAR_ARG2 (&job, 1),
4830 F77_CHAR_ARG_LEN (1)));
4837 (*current_liboctave_error_handler) (
"incorrect matrix type");
4846 solve_singularity_handler sing_handler,
4847 bool calc_cond)
const
4855 if (nr != nc || nr != b.
rows ())
4856 (*current_liboctave_error_handler)
4857 (
"matrix dimension mismatch solution of linear equations");
4859 if (nr == 0 || b.
cols () == 0)
4864 int typ = mattype.
type ();
4873 Complex *tmp_data = m_band.rwdata ();
4879 for (
F77_INT j = 0; j < ldm; j++)
4881 tmp_data[ii++] = 0.;
4889 m_band(ri - j, j) =
data (i);
4895 anorm = m_band.abs ().sum ().row (0).max ();
4897 F77_INT tmp_nr = octave::to_f77_int (nr);
4902 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4904 F77_CHAR_ARG_LEN (1)));
4922 double *piz = iz.rwdata ();
4925 (F77_CONST_CHAR_ARG2 (&job, 1),
4928 F77_CHAR_ARG_LEN (1)));
4935 if (is_singular (rcond) || octave::math::isnan (rcond))
4941 sing_handler (rcond);
4945 octave::warn_singular_matrix (rcond);
4963 retval.
xcidx (0) = 0;
4966 for (
F77_INT i = 0; i < b_nr; i++)
4967 Bx[i] = b.
elem (i, j);
4970 (F77_CONST_CHAR_ARG2 (&job, 1),
4973 F77_CHAR_ARG_LEN (1)));
4980 (*current_liboctave_error_handler)
4981 (
"SparseComplexMatrix::solve solve failed");
4995 sz = (
static_cast<double> (b_nc) - j) / b_nc
4997 sz = x_nz + (sz > 100 ? sz : 100);
5001 retval.
xdata (ii) = tmp;
5002 retval.
xridx (ii++) = i;
5005 retval.
xcidx (j+1) = ii;
5018 F77_INT ldm = n_upper + 2 * n_lower + 1;
5021 Complex *tmp_data = m_band.rwdata ();
5027 for (
F77_INT j = 0; j < ldm; j++)
5029 tmp_data[ii++] = 0.;
5034 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
5044 atmp += std::abs (
data (i));
5051 F77_INT *pipvt = ipvt.rwdata ();
5053 F77_INT tmp_nr = octave::to_f77_int (nr);
5057 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5059 ldm, pipvt, tmp_err));
5070 sing_handler (rcond);
5074 octave::warn_singular_matrix ();
5084 double *piz = iz.rwdata ();
5086 F77_INT tmp_nc = octave::to_f77_int (nc);
5089 (F77_CONST_CHAR_ARG2 (&job, 1),
5092 F77_CHAR_ARG_LEN (1)));
5099 if (is_singular (rcond) || octave::math::isnan (rcond))
5105 sing_handler (rcond);
5109 octave::warn_singular_matrix (rcond);
5122 retval.
xcidx (0) = 0;
5132 i < b.
cidx (j+1); i++)
5136 (F77_CONST_CHAR_ARG2 (&job, 1),
5139 F77_CHAR_ARG_LEN (1)));
5150 if (ii + new_nnz > x_nz)
5161 retval.
xridx (ii) = i;
5162 retval.
xdata (ii++) = work[i];
5164 retval.
xcidx (j+1) = ii;
5172 (*current_liboctave_error_handler) (
"incorrect matrix type");
5181 solve_singularity_handler sing_handler,
5182 bool calc_cond)
const
5190 if (nr != nc || nr != b.
rows ())
5191 (*current_liboctave_error_handler)
5192 (
"matrix dimension mismatch solution of linear equations");
5194 if (nr == 0 || b.
cols () == 0)
5199 int typ = mattype.
type ();
5208 Complex *tmp_data = m_band.rwdata ();
5214 for (
F77_INT j = 0; j < ldm; j++)
5216 tmp_data[ii++] = 0.;
5224 m_band(ri - j, j) =
data (i);
5230 anorm = m_band.abs ().sum ().row (0).max ();
5232 F77_INT tmp_nr = octave::to_f77_int (nr);
5237 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5239 F77_CHAR_ARG_LEN (1)));
5259 double *piz = iz.rwdata ();
5262 (F77_CONST_CHAR_ARG2 (&job, 1),
5265 F77_CHAR_ARG_LEN (1)));
5272 if (is_singular (rcond) || octave::math::isnan (rcond))
5278 sing_handler (rcond);
5282 octave::warn_singular_matrix (rcond);
5296 (F77_CONST_CHAR_ARG2 (&job, 1),
5299 F77_CHAR_ARG_LEN (1)));
5306 (*current_liboctave_error_handler)
5307 (
"SparseComplexMatrix::solve solve failed");
5319 F77_INT ldm = n_upper + 2 * n_lower + 1;
5322 Complex *tmp_data = m_band.rwdata ();
5328 for (
F77_INT j = 0; j < ldm; j++)
5330 tmp_data[ii++] = 0.;
5335 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
5345 atmp += std::abs (
data (i));
5352 F77_INT *pipvt = ipvt.rwdata ();
5354 F77_INT tmp_nr = octave::to_f77_int (nr);
5358 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5360 ldm, pipvt, tmp_err));
5371 sing_handler (rcond);
5375 octave::warn_singular_matrix ();
5385 double *piz = iz.rwdata ();
5387 F77_INT tmp_nc = octave::to_f77_int (nc);
5390 (F77_CONST_CHAR_ARG2 (&job, 1),
5393 F77_CHAR_ARG_LEN (1)));
5400 if (is_singular (rcond) || octave::math::isnan (rcond))
5406 sing_handler (rcond);
5410 octave::warn_singular_matrix (rcond);
5425 (F77_CONST_CHAR_ARG2 (&job, 1),
5428 F77_CHAR_ARG_LEN (1)));
5435 (*current_liboctave_error_handler) (
"incorrect matrix type");
5444 solve_singularity_handler sing_handler,
5445 bool calc_cond)
const
5453 if (nr != nc || nr != b.
rows ())
5454 (*current_liboctave_error_handler)
5455 (
"matrix dimension mismatch solution of linear equations");
5457 if (nr == 0 || b.
cols () == 0)
5462 int typ = mattype.
type ();
5471 Complex *tmp_data = m_band.rwdata ();
5477 for (
F77_INT j = 0; j < ldm; j++)
5479 tmp_data[ii++] = 0.;
5487 m_band(ri - j, j) =
data (i);
5493 anorm = m_band.abs ().sum ().row (0).max ();
5495 F77_INT tmp_nr = octave::to_f77_int (nr);
5500 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5502 F77_CHAR_ARG_LEN (1)));
5523 double *piz = iz.rwdata ();
5526 (F77_CONST_CHAR_ARG2 (&job, 1),
5529 F77_CHAR_ARG_LEN (1)));
5536 if (is_singular (rcond) || octave::math::isnan (rcond))
5542 sing_handler (rcond);
5546 octave::warn_singular_matrix (rcond);
5564 retval.
xcidx (0) = 0;
5568 for (
F77_INT i = 0; i < b_nr; i++)
5572 (F77_CONST_CHAR_ARG2 (&job, 1),
5575 F77_CHAR_ARG_LEN (1)));
5582 (*current_liboctave_error_handler)
5583 (
"SparseMatrix::solve solve failed");
5595 if (ii + new_nnz > x_nz)
5606 retval.
xridx (ii) = i;
5607 retval.
xdata (ii++) = Bx[i];
5610 retval.
xcidx (j+1) = ii;
5623 F77_INT ldm = n_upper + 2 * n_lower + 1;
5626 Complex *tmp_data = m_band.rwdata ();
5632 for (
F77_INT j = 0; j < ldm; j++)
5634 tmp_data[ii++] = 0.;
5639 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
5649 atmp += std::abs (
data (i));
5656 F77_INT *pipvt = ipvt.rwdata ();
5658 F77_INT tmp_nr = octave::to_f77_int (nr);
5662 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5664 ldm, pipvt, tmp_err));
5675 sing_handler (rcond);
5679 octave::warn_singular_matrix ();
5689 double *piz = iz.rwdata ();
5691 F77_INT tmp_nc = octave::to_f77_int (nc);
5694 (F77_CONST_CHAR_ARG2 (&job, 1),
5697 F77_CHAR_ARG_LEN (1)));
5704 if (is_singular (rcond) || octave::math::isnan (rcond))
5710 sing_handler (rcond);
5714 octave::warn_singular_matrix (rcond);
5727 retval.
xcidx (0) = 0;
5738 i < b.
cidx (j+1); i++)
5742 (F77_CONST_CHAR_ARG2 (&job, 1),
5745 F77_CHAR_ARG_LEN (1)));
5756 if (ii + new_nnz > x_nz)
5767 retval.
xridx (ii) = i;
5768 retval.
xdata (ii++) = Bx[i];
5770 retval.
xcidx (j+1) = ii;
5778 (*current_liboctave_error_handler) (
"incorrect matrix type");
5787 solve_singularity_handler sing_handler,
5788 bool calc_cond)
const
5791 void *Numeric =
nullptr;
5794#if defined (HAVE_UMFPACK)
5797 Control =
Matrix (UMFPACK_CONTROL, 1);
5798 double *control = Control.
rwdata ();
5801 double tmp = octave::sparse_params::get_key (
"spumoni");
5802 if (! octave::math::isnan (tmp))
5803 Control (UMFPACK_PRL) = tmp;
5804 tmp = octave::sparse_params::get_key (
"piv_tol");
5805 if (! octave::math::isnan (tmp))
5807 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
5808 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
5812 tmp = octave::sparse_params::get_key (
"autoamd");
5813 if (! octave::math::isnan (tmp))
5814 Control (UMFPACK_FIXQ) = tmp;
5825 octave::to_suitesparse_intptr (Ap),
5826 octave::to_suitesparse_intptr (Ai),
5827 reinterpret_cast<const double *
> (Ax),
5828 nullptr, 1, control);
5831 Info =
Matrix (1, UMFPACK_INFO);
5832 double *info = Info.
rwdata ();
5834 octave::to_suitesparse_intptr (Ap),
5835 octave::to_suitesparse_intptr (Ai),
5836 reinterpret_cast<const double *
> (Ax),
5837 nullptr,
nullptr, &Symbolic, control, info);
5847 (*current_liboctave_error_handler)
5848 (
"SparseComplexMatrix::solve symbolic factorization failed");
5855 status =
UMFPACK_ZNAME (numeric) (octave::to_suitesparse_intptr (Ap),
5856 octave::to_suitesparse_intptr (Ai),
5857 reinterpret_cast<const double *
> (Ax),
5858 nullptr, Symbolic, &Numeric, control, info);
5862 rcond = Info (UMFPACK_RCOND);
5866 if (status == UMFPACK_WARNING_singular_matrix
5867 || is_singular (rcond) || octave::math::isnan (rcond))
5874 sing_handler (rcond);
5876 octave::warn_singular_matrix (rcond);
5878 else if (status < 0)
5884 (*current_liboctave_error_handler)
5885 (
"SparseComplexMatrix::solve numeric factorization failed");
5900 octave_unused_parameter (rcond);
5901 octave_unused_parameter (Control);
5902 octave_unused_parameter (Info);
5903 octave_unused_parameter (sing_handler);
5904 octave_unused_parameter (calc_cond);
5906 (*current_liboctave_error_handler)
5907 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
5917 solve_singularity_handler sing_handler,
5918 bool calc_cond)
const
5926 if (nr != nc || nr != b.
rows ())
5927 (*current_liboctave_error_handler)
5928 (
"matrix dimension mismatch solution of linear equations");
5930 if (nr == 0 || b.
cols () == 0)
5935 int typ = mattype.
type ();
5940#if defined (HAVE_CHOLMOD)
5941 cholmod_common Common;
5942 cholmod_common *cm = &Common;
5946 cm->prefer_zomplex =
false;
5948 double spu = octave::sparse_params::get_key (
"spumoni");
5952 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
5956 cm->print =
static_cast<int> (spu) + 2;
5957 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5961 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5962 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5964 cm->final_ll =
true;
5966 cholmod_sparse Astore;
5967 cholmod_sparse *
A = &Astore;
5977#if defined (OCTAVE_ENABLE_64)
5978 A->itype = CHOLMOD_LONG;
5980 A->itype = CHOLMOD_INT;
5982 A->dtype = CHOLMOD_DOUBLE;
5984 A->xtype = CHOLMOD_COMPLEX;
5988 cholmod_dense Bstore;
5989 cholmod_dense *
B = &Bstore;
5990 B->nrow = b.
rows ();
5991 B->ncol = b.
cols ();
5993 B->nzmax =
B->nrow *
B->ncol;
5994 B->dtype = CHOLMOD_DOUBLE;
5995 B->xtype = CHOLMOD_REAL;
5997 B->x =
const_cast<double *
> (b.
data ());
6014 if (is_singular (rcond) || octave::math::isnan (rcond))
6020 sing_handler (rcond);
6024 octave::warn_singular_matrix (rcond);
6036 retval.
xelem (i, j) =
static_cast<Complex *
> (X->x)[jr + i];
6042 static char blank_name[] =
" ";
6046 (*current_liboctave_warning_with_id_handler)
6047 (
"Octave:missing-dependency",
6048 "support for CHOLMOD was unavailable or disabled "
6049 "when liboctave was built");
6058#if defined (HAVE_UMFPACK)
6060 void *Numeric = factorize (err, rcond, Control, Info,
6061 sing_handler, calc_cond);
6066 Control (UMFPACK_IRSTEP) = 1;
6070 double *control = Control.
rwdata ();
6071 double *info = Info.
rwdata ();
6075#if defined (UMFPACK_SEPARATE_SPLIT)
6076 const double *Bx = b.
data ();
6083 retval.
resize (b_nr, b_nc);
6088#if defined (UMFPACK_SEPARATE_SPLIT)
6090 octave::to_suitesparse_intptr (Ap),
6091 octave::to_suitesparse_intptr (Ai),
6092 reinterpret_cast<const double *
> (Ax),
6094 reinterpret_cast<double *
> (&Xx[iidx]),
6096 &Bx[iidx], Bz, Numeric,
6100 Bz[i] = b.
elem (i, j);
6103 octave::to_suitesparse_intptr (Ap),
6104 octave::to_suitesparse_intptr (Ai),
6105 reinterpret_cast<const double *
> (Ax),
6107 reinterpret_cast<double *
> (&Xx[iidx]),
6109 reinterpret_cast<const double *
> (Bz),
6119 (*current_liboctave_error_handler)
6120 (
"SparseComplexMatrix::solve solve failed");
6135 octave_unused_parameter (rcond);
6136 octave_unused_parameter (sing_handler);
6137 octave_unused_parameter (calc_cond);
6139 (*current_liboctave_error_handler)
6140 (
"support for UMFPACK was unavailable or disabled "
6141 "when liboctave was built");
6145 (*current_liboctave_error_handler) (
"incorrect matrix type");
6154 solve_singularity_handler sing_handler,
6155 bool calc_cond)
const
6163 if (nr != nc || nr != b.
rows ())
6164 (*current_liboctave_error_handler)
6165 (
"matrix dimension mismatch solution of linear equations");
6167 if (nr == 0 || b.
cols () == 0)
6172 int typ = mattype.
type ();
6177#if defined (HAVE_CHOLMOD)
6178 cholmod_common Common;
6179 cholmod_common *cm = &Common;
6183 cm->prefer_zomplex =
false;
6185 double spu = octave::sparse_params::get_key (
"spumoni");
6189 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6193 cm->print =
static_cast<int> (spu) + 2;
6194 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6198 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6199 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6201 cm->final_ll =
true;
6203 cholmod_sparse Astore;
6204 cholmod_sparse *
A = &Astore;
6214#if defined (OCTAVE_ENABLE_64)
6215 A->itype = CHOLMOD_LONG;
6217 A->itype = CHOLMOD_INT;
6219 A->dtype = CHOLMOD_DOUBLE;
6221 A->xtype = CHOLMOD_COMPLEX;
6225 cholmod_sparse Bstore;
6226 cholmod_sparse *
B = &Bstore;
6227 B->nrow = b.
rows ();
6228 B->ncol = b.
cols ();
6231 B->nzmax = b.
nnz ();
6235#if defined (OCTAVE_ENABLE_64)
6236 B->itype = CHOLMOD_LONG;
6238 B->itype = CHOLMOD_INT;
6240 B->dtype = CHOLMOD_DOUBLE;
6242 B->xtype = CHOLMOD_REAL;
6261 if (is_singular (rcond) || octave::math::isnan (rcond))
6267 sing_handler (rcond);
6271 octave::warn_singular_matrix (rcond);
6276 cholmod_sparse *X =
CHOLMOD_NAME(spsolve) (CHOLMOD_A, L,
B, cm);
6279 (octave::from_size_t (X->nrow),
6280 octave::from_size_t (X->ncol),
6284 j <= static_cast<octave_idx_type> (X->ncol); j++)
6297 static char blank_name[] =
" ";
6301 (*current_liboctave_warning_with_id_handler)
6302 (
"Octave:missing-dependency",
6303 "support for CHOLMOD was unavailable or disabled "
6304 "when liboctave was built");
6313#if defined (HAVE_UMFPACK)
6315 void *Numeric = factorize (err, rcond, Control, Info,
6316 sing_handler, calc_cond);
6321 Control (UMFPACK_IRSTEP) = 1;
6325 double *control = Control.
rwdata ();
6326 double *info = Info.
rwdata ();
6331#if defined (UMFPACK_SEPARATE_SPLIT)
6348 retval.
xcidx (0) = 0;
6352#if defined (UMFPACK_SEPARATE_SPLIT)
6354 Bx[i] = b.
elem (i, j);
6357 octave::to_suitesparse_intptr (Ap),
6358 octave::to_suitesparse_intptr (Ai),
6359 reinterpret_cast<const double *
> (Ax),
6361 reinterpret_cast<double *
> (Xx),
6363 Bx, Bz, Numeric, control,
6367 Bz[i] = b.
elem (i, j);
6370 octave::to_suitesparse_intptr (Ap),
6371 octave::to_suitesparse_intptr (Ai),
6372 reinterpret_cast<const double *
> (Ax),
6374 reinterpret_cast<double *
> (Xx),
6376 reinterpret_cast<double *
> (Bz),
6386 (*current_liboctave_error_handler)
6387 (
"SparseComplexMatrix::solve solve failed");
6402 sz = (
static_cast<double> (b_nc) - j) / b_nc
6404 sz = x_nz + (sz > 100 ? sz : 100);
6408 retval.
xdata (ii) = tmp;
6409 retval.
xridx (ii++) = i;
6412 retval.
xcidx (j+1) = ii;
6425 octave_unused_parameter (rcond);
6426 octave_unused_parameter (sing_handler);
6427 octave_unused_parameter (calc_cond);
6429 (*current_liboctave_error_handler)
6430 (
"support for UMFPACK was unavailable or disabled "
6431 "when liboctave was built");
6435 (*current_liboctave_error_handler) (
"incorrect matrix type");
6444 solve_singularity_handler sing_handler,
6445 bool calc_cond)
const
6453 if (nr != nc || nr != b.
rows ())
6454 (*current_liboctave_error_handler)
6455 (
"matrix dimension mismatch solution of linear equations");
6457 if (nr == 0 || b.
cols () == 0)
6462 int typ = mattype.
type ();
6467#if defined (HAVE_CHOLMOD)
6468 cholmod_common Common;
6469 cholmod_common *cm = &Common;
6473 cm->prefer_zomplex =
false;
6475 double spu = octave::sparse_params::get_key (
"spumoni");
6479 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6483 cm->print =
static_cast<int> (spu) + 2;
6484 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6488 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6489 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6491 cm->final_ll =
true;
6493 cholmod_sparse Astore;
6494 cholmod_sparse *
A = &Astore;
6504#if defined (OCTAVE_ENABLE_64)
6505 A->itype = CHOLMOD_LONG;
6507 A->itype = CHOLMOD_INT;
6509 A->dtype = CHOLMOD_DOUBLE;
6511 A->xtype = CHOLMOD_COMPLEX;
6515 cholmod_dense Bstore;
6516 cholmod_dense *
B = &Bstore;
6517 B->nrow = b.
rows ();
6518 B->ncol = b.
cols ();
6520 B->nzmax =
B->nrow *
B->ncol;
6521 B->dtype = CHOLMOD_DOUBLE;
6522 B->xtype = CHOLMOD_COMPLEX;
6541 if (is_singular (rcond) || octave::math::isnan (rcond))
6547 sing_handler (rcond);
6551 octave::warn_singular_matrix (rcond);
6563 retval.
xelem (i, j) =
static_cast<Complex *
> (X->x)[jr + i];
6569 static char blank_name[] =
" ";
6573 (*current_liboctave_warning_with_id_handler)
6574 (
"Octave:missing-dependency",
6575 "support for CHOLMOD was unavailable or disabled "
6576 "when liboctave was built");
6585#if defined (HAVE_UMFPACK)
6587 void *Numeric = factorize (err, rcond, Control, Info,
6588 sing_handler, calc_cond);
6593 Control (UMFPACK_IRSTEP) = 1;
6597 double *control = Control.
rwdata ();
6598 double *info = Info.
rwdata ();
6604 retval.
resize (b_nr, b_nc);
6611 octave::to_suitesparse_intptr (Ap),
6612 octave::to_suitesparse_intptr (Ai),
6613 reinterpret_cast<const double *
> (Ax),
6615 reinterpret_cast<double *
> (&Xx[iidx]),
6617 reinterpret_cast<const double *
> (&Bx[iidx]),
6618 nullptr, Numeric, control, info);
6625 (*current_liboctave_error_handler)
6626 (
"SparseComplexMatrix::solve solve failed");
6641 octave_unused_parameter (rcond);
6642 octave_unused_parameter (sing_handler);
6643 octave_unused_parameter (calc_cond);
6645 (*current_liboctave_error_handler)
6646 (
"support for UMFPACK was unavailable or disabled "
6647 "when liboctave was built");
6651 (*current_liboctave_error_handler) (
"incorrect matrix type");
6660 solve_singularity_handler sing_handler,
6661 bool calc_cond)
const
6669 if (nr != nc || nr != b.
rows ())
6670 (*current_liboctave_error_handler)
6671 (
"matrix dimension mismatch solution of linear equations");
6673 if (nr == 0 || b.
cols () == 0)
6678 int typ = mattype.
type ();
6683#if defined (HAVE_CHOLMOD)
6684 cholmod_common Common;
6685 cholmod_common *cm = &Common;
6689 cm->prefer_zomplex =
false;
6691 double spu = octave::sparse_params::get_key (
"spumoni");
6695 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6699 cm->print =
static_cast<int> (spu) + 2;
6700 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6704 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6705 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6707 cm->final_ll =
true;
6709 cholmod_sparse Astore;
6710 cholmod_sparse *
A = &Astore;
6720#if defined (OCTAVE_ENABLE_64)
6721 A->itype = CHOLMOD_LONG;
6723 A->itype = CHOLMOD_INT;
6725 A->dtype = CHOLMOD_DOUBLE;
6727 A->xtype = CHOLMOD_COMPLEX;
6731 cholmod_sparse Bstore;
6732 cholmod_sparse *
B = &Bstore;
6733 B->nrow = b.
rows ();
6734 B->ncol = b.
cols ();
6737 B->nzmax = b.
nnz ();
6741#if defined (OCTAVE_ENABLE_64)
6742 B->itype = CHOLMOD_LONG;
6744 B->itype = CHOLMOD_INT;
6746 B->dtype = CHOLMOD_DOUBLE;
6748 B->xtype = CHOLMOD_COMPLEX;
6767 if (is_singular (rcond) || octave::math::isnan (rcond))
6773 sing_handler (rcond);
6777 octave::warn_singular_matrix (rcond);
6782 cholmod_sparse *X =
CHOLMOD_NAME(spsolve) (CHOLMOD_A, L,
B, cm);
6785 (octave::from_size_t (X->nrow),
6786 octave::from_size_t (X->ncol),
6790 j <= static_cast<octave_idx_type> (X->ncol); j++)
6803 static char blank_name[] =
" ";
6807 (*current_liboctave_warning_with_id_handler)
6808 (
"Octave:missing-dependency",
6809 "support for CHOLMOD was unavailable or disabled "
6810 "when liboctave was built");
6819#if defined (HAVE_UMFPACK)
6821 void *Numeric = factorize (err, rcond, Control, Info,
6822 sing_handler, calc_cond);
6827 Control (UMFPACK_IRSTEP) = 1;
6831 double *control = Control.
rwdata ();
6832 double *info = Info.
rwdata ();
6847 retval.
xcidx (0) = 0;
6854 octave::to_suitesparse_intptr (Ap),
6855 octave::to_suitesparse_intptr (Ai),
6856 reinterpret_cast<const double *
> (Ax),
6858 reinterpret_cast<double *
> (Xx),
6860 reinterpret_cast<double *
> (Bx),
6861 nullptr, Numeric, control, info);
6868 (*current_liboctave_error_handler)
6869 (
"SparseComplexMatrix::solve solve failed");
6884 sz = (
static_cast<double> (b_nc) - j) / b_nc
6886 sz = x_nz + (sz > 100 ? sz : 100);
6890 retval.
xdata (ii) = tmp;
6891 retval.
xridx (ii++) = i;
6894 retval.
xcidx (j+1) = ii;
6899 rcond = Info (UMFPACK_RCOND);
6900 if (status == UMFPACK_WARNING_singular_matrix
6901 || is_singular (rcond) || octave::math::isnan (rcond))
6906 sing_handler (rcond);
6908 octave::warn_singular_matrix (rcond);
6919 octave_unused_parameter (rcond);
6920 octave_unused_parameter (sing_handler);
6921 octave_unused_parameter (calc_cond);
6923 (*current_liboctave_error_handler)
6924 (
"support for UMFPACK was unavailable or disabled "
6925 "when liboctave was built");
6929 (*current_liboctave_error_handler) (
"incorrect matrix type");
6940 return solve (mattype, b, info, rcond,
nullptr);
6948 return solve (mattype, b, info, rcond,
nullptr);
6955 return solve (mattype, b, info, rcond,
nullptr);
6961 solve_singularity_handler sing_handler,
6962 bool singular_fallback)
const
6965 int typ = mattype.
type (
false);
6968 typ = mattype.
type (*
this);
6971 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
6973 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
6975 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
6977 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
6980 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
6982 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
6984 (*current_liboctave_error_handler) (
"unknown matrix type");
7001 return solve (mattype, b, info, rcond,
nullptr);
7009 return solve (mattype, b, info, rcond,
nullptr);
7016 return solve (mattype, b, info, rcond,
nullptr);
7022 solve_singularity_handler sing_handler,
7023 bool singular_fallback)
const
7026 int typ = mattype.
type (
false);
7029 typ = mattype.
type (*
this);
7032 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
7034 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
7036 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
7038 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
7041 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
7043 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
7045 (*current_liboctave_error_handler) (
"unknown matrix type");
7062 return solve (mattype, b, info, rcond,
nullptr);
7070 return solve (mattype, b, info, rcond,
nullptr);
7077 return solve (mattype, b, info, rcond,
nullptr);
7083 solve_singularity_handler sing_handler,
7084 bool singular_fallback)
const
7087 int typ = mattype.
type (
false);
7090 typ = mattype.
type (*
this);
7093 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
7095 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
7097 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
7099 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
7102 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
7104 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
7106 (*current_liboctave_error_handler) (
"unknown matrix type");
7124 return solve (mattype, b, info, rcond,
nullptr);
7132 return solve (mattype, b, info, rcond,
nullptr);
7139 return solve (mattype, b, info, rcond,
nullptr);
7145 solve_singularity_handler sing_handler,
7146 bool singular_fallback)
const
7149 int typ = mattype.
type (
false);
7152 typ = mattype.
type (*
this);
7155 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
7157 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
7159 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
7161 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
7164 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
7166 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
7168 (*current_liboctave_error_handler) (
"unknown matrix type");
7184 return solve (mattype, b, info, rcond);
7192 return solve (mattype, b, info, rcond);
7199 return solve (mattype, b, info, rcond,
nullptr);
7205 solve_singularity_handler sing_handler)
const
7208 return solve (mattype, tmp, info, rcond,
7218 return solve (mattype, b, info, rcond,
nullptr);
7226 return solve (mattype, b, info, rcond,
nullptr);
7233 return solve (mattype, b, info, rcond,
nullptr);
7239 solve_singularity_handler sing_handler)
const
7242 return solve (mattype, tmp, info, rcond,
7251 return solve (b, info, rcond,
nullptr);
7258 return solve (b, info, rcond,
nullptr);
7263 double& rcond)
const
7265 return solve (b, info, rcond,
nullptr);
7271 solve_singularity_handler sing_handler)
const
7274 return solve (mattype, b, err, rcond, sing_handler);
7282 return solve (b, info, rcond,
nullptr);
7290 return solve (b, info, rcond,
nullptr);
7297 return solve (b, info, rcond,
nullptr);
7303 solve_singularity_handler sing_handler)
const
7306 return solve (mattype, b, err, rcond, sing_handler);
7314 return solve (b, info, rcond,
nullptr);
7321 return solve (b, info, rcond,
nullptr);
7327 solve_singularity_handler sing_handler)
const
7330 return solve (mattype, b, err, rcond, sing_handler);
7338 return solve (b, info, rcond,
nullptr);
7346 return solve (b, info, rcond,
nullptr);
7353 return solve (b, info, rcond,
nullptr);
7359 solve_singularity_handler sing_handler)
const
7362 return solve (mattype, b, err, rcond, sing_handler);
7369 return solve (b, info, rcond);
7376 return solve (b, info, rcond);
7381 double& rcond)
const
7383 return solve (b, info, rcond,
nullptr);
7389 solve_singularity_handler sing_handler)
const
7392 return solve (tmp, info, rcond,
7401 return solve (b, info, rcond,
nullptr);
7409 return solve (b, info, rcond,
nullptr);
7414 double& rcond)
const
7416 return solve (b, info, rcond,
nullptr);
7422 solve_singularity_handler sing_handler)
const
7425 return solve (tmp, info, rcond,
7434 octave::err_nan_to_logical_conversion ();
7450 if (jj <
cidx (i+1) &&
ridx (jj) == j)
7498 if (octave::math::isnan (val))
7513 if (octave::math::isinf (val) || octave::math::isnan (val))
7540 max_val = std::real (
data (0));
7541 min_val = std::real (
data (0));
7547 double r_val = val.real ();
7548 double i_val = val.imag ();
7550 if (r_val > max_val)
7553 if (i_val > max_val)
7556 if (r_val < min_val)
7559 if (i_val < min_val)
7562 if (octave::math::round (r_val) != r_val
7563 || octave::math::round (i_val) != i_val)
7573 return test_any (octave::too_large_for_float);
7604 if (! octave::math::isnan (data (j))) \
7625 else if ((
rows () == 1 && dim == -1) || dim == 1)
7630 Complex d = data (i); \
7631 if (nanflag && octave::math::isnan (d)) \
7632 tmp[ridx (i)] *= 1.0; \
7637 Complex d = data (i); \
7638 if (nanflag && octave::math::isnan (d)) \
7645 (
cidx (j+1) -
cidx (j) < nr ? 0.0 : 1.0), 1.0);
7656 Complex d = data (i); \
7657 if (nanflag && octave::math::isnan (d)) \
7658 tmp[ridx (i)] += 0.0; \
7663 Complex d = data (i); \
7664 if (nanflag && octave::math::isnan (d)) \
7689#define EXPR r.data (nel++) = data (i) * conj (data (i));
7692 Complex d = data (i); \
7693 if (nanflag && octave::math::isnan (d)) \
7694 tmp[ridx (i)] += 0.0; \
7696 tmp[ridx (i)] += d * conj (d)
7699 Complex d = data (i); \
7700 if (nanflag && octave::math::isnan (d)) \
7703 tmp[j] += d * conj (d)
7727 retval.
data (i) = std::abs (
data (i));
7752 os << a.
ridx (i) + 1 <<
' ' << j + 1 <<
' ';
7753 octave::write_value<Complex> (os, a.
data (i));
7766 return read_sparse_matrix<elt_type> (is, a, octave::read_value<Complex>);
7851 return do_mul_dm_sm<SparseComplexMatrix> (
d, a);
7856 return do_mul_sm_dm<SparseComplexMatrix> (a,
d);
7862 return do_mul_dm_sm<SparseComplexMatrix> (
d, a);
7867 return do_mul_sm_dm<SparseComplexMatrix> (a,
d);
7873 return do_mul_dm_sm<SparseComplexMatrix> (
d, a);
7878 return do_mul_sm_dm<SparseComplexMatrix> (a,
d);
7884 return do_add_dm_sm<SparseComplexMatrix> (
d, a);
7889 return do_add_dm_sm<SparseComplexMatrix> (
d, a);
7894 return do_add_dm_sm<SparseComplexMatrix> (
d, a);
7899 return do_add_sm_dm<SparseComplexMatrix> (a,
d);
7904 return do_add_sm_dm<SparseComplexMatrix> (a,
d);
7909 return do_add_sm_dm<SparseComplexMatrix> (a,
d);
7915 return do_sub_dm_sm<SparseComplexMatrix> (
d, a);
7920 return do_sub_dm_sm<SparseComplexMatrix> (
d, a);
7925 return do_sub_dm_sm<SparseComplexMatrix> (
d, a);
7930 return do_sub_sm_dm<SparseComplexMatrix> (a,
d);
7935 return do_sub_sm_dm<SparseComplexMatrix> (a,
d);
7940 return do_sub_sm_dm<SparseComplexMatrix> (a,
d);
7959#define EMPTY_RETURN_CHECK(T) \
7960 if (nr == 0 || nc == 0) \
7966 const bool nanflag =
true;
7967 const bool realabs =
false;
7968 return min (c, m, nanflag, realabs);
7974 const bool realabs =
false;
7975 return min (c, m, nanflag, realabs);
7980 const bool nanflag,
const bool realabs)
7997 result.
data (i) = octave::math::min (c, m.
data (i), nanflag, realabs);
8006 const bool nanflag =
true;
8007 const bool realabs =
false;
8008 return min (c, m, nanflag, realabs);
8014 const bool realabs =
false;
8015 return min (c, m, nanflag, realabs);
8020 const bool nanflag,
const bool realabs)
8022 return min (c, m, nanflag, realabs);
8028 const bool nanflag =
true;
8029 const bool realabs =
false;
8030 return min (a, b, nanflag, realabs);
8037 const bool realabs =
false;
8038 return min (a, b, nanflag, realabs);
8043 const bool nanflag,
const bool realabs)
8052 if (a_nr == b_nr && a_nc == b_nc)
8062 bool ja_lt_max = ja < ja_max;
8066 bool jb_lt_max = jb < jb_max;
8068 while (ja_lt_max || jb_lt_max)
8071 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
8082 ja_lt_max= ja < ja_max;
8084 else if ((! ja_lt_max)
8085 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
8096 jb_lt_max= jb < jb_max;
8109 ja_lt_max= ja < ja_max;
8111 jb_lt_max= jb < jb_max;
8121 if (a_nr == 0 && (b_nr == 0 || b_nr == 1))
8123 if (a_nc == 1 || b_nc == 1 || a_nc == b_nc)
8124 r.
resize (a_nr, std::max (a_nc, b_nc));
8126 octave::err_nonconformant (
"min", a_nr, a_nc, b_nr, b_nc);
8128 else if (a_nc == 0 && (b_nc == 0 || b_nc == 1))
8130 if (a_nr == 1 || b_nr == 1 || a_nr == b_nr)
8131 r.
resize (std::max (a_nr, b_nr), a_nc);
8133 octave::err_nonconformant (
"min", a_nr, a_nc, b_nr, b_nc);
8135 else if (b_nr == 0 && (a_nr == 0 || a_nr == 1))
8137 if (b_nc == 1 || a_nc == 1 || b_nc == a_nc)
8138 r.
resize (b_nr, std::max (a_nc, b_nc));
8140 octave::err_nonconformant (
"min", a_nr, a_nc, b_nr, b_nc);
8142 else if (b_nc == 0 && (a_nc == 0 || a_nc == 1))
8144 if (b_nr == 1 || a_nr == 1 || b_nr == a_nr)
8145 r.
resize (std::max (a_nr, b_nr), b_nc);
8147 octave::err_nonconformant (
"min", a_nr, a_nc, b_nr, b_nc);
8150 octave::err_nonconformant (
"min", a_nr, a_nc, b_nr, b_nc);
8159 const bool nanflag =
true;
8160 const bool realabs =
false;
8161 return max (c, m, nanflag, realabs);
8167 const bool realabs =
false;
8168 return max (c, m, nanflag, realabs);
8183 if (octave::math::max (c, 0., nanflag, realabs) != 0.)
8188 result.
xdata (m.
ridx (i) + j * nr) = octave::math::max
8189 (c, m.
data (i), nanflag, realabs);
8200 const bool nanflag =
true;
8201 const bool realabs =
false;
8202 return max (c, m, nanflag, realabs);
8208 const bool realabs =
false;
8209 return max (c, m, nanflag, realabs);
8216 return max (c, m, nanflag, realabs);
8222 const bool nanflag =
true;
8223 const bool realabs =
false;
8224 return max (a, b, nanflag, realabs);
8231 const bool realabs =
false;
8232 return max (a, b, nanflag, realabs);
8237 const bool nanflag,
const bool realabs)
8246 if (a_nr == b_nr && a_nc == b_nc)
8256 bool ja_lt_max = ja < ja_max;
8260 bool jb_lt_max = jb < jb_max;
8262 while (ja_lt_max || jb_lt_max)
8265 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
8276 ja_lt_max= ja < ja_max;
8278 else if ((! ja_lt_max)
8279 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
8290 jb_lt_max= jb < jb_max;
8303 ja_lt_max= ja < ja_max;
8305 jb_lt_max= jb < jb_max;
8315 if (a_nr == 0 && (b_nr == 0 || b_nr == 1))
8317 if (a_nc == 1 || b_nc == 1 || a_nc == b_nc)
8318 r.
resize (a_nr, std::max (a_nc, b_nc));
8320 octave::err_nonconformant (
"max", a_nr, a_nc, b_nr, b_nc);
8322 else if (a_nc == 0 && (b_nc == 0 || b_nc == 1))
8324 if (a_nr == 1 || b_nr == 1 || a_nr == b_nr)
8325 r.
resize (std::max (a_nr, b_nr), a_nc);
8327 octave::err_nonconformant (
"max", a_nr, a_nc, b_nr, b_nc);
8329 else if (b_nr == 0 && (a_nr == 0 || a_nr == 1))
8331 if (b_nc == 1 || a_nc == 1 || b_nc == a_nc)
8332 r.
resize (b_nr, std::max (a_nc, b_nc));
8334 octave::err_nonconformant (
"max", a_nr, a_nc, b_nr, b_nc);
8336 else if (b_nc == 0 && (a_nc == 0 || a_nc == 1))
8338 if (b_nr == 1 || a_nr == 1 || b_nr == a_nr)
8339 r.
resize (std::max (a_nr, b_nr), b_nc);
8341 octave::err_nonconformant (
"max", a_nr, a_nc, b_nr, b_nc);
8344 octave::err_nonconformant (
"max", a_nr, a_nc, b_nr, b_nc);
SparseComplexMatrix max(const Complex &c, const SparseComplexMatrix &m)
SparseComplexMatrix operator-(const ComplexDiagMatrix &d, const SparseMatrix &a)
SparseComplexMatrix conj(const SparseComplexMatrix &a)
ComplexMatrix mul_trans(const ComplexMatrix &m, const SparseComplexMatrix &a)
SparseComplexMatrix operator*(const SparseComplexMatrix &m, const SparseMatrix &a)
#define EMPTY_RETURN_CHECK(T)
ComplexMatrix mul_herm(const ComplexMatrix &m, const SparseComplexMatrix &a)
SparseComplexMatrix min(const Complex &c, const SparseComplexMatrix &m)
std::ostream & operator<<(std::ostream &os, const SparseComplexMatrix &a)
SparseComplexMatrix operator+(const ComplexDiagMatrix &d, const SparseMatrix &a)
ComplexMatrix trans_mul(const SparseComplexMatrix &m, const ComplexMatrix &a)
std::istream & operator>>(std::istream &is, SparseComplexMatrix &a)
ComplexMatrix herm_mul(const SparseComplexMatrix &m, const ComplexMatrix &a)
base_det< Complex > ComplexDET
#define SPARSE_SSM_BOOL_OPS(S, M)
#define SPARSE_ANY_ALL_HEADER
#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_CUMSUM(RET_TYPE, ELT_TYPE, FCN, NAN_EXPR, EXPR)
#define SPARSE_XSUM_REDUCTION_OP(RET_TYPE, EL_TYPE)
#define SPARSE_SUMSQ_HEADER(RET_TYPE, EXPR)
#define SPARSE_SSM_CMP_OPS(S, M)
#define SPARSE_SMSM_CMP_OPS(M1, M2)
#define SPARSE_SMS_CMP_OPS(M, S)
#define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE)
#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)
N Dimensional Array with copy-on-write semantics.
T & xelem(octave_idx_type n)
Size of the specified dimension.
T & elem(octave_idx_type n)
Size of the specified dimension.
octave_idx_type rows() const
void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
octave_idx_type cols() const
const T * data() const
Size of the specified dimension.
T * rwdata()
Size of the specified dimension.
friend ColumnVector imag(const ComplexColumnVector &a)
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
ComplexColumnVector column(octave_idx_type i) const
octave_idx_type length() const
octave_idx_type cols() const
MSparse< T > diag(octave_idx_type k=0) const
MSparse< T > & insert(const Sparse< T > &a, octave_idx_type r, octave_idx_type c)
MSparse< T > permute(const Array< octave_idx_type > &vec, bool inv=false) const
MSparse< T > reshape(const dim_vector &new_dims) const
MSparse< T > squeeze() const
MSparse< T > ipermute(const Array< octave_idx_type > &vec) const
void mark_as_unsymmetric()
octave_idx_type * triangular_perm() const
MatrixType transpose() const
int type(bool quiet=true)
void mark_as_rectangular()
SparseComplexMatrix cumprod(int dim=-1, bool nanflag=false) const
ComplexRowVector row(octave_idx_type i) const
SparseComplexMatrix cumsum(int dim=-1, bool nanflag=false) const
SparseComplexMatrix xsum(int dim=-1, bool nanflag=false) const
bool too_large_for_float() const
bool any_element_is_inf_or_nan() const
friend SparseComplexMatrix conj(const SparseComplexMatrix &a)
SparseComplexMatrix hermitian() const
SparseComplexMatrix min(int dim=-1, bool nanflag=true, bool realabs=false) const
SparseComplexMatrix diag(octave_idx_type k=0) const
ComplexColumnVector column(octave_idx_type i) const
ComplexMatrix matrix_value() const
bool all_elements_are_real() const
SparseComplexMatrix & insert(const SparseComplexMatrix &a, octave_idx_type r, octave_idx_type c)
bool any_element_is_nan() const
SparseComplexMatrix inverse() const
bool all_integers(double &max_val, double &min_val) const
SparseComplexMatrix reshape(const dim_vector &new_dims) const
SparseComplexMatrix permute(const Array< octave_idx_type > &vec, bool inv=false) const
bool operator!=(const SparseComplexMatrix &a) const
SparseComplexMatrix max(int dim=-1, bool nanflag=true, bool realabs=false) const
SparseBoolMatrix any(int dim=-1) const
SparseComplexMatrix transpose() const
SparseComplexMatrix ipermute(const Array< octave_idx_type > &vec) const
bool operator==(const SparseComplexMatrix &a) const
SparseComplexMatrix concat(const SparseComplexMatrix &rb, const Array< octave_idx_type > &ra_idx)
SparseComplexMatrix squeeze() const
ComplexDET determinant() const
SparseBoolMatrix operator!() const
ComplexMatrix solve(MatrixType &mattype, const Matrix &b) const
SparseComplexMatrix prod(int dim=-1, bool nanflag=false) const
SparseComplexMatrix sumsq(int dim=-1, bool nanflag=false) const
SparseComplexMatrix sum(int dim=-1, bool nanflag=false) const
SparseBoolMatrix all(int dim=-1) const
octave_idx_type cols() const
Array< T > array_value() const
octave_idx_type columns() const
void resize(octave_idx_type r, octave_idx_type c)
bool test_any(F fcn) const
T & elem(octave_idx_type n)
Sparse< T, Alloc > maybe_compress(bool remove_zeros=false)
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 * xcidx()
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
ColumnVector real(const ComplexColumnVector &a)
ColumnVector imag(const ComplexColumnVector &a)
#define F77_DBLE_CMPLX_ARG(x)
#define F77_XFCN(f, F, args)
octave_f77_int_type F77_INT
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
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
bool mx_inline_all_real(std::size_t n, const std::complex< T > *x)
std::complex< double > Complex
#define liboctave_panic_unless(cond)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
#define UMFPACK_ZNAME(name)
#define CHOLMOD_NAME(name)
const octave_base_value const Array< octave_idx_type > & ra_idx
template int8_t abs(int8_t)
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
template ComplexMatrix dmsolve< ComplexMatrix, SparseComplexMatrix, ComplexMatrix >(const SparseComplexMatrix &, const ComplexMatrix &, octave_idx_type &)
template ComplexMatrix dmsolve< ComplexMatrix, SparseComplexMatrix, Matrix >(const SparseComplexMatrix &, const Matrix &, octave_idx_type &)
RT dmsolve(const ST &a, const T &b, octave_idx_type &info)
template SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseComplexMatrix, SparseMatrix >(const SparseComplexMatrix &, const SparseMatrix &, octave_idx_type &)
void SparseCholError(int status, char *file, int line, char *message)
int SparseCholPrint(const char *fmt,...)