133 if (nr == nc && nr > 0)
186 return max (dummy_idx, dim, nanflag, realabs);
191 bool nanflag,
bool realabs)
const
198 if (dim >= dv.
ndims ())
211 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
219 double tmp_max = octave::numeric_limits<double>::NaN ();
223 if (
ridx (i) != idx_j)
234 double tmp =
data (i);
236 if (octave::math::isnan (tmp) && ! nanflag)
242 else if (octave::math::isnan (tmp))
244 else if (octave::math::isnan (tmp_max) || tmp > tmp_max)
252 if (octave::math::isnan (tmp_max) && ! nanflag)
253 idx_arg.
elem (j) = idx_j;
255 idx_arg.
elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
264 double tmp_max = octave::numeric_limits<double>::NaN ();
265 double abs_max = octave::numeric_limits<double>::NaN ();
269 if (
ridx (i) != idx_j)
282 double tmp =
data (i);
284 if (octave::math::isnan (tmp) && ! nanflag)
290 else if (octave::math::isnan (tmp))
293 double abs_tmp = std::abs (tmp);
294 if (octave::math::isnan (abs_max) || abs_tmp > abs_max)
300 else if (abs_tmp == abs_max && tmp > tmp_max)
309 if (octave::math::isnan (tmp_max) && ! nanflag)
310 idx_arg.
elem (j) = idx_j;
312 idx_arg.
elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
321 result.
xcidx (0) = 0;
324 double tmp =
elem (idx_arg(j), j);
327 result.
xdata (ii) = tmp;
328 result.
xridx (ii++) = 0;
330 result.
xcidx (j+1) = ii;
338 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
348 if (found[
ridx (i)] == -j)
349 found[
ridx (i)] = -j - 1;
352 if (found[i] > -nc && found[i] < 0)
353 idx_arg.
elem (i) = -found[i];
363 double tmp =
data (i);
365 if (octave::math::isnan (tmp) && ! nanflag)
367 idx_arg.
elem (ir) = j;
370 else if (octave::math::isnan (tmp))
372 else if (ix == -1 || octave::math::isnan (
elem (ir, ix)) ||
374 idx_arg.
elem (ir) = j;
386 double tmp =
data (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 std::abs (tmp) > std::abs (
elem (ir, ix)))
397 idx_arg.
elem (ir) = j;
398 else if (std::abs (tmp) == std::abs (
elem (ir, ix)) &&
400 idx_arg.
elem (ir) = j;
407 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
413 result.
xcidx (0) = 0;
414 result.
xcidx (1) = nel;
417 if (idx_arg(j) == -1)
420 result.
xdata (ii) = octave::numeric_limits<double>::NaN ();
421 result.
xridx (ii++) = j;
425 double tmp =
elem (j, idx_arg(j));
428 result.
xdata (ii) = tmp;
429 result.
xridx (ii++) = j;
442 return min (dummy_idx, dim, nanflag, realabs);
447 bool nanflag,
bool realabs)
const
454 if (dim >= dv.
ndims ())
467 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
475 double tmp_max = octave::numeric_limits<double>::NaN ();
479 if (
ridx (i) != idx_j)
490 double tmp =
data (i);
492 if (octave::math::isnan (tmp) && ! nanflag)
498 else if (octave::math::isnan (tmp))
500 else if (octave::math::isnan (tmp_max) || tmp < tmp_max)
508 if (octave::math::isnan (tmp_max) && ! nanflag)
509 idx_arg.
elem (j) = idx_j;
511 idx_arg.
elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
520 double tmp_max = octave::numeric_limits<double>::NaN ();
521 double abs_max = octave::numeric_limits<double>::NaN ();
525 if (
ridx (i) != idx_j)
538 double tmp =
data (i);
540 if (octave::math::isnan (tmp) && ! nanflag)
546 else if (octave::math::isnan (tmp))
549 double abs_tmp = std::abs (tmp);
550 if (octave::math::isnan (abs_max) || abs_tmp < abs_max)
556 else if (abs_tmp == abs_max && tmp < tmp_max)
565 if (octave::math::isnan (tmp_max) && ! nanflag)
566 idx_arg.
elem (j) = idx_j;
568 idx_arg.
elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
577 result.
xcidx (0) = 0;
580 double tmp =
elem (idx_arg(j), j);
583 result.
xdata (ii) = tmp;
584 result.
xridx (ii++) = 0;
586 result.
xcidx (j+1) = ii;
594 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
604 if (found[
ridx (i)] == -j)
605 found[
ridx (i)] = -j - 1;
608 if (found[i] > -nc && found[i] < 0)
609 idx_arg.
elem (i) = -found[i];
619 double tmp =
data (i);
621 if (octave::math::isnan (tmp) && ! nanflag)
623 idx_arg.
elem (ir) = j;
626 else if (octave::math::isnan (tmp))
628 else if (ix == -1 || octave::math::isnan (
elem (ir, ix)) ||
630 idx_arg.
elem (ir) = j;
642 double tmp =
data (i);
644 if (octave::math::isnan (tmp) && ! nanflag)
646 idx_arg.
elem (ir) = j;
649 else if (octave::math::isnan (tmp))
651 else if (ix == -1 || octave::math::isnan (
elem (ir, ix)) ||
652 std::abs (tmp) < std::abs (
elem (ir, ix)))
653 idx_arg.
elem (ir) = j;
654 else if (std::abs (tmp) == std::abs (
elem (ir, ix)) &&
656 idx_arg.
elem (ir) = j;
663 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
669 result.
xcidx (0) = 0;
670 result.
xcidx (1) = nel;
673 if (idx_arg(j) == -1)
676 result.
xdata (ii) = octave::numeric_limits<double>::NaN ();
677 result.
xridx (ii++) = j;
681 double tmp =
elem (j, idx_arg(j));
684 result.
xdata (ii) = tmp;
685 result.
xridx (ii++) = j;
720 retval(j) =
data (k);
745 if (rb.
rows () > 0 && rb.
cols () > 0)
755 if (rb.
rows () > 0 && rb.
cols () > 0)
773 r.
data (i) = std::real (a.
data (i));
794 r.
data (i) = std::imag (a.
data (i));
812is_singular (
const double rcond)
814 return (std::abs (rcond) <= std::numeric_limits<double>::epsilon ());
823 return inverse (mattype, info, rcond, 0, 0);
831 return inverse (mattype, info, rcond, 0, 0);
838 return inverse (mattype, info, rcond, 0, 0);
843 double& rcond,
const bool,
844 const bool calccond)
const
852 if (nr == 0 || nc == 0 || nr != nc)
853 (*current_liboctave_error_handler) (
"inverse requires square matrix");
856 int typ = mattype.
type ();
860 (*current_liboctave_error_handler) (
"incorrect matrix type");
868 double *v = retval.
data ();
873 double dmin = octave::numeric_limits<double>::Inf ();
876 double tmp = fabs (v[i]);
893 double& rcond,
const bool,
894 const bool calccond)
const
902 if (nr == 0 || nc == 0 || nr != nc)
903 (*current_liboctave_error_handler) (
"inverse requires square matrix");
906 int typ = mattype.
type ();
911 (*current_liboctave_error_handler) (
"incorrect matrix type");
914 double ainvnorm = 0.;
923 atmp += fabs (
data (i));
948 retval.
xcidx (i) = cx;
949 retval.
xridx (cx) = i;
950 retval.
xdata (cx) = 1.0;
963 (*current_liboctave_error_handler) (
"division by zero");
968 rpX = retval.
xridx (colXp);
977 v -= retval.
xdata (colXp) *
data (colUp);
982 while (rpX < j && rpU < j && colXp < cx && colUp < nz);
986 colUp =
cidx (j+1) - 1;
989 double pivot =
data (colUp);
990 if (pivot == 0. ||
ridx (colUp) != j)
991 (*current_liboctave_error_handler) (
"division by zero");
1001 retval.
xridx (cx) = j;
1002 retval.
xdata (cx) = v / pivot;
1010 colUp =
cidx (i+1) - 1;
1013 double pivot =
data (colUp);
1014 if (pivot == 0. ||
ridx (colUp) != i)
1015 (*current_liboctave_error_handler) (
"division by zero");
1019 retval.
xdata (j) /= pivot;
1021 retval.
xcidx (nr) = cx;
1066 k <
cidx (jidx+1); k++)
1079 (*current_liboctave_error_handler) (
"division by zero");
1081 work[j] = v / pivot;
1087 colUp =
cidx (perm[iidx]+1) - 1;
1089 colUp =
cidx (perm[iidx]);
1091 double pivot =
data (colUp);
1093 (*current_liboctave_error_handler) (
"division by zero");
1106 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
1110 retval.
xcidx (i) = cx;
1114 retval.
xridx (cx) = j;
1115 retval.
xdata (cx++) = work[j];
1119 retval.
xcidx (nr) = cx;
1130 i < retval.
cidx (j+1); i++)
1131 atmp += fabs (retval.
data (i));
1132 if (atmp > ainvnorm)
1136 rcond = 1. / ainvnorm / anorm;
1144 double& rcond,
bool,
bool calc_cond)
const
1148 (*current_liboctave_error_handler)
1149 (
"inverse of the null matrix not defined");
1152 int typ = mattype.
type (
false);
1156 typ = mattype.
type (*
this);
1159 ret = dinverse (mattype, info, rcond,
true, calc_cond);
1161 ret = tinverse (mattype, info, rcond,
true, calc_cond).
transpose ();
1165 ret =
transpose ().tinverse (newtype, info, rcond,
true, calc_cond);
1172 octave::math::sparse_chol<SparseMatrix> fact (*
this, info,
false);
1173 rcond = fact.rcond ();
1179 info, rcond2,
true,
false);
1197 octave::math::sparse_lu<SparseMatrix> fact (*
this,
1200 rcond = fact.rcond ();
1207 octave::numeric_limits<double>::Inf ());
1208 std::copy_n (
ridx (), nz, ret.
xridx ());
1216 info, rcond2,
true,
false);
1217 SparseMatrix InvU = fact.U ().tinverse (tmp_typ, info, rcond2,
1219 ret = fact.Pc ().
transpose () * InvU * InvL * fact.Pr ();
1246#if defined (HAVE_UMFPACK)
1251 if (nr == 0 || nc == 0 || nr != nc)
1260 Matrix Control (UMFPACK_CONTROL, 1);
1261 double *control = Control.
rwdata ();
1264 double tmp = octave::sparse_params::get_key (
"spumoni");
1265 if (! octave::math::isnan (tmp))
1266 Control (UMFPACK_PRL) = tmp;
1268 tmp = octave::sparse_params::get_key (
"piv_tol");
1269 if (! octave::math::isnan (tmp))
1271 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
1272 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
1276 tmp = octave::sparse_params::get_key (
"autoamd");
1277 if (! octave::math::isnan (tmp))
1278 Control (UMFPACK_FIXQ) = tmp;
1281 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1287 const double *Ax =
data ();
1290 octave::to_suitesparse_intptr (Ap),
1291 octave::to_suitesparse_intptr (Ai),
1295 Matrix Info (1, UMFPACK_INFO);
1296 double *info = Info.
rwdata ();
1298 octave::to_suitesparse_intptr (Ap),
1299 octave::to_suitesparse_intptr (Ai),
1300 Ax,
nullptr, &Symbolic, control, info);
1309 (*current_liboctave_error_handler)
1310 (
"SparseMatrix::determinant symbolic factorization failed");
1317 status =
UMFPACK_DNAME (numeric) (octave::to_suitesparse_intptr (Ap),
1318 octave::to_suitesparse_intptr (Ai),
1320 &Numeric, control, info);
1323 rcond = Info (UMFPACK_RCOND);
1331 (*current_liboctave_error_handler)
1332 (
"SparseMatrix::determinant numeric factorization failed");
1340 status =
UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric,
1348 (*current_liboctave_error_handler)
1349 (
"SparseMatrix::determinant error calculating determinant");
1352 retval =
DET (c10, e10, 10);
1361 octave_unused_parameter (err);
1362 octave_unused_parameter (rcond);
1364 (*current_liboctave_error_handler)
1365 (
"support for UMFPACK was unavailable or disabled "
1366 "when liboctave was built");
1376 double& rcond, solve_singularity_handler,
1377 bool calc_cond)
const
1386 if (nr != b.
rows ())
1388 (
"matrix dimension mismatch solution of linear equations");
1390 if (nr == 0 || nc == 0 || b.
cols () == 0)
1395 int typ = mattype.
type ();
1399 (*current_liboctave_error_handler) (
"incorrect matrix type");
1405 retval(i, j) = b(i, j) /
data (i);
1410 retval(k, j) = b(
ridx (i), j) /
data (i);
1415 double dmin = octave::numeric_limits<double>::Inf ();
1418 double tmp = fabs (
data (i));
1424 rcond = dmin / dmax;
1436 solve_singularity_handler,
bool calc_cond)
const
1445 if (nr != b.
rows ())
1447 (
"matrix dimension mismatch solution of linear equations");
1449 if (nr == 0 || nc == 0 || b.
cols () == 0)
1454 int typ = mattype.
type ();
1458 (*current_liboctave_error_handler) (
"incorrect matrix type");
1464 retval.
xcidx (0) = 0;
1471 if (b.
ridx (i) >= nm)
1476 retval.
xcidx (j+1) = ii;
1486 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1494 retval.
xridx (ii) = l;
1498 retval.
xcidx (j+1) = ii;
1504 double dmin = octave::numeric_limits<double>::Inf ();
1507 double tmp = fabs (
data (i));
1513 rcond = dmin / dmax;
1525 solve_singularity_handler,
bool calc_cond)
const
1534 if (nr != b.
rows ())
1536 (
"matrix dimension mismatch solution of linear equations");
1538 if (nr == 0 || nc == 0 || b.
cols () == 0)
1543 int typ = mattype.
type ();
1547 (*current_liboctave_error_handler) (
"incorrect matrix type");
1553 retval(i, j) = b(i, j) /
data (i);
1558 retval(k, j) = b(
ridx (i), j) /
data (i);
1563 double dmin = octave::numeric_limits<double>::Inf ();
1566 double tmp = fabs (
data (i));
1572 rcond = dmin / dmax;
1584 solve_singularity_handler,
bool calc_cond)
const
1593 if (nr != b.
rows ())
1595 (
"matrix dimension mismatch solution of linear equations");
1597 if (nr == 0 || nc == 0 || b.
cols () == 0)
1602 int typ = mattype.
type ();
1606 (*current_liboctave_error_handler) (
"incorrect matrix type");
1612 retval.
xcidx (0) = 0;
1619 if (b.
ridx (i) >= nm)
1624 retval.
xcidx (j+1) = ii;
1634 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1642 retval.
xridx (ii) = l;
1646 retval.
xcidx (j+1) = ii;
1652 double dmin = octave::numeric_limits<double>::Inf ();
1655 double tmp = fabs (
data (i));
1661 rcond = dmin / dmax;
1673 solve_singularity_handler sing_handler,
1674 bool calc_cond)
const
1683 if (nr != b.
rows ())
1684 (*current_liboctave_error_handler)
1685 (
"matrix dimension mismatch solution of linear equations");
1687 if (nr == 0 || nc == 0 || b.
cols () == 0)
1692 int typ = mattype.
type ();
1696 (*current_liboctave_error_handler) (
"incorrect matrix type");
1699 double ainvnorm = 0.;
1710 atmp += fabs (
data (i));
1718 retval.
resize (nc, b_nc);
1739 goto triangular_error;
1742 double tmp = work[k] /
data (
cidx (kidx+1)-1);
1745 i <
cidx (kidx+1)-1; i++)
1748 work[iidx] = work[iidx] - tmp *
data (i);
1754 retval.
xelem (perm[i], j) = work[i];
1773 double tmp = work[k] /
data (
cidx (iidx+1)-1);
1776 i <
cidx (iidx+1)-1; i++)
1779 work[idx2] = work[idx2] - tmp *
data (i);
1786 atmp += fabs (work[i]);
1789 if (atmp > ainvnorm)
1792 rcond = 1. / ainvnorm / anorm;
1798 retval.
resize (nc, b_nc);
1815 goto triangular_error;
1818 double tmp = work[k] /
data (
cidx (k+1)-1);
1823 work[iidx] = work[iidx] - tmp *
data (i);
1829 retval.
xelem (i, j) = work[i];
1846 double tmp = work[k] /
data (
cidx (k+1)-1);
1851 work[iidx] = work[iidx] - tmp *
data (i);
1858 atmp += fabs (work[i]);
1861 if (atmp > ainvnorm)
1864 rcond = 1. / ainvnorm / anorm;
1873 sing_handler (rcond);
1877 octave::warn_singular_matrix (rcond);
1880 if (is_singular (rcond) || octave::math::isnan (rcond))
1886 sing_handler (rcond);
1890 octave::warn_singular_matrix (rcond);
1900 solve_singularity_handler sing_handler,
1901 bool calc_cond)
const
1910 if (nr != b.
rows ())
1911 (*current_liboctave_error_handler)
1912 (
"matrix dimension mismatch solution of linear equations");
1914 if (nr == 0 || nc == 0 || b.
cols () == 0)
1919 int typ = mattype.
type ();
1923 (*current_liboctave_error_handler) (
"incorrect matrix type");
1926 double ainvnorm = 0.;
1936 atmp += fabs (
data (i));
1945 retval.
xcidx (0) = 0;
1975 goto triangular_error;
1978 double tmp = work[k] /
data (
cidx (kidx+1)-1);
1981 i <
cidx (kidx+1)-1; i++)
1984 work[iidx] = work[iidx] - tmp *
data (i);
1996 if (ii + new_nnz > x_nz)
2005 if (work[rperm[i]] != 0.)
2007 retval.
xridx (ii) = i;
2008 retval.
xdata (ii++) = work[rperm[i]];
2010 retval.
xcidx (j+1) = ii;
2031 double tmp = work[k] /
data (
cidx (iidx+1)-1);
2034 i <
cidx (iidx+1)-1; i++)
2037 work[idx2] = work[idx2] - tmp *
data (i);
2044 atmp += fabs (work[i]);
2047 if (atmp > ainvnorm)
2050 rcond = 1. / ainvnorm / anorm;
2072 goto triangular_error;
2075 double tmp = work[k] /
data (
cidx (k+1)-1);
2080 work[iidx] = work[iidx] - tmp *
data (i);
2092 if (ii + new_nnz > x_nz)
2103 retval.
xridx (ii) = i;
2104 retval.
xdata (ii++) = work[i];
2106 retval.
xcidx (j+1) = ii;
2125 double tmp = work[k] /
data (
cidx (k+1)-1);
2128 i <
cidx (k+1)-1; i++)
2131 work[iidx] = work[iidx] - tmp *
data (i);
2138 atmp += fabs (work[i]);
2141 if (atmp > ainvnorm)
2144 rcond = 1. / ainvnorm / anorm;
2153 sing_handler (rcond);
2157 octave::warn_singular_matrix (rcond);
2160 if (is_singular (rcond) || octave::math::isnan (rcond))
2166 sing_handler (rcond);
2170 octave::warn_singular_matrix (rcond);
2179 solve_singularity_handler sing_handler,
2180 bool calc_cond)
const
2189 if (nr != b.
rows ())
2190 (*current_liboctave_error_handler)
2191 (
"matrix dimension mismatch solution of linear equations");
2193 if (nr == 0 || nc == 0 || b.
cols () == 0)
2198 int typ = mattype.
type ();
2202 (*current_liboctave_error_handler) (
"incorrect matrix type");
2205 double ainvnorm = 0.;
2216 atmp += fabs (
data (i));
2224 retval.
resize (nc, b_nc);
2245 goto triangular_error;
2251 i <
cidx (kidx+1)-1; i++)
2254 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2260 retval.
xelem (perm[i], j) = cwork[i];
2280 double tmp = work[k] /
data (
cidx (iidx+1)-1);
2283 i <
cidx (iidx+1)-1; i++)
2286 work[idx2] = work[idx2] - tmp *
data (i);
2293 atmp += fabs (work[i]);
2296 if (atmp > ainvnorm)
2299 rcond = 1. / ainvnorm / anorm;
2305 retval.
resize (nc, b_nc);
2322 goto triangular_error;
2330 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2336 retval.
xelem (i, j) = cwork[i];
2354 double tmp = work[k] /
data (
cidx (k+1)-1);
2357 i <
cidx (k+1)-1; i++)
2360 work[iidx] = work[iidx] - tmp *
data (i);
2367 atmp += fabs (work[i]);
2370 if (atmp > ainvnorm)
2373 rcond = 1. / ainvnorm / anorm;
2382 sing_handler (rcond);
2386 octave::warn_singular_matrix (rcond);
2389 if (is_singular (rcond) || octave::math::isnan (rcond))
2395 sing_handler (rcond);
2399 octave::warn_singular_matrix (rcond);
2409 solve_singularity_handler sing_handler,
2410 bool calc_cond)
const
2419 if (nr != b.
rows ())
2420 (*current_liboctave_error_handler)
2421 (
"matrix dimension mismatch solution of linear equations");
2423 if (nr == 0 || nc == 0 || b.
cols () == 0)
2428 int typ = mattype.
type ();
2432 (*current_liboctave_error_handler) (
"incorrect matrix type");
2435 double ainvnorm = 0.;
2445 atmp += fabs (
data (i));
2454 retval.
xcidx (0) = 0;
2484 goto triangular_error;
2490 i <
cidx (kidx+1)-1; i++)
2493 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2505 if (ii + new_nnz > x_nz)
2514 if (cwork[rperm[i]] != 0.)
2516 retval.
xridx (ii) = i;
2517 retval.
xdata (ii++) = cwork[rperm[i]];
2519 retval.
xcidx (j+1) = ii;
2541 double tmp = work[k] /
data (
cidx (iidx+1)-1);
2544 i <
cidx (iidx+1)-1; i++)
2547 work[idx2] = work[idx2] - tmp *
data (i);
2554 atmp += fabs (work[i]);
2557 if (atmp > ainvnorm)
2560 rcond = 1. / ainvnorm / anorm;
2582 goto triangular_error;
2590 cwork[iidx] = cwork[iidx] - tmp *
data (i);
2602 if (ii + new_nnz > x_nz)
2613 retval.
xridx (ii) = i;
2614 retval.
xdata (ii++) = cwork[i];
2616 retval.
xcidx (j+1) = ii;
2636 double tmp = work[k] /
data (
cidx (k+1)-1);
2639 i <
cidx (k+1)-1; i++)
2642 work[iidx] = work[iidx] - tmp *
data (i);
2649 atmp += fabs (work[i]);
2652 if (atmp > ainvnorm)
2655 rcond = 1. / ainvnorm / anorm;
2664 sing_handler (rcond);
2668 octave::warn_singular_matrix (rcond);
2671 if (is_singular (rcond) || octave::math::isnan (rcond))
2677 sing_handler (rcond);
2681 octave::warn_singular_matrix (rcond);
2691 solve_singularity_handler sing_handler,
2692 bool calc_cond)
const
2701 if (nr != b.
rows ())
2702 (*current_liboctave_error_handler)
2703 (
"matrix dimension mismatch solution of linear equations");
2705 if (nr == 0 || nc == 0 || b.
cols () == 0)
2710 int typ = mattype.
type ();
2714 (*current_liboctave_error_handler) (
"incorrect matrix type");
2717 double ainvnorm = 0.;
2728 atmp += fabs (
data (i));
2736 retval.
resize (nc, b_nc);
2742 inv_perm[perm[i]] = i;
2750 work[inv_perm[i]] = b(i, j);
2761 if (inv_perm[
ridx (i)] < minr)
2763 minr = inv_perm[
ridx (i)];
2767 if (minr != k ||
data (mini) == 0.)
2770 goto triangular_error;
2773 double tmp = work[k] /
data (mini);
2781 work[iidx] = work[iidx] - tmp *
data (i);
2787 retval(i, j) = work[i];
2808 i <
cidx (k+1); i++)
2809 if (inv_perm[
ridx (i)] < minr)
2811 minr = inv_perm[
ridx (i)];
2815 double tmp = work[k] /
data (mini);
2818 i <
cidx (k+1); i++)
2824 work[iidx] = work[iidx] - tmp *
data (i);
2832 atmp += fabs (work[i]);
2835 if (atmp > ainvnorm)
2838 rcond = 1. / ainvnorm / anorm;
2844 retval.
resize (nc, b_nc, 0.);
2859 goto triangular_error;
2862 double tmp = work[k] /
data (
cidx (k));
2867 work[iidx] = work[iidx] - tmp *
data (i);
2873 retval.
xelem (i, j) = work[i];
2891 double tmp = work[k] /
data (
cidx (k));
2894 i <
cidx (k+1); i++)
2897 work[iidx] = work[iidx] - tmp *
data (i);
2904 atmp += fabs (work[i]);
2907 if (atmp > ainvnorm)
2910 rcond = 1. / ainvnorm / anorm;
2919 sing_handler (rcond);
2923 octave::warn_singular_matrix (rcond);
2926 if (is_singular (rcond) || octave::math::isnan (rcond))
2932 sing_handler (rcond);
2936 octave::warn_singular_matrix (rcond);
2946 solve_singularity_handler sing_handler,
2947 bool calc_cond)
const
2956 if (nr != b.
rows ())
2957 (*current_liboctave_error_handler)
2958 (
"matrix dimension mismatch solution of linear equations");
2960 if (nr == 0 || nc == 0 || b.
cols () == 0)
2965 int typ = mattype.
type ();
2969 (*current_liboctave_error_handler) (
"incorrect matrix type");
2972 double ainvnorm = 0.;
2982 atmp += fabs (
data (i));
2991 retval.
xcidx (0) = 0;
3002 inv_perm[perm[i]] = i;
3009 work[inv_perm[b.
ridx (i)]] = b.
data (i);
3020 if (inv_perm[
ridx (i)] < minr)
3022 minr = inv_perm[
ridx (i)];
3026 if (minr != k ||
data (mini) == 0.)
3029 goto triangular_error;
3032 double tmp = work[k] /
data (mini);
3040 work[iidx] = work[iidx] - tmp *
data (i);
3052 if (ii + new_nnz > x_nz)
3063 retval.
xridx (ii) = i;
3064 retval.
xdata (ii++) = work[i];
3066 retval.
xcidx (j+1) = ii;
3089 i <
cidx (k+1); i++)
3090 if (inv_perm[
ridx (i)] < minr)
3092 minr = inv_perm[
ridx (i)];
3096 double tmp = work[k] /
data (mini);
3099 i <
cidx (k+1); i++)
3105 work[iidx] = work[iidx] - tmp *
data (i);
3113 atmp += fabs (work[i]);
3116 if (atmp > ainvnorm)
3119 rcond = 1. / ainvnorm / anorm;
3140 goto triangular_error;
3143 double tmp = work[k] /
data (
cidx (k));
3148 work[iidx] = work[iidx] - tmp *
data (i);
3160 if (ii + new_nnz > x_nz)
3171 retval.
xridx (ii) = i;
3172 retval.
xdata (ii++) = work[i];
3174 retval.
xcidx (j+1) = ii;
3194 double tmp = work[k] /
data (
cidx (k));
3197 i <
cidx (k+1); i++)
3200 work[iidx] = work[iidx] - tmp *
data (i);
3207 atmp += fabs (work[i]);
3210 if (atmp > ainvnorm)
3213 rcond = 1. / ainvnorm / anorm;
3222 sing_handler (rcond);
3226 octave::warn_singular_matrix (rcond);
3229 if (is_singular (rcond) || octave::math::isnan (rcond))
3235 sing_handler (rcond);
3239 octave::warn_singular_matrix (rcond);
3249 solve_singularity_handler sing_handler,
3250 bool calc_cond)
const
3259 if (nr != b.
rows ())
3260 (*current_liboctave_error_handler)
3261 (
"matrix dimension mismatch solution of linear equations");
3263 if (nr == 0 || nc == 0 || b.
cols () == 0)
3268 int typ = mattype.
type ();
3272 (*current_liboctave_error_handler) (
"incorrect matrix type");
3275 double ainvnorm = 0.;
3286 atmp += fabs (
data (i));
3294 retval.
resize (nc, b_nc);
3300 inv_perm[perm[i]] = i;
3307 cwork[inv_perm[i]] = b(i, j);
3317 if (inv_perm[
ridx (i)] < minr)
3319 minr = inv_perm[
ridx (i)];
3323 if (minr != k ||
data (mini) == 0)
3326 goto triangular_error;
3337 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3343 retval(i, j) = cwork[i];
3365 i <
cidx (k+1); i++)
3366 if (inv_perm[
ridx (i)] < minr)
3368 minr = inv_perm[
ridx (i)];
3372 double tmp = work[k] /
data (mini);
3375 i <
cidx (k+1); i++)
3381 work[iidx] = work[iidx] - tmp *
data (i);
3389 atmp += fabs (work[i]);
3392 if (atmp > ainvnorm)
3395 rcond = 1. / ainvnorm / anorm;
3401 retval.
resize (nc, b_nc, 0.);
3417 goto triangular_error;
3425 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3431 retval.
xelem (i, j) = cwork[i];
3450 double tmp = work[k] /
data (
cidx (k));
3453 i <
cidx (k+1); i++)
3456 work[iidx] = work[iidx] - tmp *
data (i);
3463 atmp += fabs (work[i]);
3466 if (atmp > ainvnorm)
3469 rcond = 1. / ainvnorm / anorm;
3478 sing_handler (rcond);
3482 octave::warn_singular_matrix (rcond);
3485 if (is_singular (rcond) || octave::math::isnan (rcond))
3491 sing_handler (rcond);
3495 octave::warn_singular_matrix (rcond);
3505 solve_singularity_handler sing_handler,
3506 bool calc_cond)
const
3515 if (nr != b.
rows ())
3516 (*current_liboctave_error_handler)
3517 (
"matrix dimension mismatch solution of linear equations");
3519 if (nr == 0 || nc == 0 || b.
cols () == 0)
3524 int typ = mattype.
type ();
3528 (*current_liboctave_error_handler) (
"incorrect matrix type");
3531 double ainvnorm = 0.;
3541 atmp += fabs (
data (i));
3550 retval.
xcidx (0) = 0;
3561 inv_perm[perm[i]] = i;
3568 cwork[inv_perm[b.
ridx (i)]] = b.
data (i);
3578 if (inv_perm[
ridx (i)] < minr)
3580 minr = inv_perm[
ridx (i)];
3584 if (minr != k ||
data (mini) == 0.)
3587 goto triangular_error;
3598 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3610 if (ii + new_nnz > x_nz)
3621 retval.
xridx (ii) = i;
3622 retval.
xdata (ii++) = cwork[i];
3624 retval.
xcidx (j+1) = ii;
3648 i <
cidx (k+1); i++)
3649 if (inv_perm[
ridx (i)] < minr)
3651 minr = inv_perm[
ridx (i)];
3655 double tmp = work[k] /
data (mini);
3658 i <
cidx (k+1); i++)
3664 work[iidx] = work[iidx] - tmp *
data (i);
3672 atmp += fabs (work[i]);
3675 if (atmp > ainvnorm)
3678 rcond = 1. / ainvnorm / anorm;
3699 goto triangular_error;
3707 cwork[iidx] = cwork[iidx] - tmp *
data (i);
3719 if (ii + new_nnz > x_nz)
3730 retval.
xridx (ii) = i;
3731 retval.
xdata (ii++) = cwork[i];
3733 retval.
xcidx (j+1) = ii;
3754 double tmp = work[k] /
data (
cidx (k));
3757 i <
cidx (k+1); i++)
3760 work[iidx] = work[iidx] - tmp *
data (i);
3767 atmp += fabs (work[i]);
3770 if (atmp > ainvnorm)
3773 rcond = 1. / ainvnorm / anorm;
3782 sing_handler (rcond);
3786 octave::warn_singular_matrix (rcond);
3789 if (is_singular (rcond) || octave::math::isnan (rcond))
3795 sing_handler (rcond);
3799 octave::warn_singular_matrix (rcond);
3809 solve_singularity_handler sing_handler,
3810 bool calc_cond)
const
3818 if (nr != nc || nr != b.
rows ())
3819 (*current_liboctave_error_handler)
3820 (
"matrix dimension mismatch solution of linear equations");
3822 if (nr == 0 || b.
cols () == 0)
3825 (*current_liboctave_error_handler)
3826 (
"calculation of condition number not implemented");
3830 int typ = mattype.
type ();
3848 D[nc-1] =
data (ii);
3864 else if (
ridx (i) == j + 1)
3869 F77_INT tmp_nr = octave::to_f77_int (nr);
3875 double *result = retval.
rwdata ();
3879 F77_XFCN (dptsv, DPTSV, (tmp_nr, b_nc, D, DL, result,
3907 DL[j] =
data (ii++);
3908 DU[j] =
data (ii++);
3910 D[nc-1] =
data (ii);
3927 else if (
ridx (i) == j + 1)
3929 else if (
ridx (i) == j - 1)
3934 F77_INT tmp_nr = octave::to_f77_int (nr);
3940 double *result = retval.
rwdata ();
3944 F77_XFCN (dgtsv, DGTSV, (tmp_nr, b_nc, DL, D, DU, result,
3956 sing_handler (rcond);
3960 octave::warn_singular_matrix ();
3966 (*current_liboctave_error_handler) (
"incorrect matrix type");
3975 solve_singularity_handler sing_handler,
3976 bool calc_cond)
const
3984 if (nr != nc || nr != b.
rows ())
3985 (*current_liboctave_error_handler)
3986 (
"matrix dimension mismatch solution of linear equations");
3988 if (nr == 0 || b.
cols () == 0)
3991 (*current_liboctave_error_handler)
3992 (
"calculation of condition number not implemented");
3996 int typ = mattype.
type ();
4008 F77_INT *pipvt = ipvt.rwdata ();
4017 DL[j] =
data (ii++);
4018 DU[j] =
data (ii++);
4020 D[nc-1] =
data (ii);
4037 else if (
ridx (i) == j + 1)
4039 else if (
ridx (i) == j - 1)
4044 F77_INT tmp_nr = octave::to_f77_int (nr);
4048 F77_XFCN (dgttrf, DGTTRF, (tmp_nr, DL, D, DU, DU2, pipvt, tmp_err));
4057 sing_handler (rcond);
4061 octave::warn_singular_matrix ();
4070 retval.
xcidx (0) = 0;
4085 (F77_CONST_CHAR_ARG2 (&job, 1),
4086 tmp_nr, 1, DL, D, DU, DU2, pipvt,
4088 F77_CHAR_ARG_LEN (1)));
4099 if (ii + new_nnz > x_nz)
4110 retval.
xridx (ii) = i;
4111 retval.
xdata (ii++) = work[i];
4113 retval.
xcidx (j+1) = ii;
4120 (*current_liboctave_error_handler) (
"incorrect matrix type");
4129 solve_singularity_handler sing_handler,
4130 bool calc_cond)
const
4138 if (nr != nc || nr != b.
rows ())
4139 (*current_liboctave_error_handler)
4140 (
"matrix dimension mismatch solution of linear equations");
4142 if (nr == 0 || b.
cols () == 0)
4145 (*current_liboctave_error_handler)
4146 (
"calculation of condition number not implemented");
4150 int typ = mattype.
type ();
4168 D[nc-1] =
data (ii);
4184 else if (
ridx (i) == j + 1)
4189 F77_INT tmp_nr = octave::to_f77_int (nr);
4228 DL[j] =
data (ii++);
4229 DU[j] =
data (ii++);
4231 D[nc-1] =
data (ii);
4248 else if (
ridx (i) == j + 1)
4250 else if (
ridx (i) == j - 1)
4255 F77_INT tmp_nr = octave::to_f77_int (nr);
4282 sing_handler (rcond);
4286 octave::warn_singular_matrix ();
4290 (*current_liboctave_error_handler) (
"incorrect matrix type");
4299 solve_singularity_handler sing_handler,
4300 bool calc_cond)
const
4308 if (nr != nc || nr != b.
rows ())
4309 (*current_liboctave_error_handler)
4310 (
"matrix dimension mismatch solution of linear equations");
4312 if (nr == 0 || b.
cols () == 0)
4315 (*current_liboctave_error_handler)
4316 (
"calculation of condition number not implemented");
4320 int typ = mattype.
type ();
4332 F77_INT *pipvt = ipvt.rwdata ();
4341 DL[j] =
data (ii++);
4342 DU[j] =
data (ii++);
4344 D[nc-1] =
data (ii);
4361 else if (
ridx (i) == j + 1)
4363 else if (
ridx (i) == j - 1)
4368 F77_INT tmp_nr = octave::to_f77_int (nr);
4372 F77_XFCN (dgttrf, DGTTRF, (tmp_nr, DL, D, DU, DU2, pipvt, tmp_err));
4383 sing_handler (rcond);
4387 octave::warn_singular_matrix ();
4404 retval.
xcidx (0) = 0;
4408 for (
F77_INT i = 0; i < b_nr; i++)
4416 (F77_CONST_CHAR_ARG2 (&job, 1),
4417 tmp_nr, 1, DL, D, DU, DU2, pipvt,
4419 F77_CHAR_ARG_LEN (1)));
4426 (*current_liboctave_error_handler)
4427 (
"SparseMatrix::solve solve failed");
4434 (F77_CONST_CHAR_ARG2 (&job, 1),
4435 tmp_nr, 1, DL, D, DU, DU2, pipvt,
4437 F77_CHAR_ARG_LEN (1)));
4444 (*current_liboctave_error_handler)
4445 (
"SparseMatrix::solve solve failed");
4455 if (Bx[i] != 0. || Bz[i] != 0.)
4458 if (ii + new_nnz > x_nz)
4467 if (Bx[i] != 0. || Bz[i] != 0.)
4469 retval.
xridx (ii) = i;
4473 retval.
xcidx (j+1) = ii;
4480 (*current_liboctave_error_handler) (
"incorrect matrix type");
4489 solve_singularity_handler sing_handler,
4490 bool calc_cond)
const
4498 if (nr != nc || nr != b.
rows ())
4499 (*current_liboctave_error_handler)
4500 (
"matrix dimension mismatch solution of linear equations");
4502 if (nr == 0 || b.
cols () == 0)
4507 int typ = mattype.
type ();
4515 double *tmp_data = m_band.rwdata ();
4523 tmp_data[ii++] = 0.;
4531 m_band(ri - j, j) =
data (i);
4537 anorm = m_band.abs ().sum ().row (0).max ();
4539 F77_INT tmp_nr = octave::to_f77_int (nr);
4544 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4545 tmp_nr, n_lower, tmp_data, ldm, tmp_err
4546 F77_CHAR_ARG_LEN (1)));
4564 double *pz = z.rwdata ();
4569 (F77_CONST_CHAR_ARG2 (&job, 1),
4570 tmp_nr, n_lower, tmp_data, ldm,
4571 anorm, rcond, pz, piz, tmp_err
4572 F77_CHAR_ARG_LEN (1)));
4579 if (is_singular (rcond) || octave::math::isnan (rcond))
4585 sing_handler (rcond);
4589 octave::warn_singular_matrix (rcond);
4598 double *result = retval.
rwdata ();
4604 (F77_CONST_CHAR_ARG2 (&job, 1),
4605 tmp_nr, n_lower, b_nc, tmp_data,
4606 ldm, result, b_nr, tmp_err
4607 F77_CHAR_ARG_LEN (1)));
4614 (*current_liboctave_error_handler)
4615 (
"SparseMatrix::solve solve failed");
4627 F77_INT ldm = n_upper + 2 * n_lower + 1;
4630 double *tmp_data = m_band.rwdata ();
4636 for (
F77_INT j = 0; j < ldm; j++)
4638 tmp_data[ii++] = 0.;
4643 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
4653 atmp += fabs (
data (i));
4659 F77_INT tmp_nr = octave::to_f77_int (nr);
4662 F77_INT *pipvt = ipvt.rwdata ();
4666 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4667 tmp_data, ldm, pipvt, tmp_err));
4680 sing_handler (rcond);
4684 octave::warn_singular_matrix ();
4692 double *pz = z.rwdata ();
4696 F77_INT tmp_nc = octave::to_f77_int (nc);
4699 (F77_CONST_CHAR_ARG2 (&job, 1),
4700 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4701 anorm, rcond, pz, piz, tmp_err
4702 F77_CHAR_ARG_LEN (1)));
4709 if (is_singular (rcond) || octave::math::isnan (rcond))
4715 sing_handler (rcond);
4719 octave::warn_singular_matrix (rcond);
4728 double *result = retval.
rwdata ();
4735 (F77_CONST_CHAR_ARG2 (&job, 1),
4736 tmp_nr, n_lower, n_upper, b_nc, tmp_data,
4737 ldm, pipvt, result, b_nr, tmp_err
4738 F77_CHAR_ARG_LEN (1)));
4745 (*current_liboctave_error_handler) (
"incorrect matrix type");
4754 solve_singularity_handler sing_handler,
4755 bool calc_cond)
const
4763 if (nr != nc || nr != b.
rows ())
4764 (*current_liboctave_error_handler)
4765 (
"matrix dimension mismatch solution of linear equations");
4767 if (nr == 0 || b.
cols () == 0)
4772 int typ = mattype.
type ();
4778 F77_INT ldm = octave::to_f77_int (n_lower + 1);
4781 double *tmp_data = m_band.rwdata ();
4787 for (
F77_INT j = 0; j < ldm; j++)
4789 tmp_data[ii++] = 0.;
4797 m_band(ri - j, j) =
data (i);
4803 anorm = m_band.abs ().sum ().row (0).max ();
4805 F77_INT tmp_nr = octave::to_f77_int (nr);
4810 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4811 tmp_nr, n_lower, tmp_data, ldm, tmp_err
4812 F77_CHAR_ARG_LEN (1)));
4828 double *pz = z.rwdata ();
4833 (F77_CONST_CHAR_ARG2 (&job, 1),
4834 tmp_nr, n_lower, tmp_data, ldm,
4835 anorm, rcond, pz, piz, tmp_err
4836 F77_CHAR_ARG_LEN (1)));
4843 if (is_singular (rcond) || octave::math::isnan (rcond))
4849 sing_handler (rcond);
4853 octave::warn_singular_matrix (rcond);
4871 retval.
xcidx (0) = 0;
4874 for (
F77_INT i = 0; i < b_nr; i++)
4875 Bx[i] = b.
elem (i, j);
4878 (F77_CONST_CHAR_ARG2 (&job, 1),
4879 tmp_nr, n_lower, 1, tmp_data,
4880 ldm, Bx, b_nr, tmp_err
4881 F77_CHAR_ARG_LEN (1)));
4888 (*current_liboctave_error_handler)
4889 (
"SparseMatrix::solve solve failed");
4894 for (
F77_INT i = 0; i < b_nr; i++)
4903 sz = (
static_cast<double> (b_nc) - j) / b_nc
4905 sz = x_nz + (sz > 100 ? sz : 100);
4909 retval.
xdata (ii) = tmp;
4910 retval.
xridx (ii++) = i;
4913 retval.
xcidx (j+1) = ii;
4926 F77_INT ldm = octave::to_f77_int (n_upper + 2 * n_lower + 1);
4929 double *tmp_data = m_band.rwdata ();
4937 tmp_data[ii++] = 0.;
4942 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
4953 atmp += fabs (
data (i));
4959 F77_INT tmp_nr = octave::to_f77_int (nr);
4962 F77_INT *pipvt = ipvt.rwdata ();
4966 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4967 tmp_data, ldm, pipvt, tmp_err));
4978 sing_handler (rcond);
4982 octave::warn_singular_matrix ();
4990 double *pz = z.rwdata ();
4994 F77_INT tmp_nc = octave::to_f77_int (nc);
4997 (F77_CONST_CHAR_ARG2 (&job, 1),
4998 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4999 anorm, rcond, pz, piz, tmp_err
5000 F77_CHAR_ARG_LEN (1)));
5007 if (is_singular (rcond) || octave::math::isnan (rcond))
5013 sing_handler (rcond);
5017 octave::warn_singular_matrix (rcond);
5029 retval.
xcidx (0) = 0;
5039 i < b.
cidx (j+1); i++)
5045 (F77_CONST_CHAR_ARG2 (&job, 1),
5046 tmp_nr, n_lower, n_upper, 1, tmp_data,
5047 ldm, pipvt, work, b_nr, tmp_err
5048 F77_CHAR_ARG_LEN (1)));
5059 if (ii + new_nnz > x_nz)
5070 retval.
xridx (ii) = i;
5071 retval.
xdata (ii++) = work[i];
5073 retval.
xcidx (j+1) = ii;
5081 (*current_liboctave_error_handler) (
"incorrect matrix type");
5090 solve_singularity_handler sing_handler,
5091 bool calc_cond)
const
5099 if (nr != nc || nr != b.
rows ())
5100 (*current_liboctave_error_handler)
5101 (
"matrix dimension mismatch solution of linear equations");
5103 if (nr == 0 || b.
cols () == 0)
5108 int typ = mattype.
type ();
5117 double *tmp_data = m_band.rwdata ();
5123 for (
F77_INT j = 0; j < ldm; j++)
5125 tmp_data[ii++] = 0.;
5133 m_band(ri - j, j) =
data (i);
5139 anorm = m_band.abs ().sum ().row (0).max ();
5141 F77_INT tmp_nr = octave::to_f77_int (nr);
5146 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5147 tmp_nr, n_lower, tmp_data, ldm, tmp_err
5148 F77_CHAR_ARG_LEN (1)));
5166 double *pz = z.rwdata ();
5171 (F77_CONST_CHAR_ARG2 (&job, 1),
5172 tmp_nr, n_lower, tmp_data, ldm,
5173 anorm, rcond, pz, piz, tmp_err
5174 F77_CHAR_ARG_LEN (1)));
5181 if (is_singular (rcond) || octave::math::isnan (rcond))
5187 sing_handler (rcond);
5191 octave::warn_singular_matrix (rcond);
5205 retval.
resize (b_nr, b_nc);
5209 for (
F77_INT i = 0; i < b_nr; i++)
5217 (F77_CONST_CHAR_ARG2 (&job, 1),
5218 tmp_nr, n_lower, 1, tmp_data,
5219 ldm, Bx, b_nr, tmp_err
5220 F77_CHAR_ARG_LEN (1)));
5227 (*current_liboctave_error_handler)
5228 (
"SparseMatrix::solve solve failed");
5234 (F77_CONST_CHAR_ARG2 (&job, 1),
5235 tmp_nr, n_lower, 1, tmp_data,
5236 ldm, Bz, b_nr, tmp_err
5237 F77_CHAR_ARG_LEN (1)));
5244 (*current_liboctave_error_handler)
5245 (
"SparseMatrix::solve solve failed");
5251 retval(i, j) =
Complex (Bx[i], Bz[i]);
5262 F77_INT ldm = n_upper + 2 * n_lower + 1;
5265 double *tmp_data = m_band.rwdata ();
5271 for (
F77_INT j = 0; j < ldm; j++)
5273 tmp_data[ii++] = 0.;
5278 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
5289 atmp += fabs (
data (i));
5295 F77_INT tmp_nr = octave::to_f77_int (nr);
5298 F77_INT *pipvt = ipvt.rwdata ();
5302 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5303 tmp_data, ldm, pipvt, tmp_err));
5314 sing_handler (rcond);
5318 octave::warn_singular_matrix ();
5326 double *pz = z.rwdata ();
5330 F77_INT tmp_nc = octave::to_f77_int (nc);
5333 (F77_CONST_CHAR_ARG2 (&job, 1),
5334 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5335 anorm, rcond, pz, piz, tmp_err
5336 F77_CHAR_ARG_LEN (1)));
5343 if (is_singular (rcond) || octave::math::isnan (rcond))
5349 sing_handler (rcond);
5353 octave::warn_singular_matrix (rcond);
5364 retval.
resize (nr, b_nc);
5379 (F77_CONST_CHAR_ARG2 (&job, 1),
5380 tmp_nr, n_lower, n_upper, 1, tmp_data,
5381 ldm, pipvt, Bx, b_nr, tmp_err
5382 F77_CHAR_ARG_LEN (1)));
5387 (F77_CONST_CHAR_ARG2 (&job, 1),
5388 tmp_nr, n_lower, n_upper, 1, tmp_data,
5389 ldm, pipvt, Bz, b_nr, tmp_err
5390 F77_CHAR_ARG_LEN (1)));
5395 retval(i, j) =
Complex (Bx[i], Bz[i]);
5401 (*current_liboctave_error_handler) (
"incorrect matrix type");
5410 solve_singularity_handler sing_handler,
5411 bool calc_cond)
const
5419 if (nr != nc || nr != b.
rows ())
5420 (*current_liboctave_error_handler)
5421 (
"matrix dimension mismatch solution of linear equations");
5423 if (nr == 0 || b.
cols () == 0)
5428 int typ = mattype.
type ();
5437 double *tmp_data = m_band.rwdata ();
5443 for (
F77_INT j = 0; j < ldm; j++)
5445 tmp_data[ii++] = 0.;
5453 m_band(ri - j, j) =
data (i);
5459 anorm = m_band.abs ().sum ().row (0).max ();
5461 F77_INT tmp_nr = octave::to_f77_int (nr);
5466 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5467 tmp_nr, n_lower, tmp_data, ldm, tmp_err
5468 F77_CHAR_ARG_LEN (1)));
5487 double *pz = z.rwdata ();
5492 (F77_CONST_CHAR_ARG2 (&job, 1),
5493 tmp_nr, n_lower, tmp_data, ldm,
5494 anorm, rcond, pz, piz, tmp_err
5495 F77_CHAR_ARG_LEN (1)));
5502 if (is_singular (rcond) || octave::math::isnan (rcond))
5508 sing_handler (rcond);
5512 octave::warn_singular_matrix (rcond);
5531 retval.
xcidx (0) = 0;
5535 for (
F77_INT i = 0; i < b_nr; i++)
5543 (F77_CONST_CHAR_ARG2 (&job, 1),
5544 tmp_nr, n_lower, 1, tmp_data,
5545 ldm, Bx, b_nr, tmp_err
5546 F77_CHAR_ARG_LEN (1)));
5553 (*current_liboctave_error_handler)
5554 (
"SparseMatrix::solve solve failed");
5560 (F77_CONST_CHAR_ARG2 (&job, 1),
5561 tmp_nr, n_lower, 1, tmp_data,
5562 ldm, Bz, b_nr, tmp_err
5563 F77_CHAR_ARG_LEN (1)));
5570 (*current_liboctave_error_handler)
5571 (
"SparseMatrix::solve solve failed");
5581 if (Bx[i] != 0. || Bz[i] != 0.)
5584 if (ii + new_nnz > x_nz)
5593 if (Bx[i] != 0. || Bz[i] != 0.)
5595 retval.
xridx (ii) = i;
5599 retval.
xcidx (j+1) = ii;
5612 F77_INT ldm = n_upper + 2 * n_lower + 1;
5615 double *tmp_data = m_band.rwdata ();
5621 for (
F77_INT j = 0; j < ldm; j++)
5623 tmp_data[ii++] = 0.;
5628 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
5639 atmp += fabs (
data (i));
5645 F77_INT tmp_nr = octave::to_f77_int (nr);
5648 F77_INT *pipvt = ipvt.rwdata ();
5652 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5653 tmp_data, ldm, pipvt, tmp_err));
5664 sing_handler (rcond);
5668 octave::warn_singular_matrix ();
5676 double *pz = z.rwdata ();
5680 F77_INT tmp_nc = octave::to_f77_int (nc);
5683 (F77_CONST_CHAR_ARG2 (&job, 1),
5684 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5685 anorm, rcond, pz, piz, tmp_err
5686 F77_CHAR_ARG_LEN (1)));
5693 if (is_singular (rcond) || octave::math::isnan (rcond))
5699 sing_handler (rcond);
5703 octave::warn_singular_matrix (rcond);
5716 retval.
xcidx (0) = 0;
5730 i < b.
cidx (j+1); i++)
5733 Bx[b.
ridx (i)] = c.real ();
5734 Bz[b.
ridx (i)] = c.imag ();
5738 (F77_CONST_CHAR_ARG2 (&job, 1),
5739 tmp_nr, n_lower, n_upper, 1, tmp_data,
5740 ldm, pipvt, Bx, b_nr, tmp_err
5741 F77_CHAR_ARG_LEN (1)));
5746 (F77_CONST_CHAR_ARG2 (&job, 1),
5747 tmp_nr, n_lower, n_upper, 1, tmp_data,
5748 ldm, pipvt, Bz, b_nr, tmp_err
5749 F77_CHAR_ARG_LEN (1)));
5757 if (Bx[i] != 0. || Bz[i] != 0.)
5760 if (ii + new_nnz > x_nz)
5769 if (Bx[i] != 0. || Bz[i] != 0.)
5771 retval.
xridx (ii) = i;
5774 retval.
xcidx (j+1) = ii;
5782 (*current_liboctave_error_handler) (
"incorrect matrix type");
5790 Matrix& Info, solve_singularity_handler sing_handler,
5791 bool calc_cond)
const
5794 void *Numeric =
nullptr;
5797#if defined (HAVE_UMFPACK)
5800 Control =
Matrix (UMFPACK_CONTROL, 1);
5801 double *control = Control.
rwdata ();
5804 double tmp = octave::sparse_params::get_key (
"spumoni");
5805 if (! octave::math::isnan (tmp))
5806 Control (UMFPACK_PRL) = tmp;
5807 tmp = octave::sparse_params::get_key (
"piv_tol");
5808 if (! octave::math::isnan (tmp))
5810 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
5811 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
5815 tmp = octave::sparse_params::get_key (
"autoamd");
5816 if (! octave::math::isnan (tmp))
5817 Control (UMFPACK_FIXQ) = tmp;
5823 const double *Ax =
data ();
5828 octave::to_suitesparse_intptr (Ap),
5829 octave::to_suitesparse_intptr (Ai),
5833 Info =
Matrix (1, UMFPACK_INFO);
5834 double *info = Info.
rwdata ();
5836 octave::to_suitesparse_intptr (Ap),
5837 octave::to_suitesparse_intptr (Ai),
5838 Ax,
nullptr, &Symbolic, control, info);
5848 (*current_liboctave_error_handler)
5849 (
"SparseMatrix::solve symbolic factorization failed");
5856 status =
UMFPACK_DNAME (numeric) (octave::to_suitesparse_intptr (Ap),
5857 octave::to_suitesparse_intptr (Ai),
5858 Ax, 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 (
"SparseMatrix::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 "
5908 "when liboctave was built");
5918 solve_singularity_handler sing_handler,
5919 bool calc_cond)
const
5927 if (nr != nc || nr != b.
rows ())
5928 (*current_liboctave_error_handler)
5929 (
"matrix dimension mismatch solution of linear equations");
5931 if (nr == 0 || b.
cols () == 0)
5936 int typ = mattype.
type ();
5941#if defined (HAVE_CHOLMOD)
5942 cholmod_common Common;
5943 cholmod_common *cm = &Common;
5947 cm->prefer_zomplex =
false;
5949 double spu = octave::sparse_params::get_key (
"spumoni");
5953 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
5957 cm->print =
static_cast<int> (spu) + 2;
5958 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5962 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5963 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5965 cm->final_ll =
true;
5967 cholmod_sparse Astore;
5968 cholmod_sparse *
A = &Astore;
5978#if defined (OCTAVE_ENABLE_64)
5979 A->itype = CHOLMOD_LONG;
5981 A->itype = CHOLMOD_INT;
5983 A->dtype = CHOLMOD_DOUBLE;
5985 A->xtype = CHOLMOD_REAL;
5989 cholmod_dense Bstore;
5990 cholmod_dense *
B = &Bstore;
5991 B->nrow = b.
rows ();
5992 B->ncol = b.
cols ();
5994 B->nzmax =
B->nrow *
B->ncol;
5995 B->dtype = CHOLMOD_DOUBLE;
5996 B->xtype = CHOLMOD_REAL;
5998 B->x =
const_cast<double *
> (b.
data ());
6015 if (is_singular (rcond) || octave::math::isnan (rcond))
6021 sing_handler (rcond);
6025 octave::warn_singular_matrix (rcond);
6037 retval.
xelem (i, j) =
static_cast<double *
> (X->x)[jr + i];
6043 static char blank_name[] =
" ";
6047 (*current_liboctave_warning_with_id_handler)
6048 (
"Octave:missing-dependency",
6049 "support for CHOLMOD was unavailable or disabled "
6050 "when liboctave was built");
6059#if defined (HAVE_UMFPACK)
6062 = factorize (err, rcond, Control, Info, sing_handler, calc_cond);
6067 Control (UMFPACK_IRSTEP) = 1;
6068 const double *Bx = b.
data ();
6070 double *result = retval.
rwdata ();
6074 double *control = Control.
rwdata ();
6075 double *info = Info.
rwdata ();
6078 const double *Ax =
data ();
6083 octave::to_suitesparse_intptr (Ap),
6084 octave::to_suitesparse_intptr (Ai),
6085 Ax, &result[iidx], &Bx[iidx],
6086 Numeric, control, info);
6092 (*current_liboctave_error_handler)
6093 (
"SparseMatrix::solve solve failed");
6108 octave_unused_parameter (rcond);
6109 octave_unused_parameter (sing_handler);
6110 octave_unused_parameter (calc_cond);
6112 (*current_liboctave_error_handler)
6113 (
"support for UMFPACK was unavailable or disabled "
6114 "when liboctave was built");
6118 (*current_liboctave_error_handler) (
"incorrect matrix type");
6127 solve_singularity_handler sing_handler,
6128 bool calc_cond)
const
6136 if (nr != nc || nr != b.
rows ())
6137 (*current_liboctave_error_handler)
6138 (
"matrix dimension mismatch solution of linear equations");
6140 if (nr == 0 || b.
cols () == 0)
6145 int typ = mattype.
type ();
6150#if defined (HAVE_CHOLMOD)
6151 cholmod_common Common;
6152 cholmod_common *cm = &Common;
6156 cm->prefer_zomplex =
false;
6158 double spu = octave::sparse_params::get_key (
"spumoni");
6162 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6166 cm->print =
static_cast<int> (spu) + 2;
6167 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6171 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6172 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6174 cm->final_ll =
true;
6176 cholmod_sparse Astore;
6177 cholmod_sparse *
A = &Astore;
6187#if defined (OCTAVE_ENABLE_64)
6188 A->itype = CHOLMOD_LONG;
6190 A->itype = CHOLMOD_INT;
6192 A->dtype = CHOLMOD_DOUBLE;
6194 A->xtype = CHOLMOD_REAL;
6198 cholmod_sparse Bstore;
6199 cholmod_sparse *
B = &Bstore;
6200 B->nrow = b.
rows ();
6201 B->ncol = b.
cols ();
6204 B->nzmax = b.
nnz ();
6208#if defined (OCTAVE_ENABLE_64)
6209 B->itype = CHOLMOD_LONG;
6211 B->itype = CHOLMOD_INT;
6213 B->dtype = CHOLMOD_DOUBLE;
6215 B->xtype = CHOLMOD_REAL;
6234 if (is_singular (rcond) || octave::math::isnan (rcond))
6240 sing_handler (rcond);
6244 octave::warn_singular_matrix (rcond);
6249 cholmod_sparse *X =
CHOLMOD_NAME(spsolve) (CHOLMOD_A, L,
B, cm);
6252 (octave::from_size_t (X->nrow),
6253 octave::from_size_t (X->ncol),
6257 j <= static_cast<octave_idx_type> (X->ncol); j++)
6264 retval.
xdata (j) =
static_cast<double *
> (X->x)[j];
6270 static char blank_name[] =
" ";
6274 (*current_liboctave_warning_with_id_handler)
6275 (
"Octave:missing-dependency",
6276 "support for CHOLMOD was unavailable or disabled "
6277 "when liboctave was built");
6286#if defined (HAVE_UMFPACK)
6288 void *Numeric = factorize (err, rcond, Control, Info,
6289 sing_handler, calc_cond);
6294 Control (UMFPACK_IRSTEP) = 1;
6298 double *control = Control.
rwdata ();
6299 double *info = Info.
rwdata ();
6302 const double *Ax =
data ();
6313 retval.
xcidx (0) = 0;
6318 Bx[i] = b.
elem (i, j);
6321 octave::to_suitesparse_intptr (Ap),
6322 octave::to_suitesparse_intptr (Ai),
6323 Ax, Xx, Bx, Numeric,
6330 (*current_liboctave_error_handler)
6331 (
"SparseMatrix::solve solve failed");
6346 sz = (
static_cast<double> (b_nc) - j) / b_nc
6348 sz = x_nz + (sz > 100 ? sz : 100);
6352 retval.
xdata (ii) = tmp;
6353 retval.
xridx (ii++) = i;
6356 retval.
xcidx (j+1) = ii;
6369 octave_unused_parameter (rcond);
6370 octave_unused_parameter (sing_handler);
6371 octave_unused_parameter (calc_cond);
6373 (*current_liboctave_error_handler)
6374 (
"support for UMFPACK was unavailable or disabled "
6375 "when liboctave was built");
6379 (*current_liboctave_error_handler) (
"incorrect matrix type");
6388 solve_singularity_handler sing_handler,
6389 bool calc_cond)
const
6397 if (nr != nc || nr != b.
rows ())
6398 (*current_liboctave_error_handler)
6399 (
"matrix dimension mismatch solution of linear equations");
6401 if (nr == 0 || b.
cols () == 0)
6406 int typ = mattype.
type ();
6411#if defined (HAVE_CHOLMOD)
6412 cholmod_common Common;
6413 cholmod_common *cm = &Common;
6417 cm->prefer_zomplex =
false;
6419 double spu = octave::sparse_params::get_key (
"spumoni");
6423 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6427 cm->print =
static_cast<int> (spu) + 2;
6428 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6432 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6433 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6435 cm->final_ll =
true;
6437 cholmod_sparse Astore;
6438 cholmod_sparse *
A = &Astore;
6448#if defined (OCTAVE_ENABLE_64)
6449 A->itype = CHOLMOD_LONG;
6451 A->itype = CHOLMOD_INT;
6453 A->dtype = CHOLMOD_DOUBLE;
6455 A->xtype = CHOLMOD_REAL;
6459 cholmod_dense Bstore;
6460 cholmod_dense *
B = &Bstore;
6461 B->nrow = b.
rows ();
6462 B->ncol = b.
cols ();
6464 B->nzmax =
B->nrow *
B->ncol;
6465 B->dtype = CHOLMOD_DOUBLE;
6466 B->xtype = CHOLMOD_COMPLEX;
6485 if (is_singular (rcond) || octave::math::isnan (rcond))
6491 sing_handler (rcond);
6495 octave::warn_singular_matrix (rcond);
6507 retval.
xelem (i, j) =
static_cast<Complex *
> (X->x)[jr + i];
6513 static char blank_name[] =
" ";
6517 (*current_liboctave_warning_with_id_handler)
6518 (
"Octave:missing-dependency",
6519 "support for CHOLMOD was unavailable or disabled "
6520 "when liboctave was built");
6529#if defined (HAVE_UMFPACK)
6531 void *Numeric = factorize (err, rcond, Control, Info,
6532 sing_handler, calc_cond);
6537 Control (UMFPACK_IRSTEP) = 1;
6541 double *control = Control.
rwdata ();
6542 double *info = Info.
rwdata ();
6545 const double *Ax =
data ();
6550 retval.
resize (b_nr, b_nc);
6565 octave::to_suitesparse_intptr (Ap),
6566 octave::to_suitesparse_intptr (Ai),
6567 Ax, Xx, Bx, Numeric,
6570 octave::to_suitesparse_intptr (Ap),
6571 octave::to_suitesparse_intptr (Ai),
6573 Numeric, control, info);
6575 if (status < 0 || status2 < 0)
6580 (*current_liboctave_error_handler)
6581 (
"SparseMatrix::solve solve failed");
6588 retval(i, j) =
Complex (Xx[i], Xz[i]);
6599 octave_unused_parameter (rcond);
6600 octave_unused_parameter (sing_handler);
6601 octave_unused_parameter (calc_cond);
6603 (*current_liboctave_error_handler)
6604 (
"support for UMFPACK was unavailable or disabled "
6605 "when liboctave was built");
6609 (*current_liboctave_error_handler) (
"incorrect matrix type");
6618 solve_singularity_handler sing_handler,
6619 bool calc_cond)
const
6627 if (nr != nc || nr != b.
rows ())
6628 (*current_liboctave_error_handler)
6629 (
"matrix dimension mismatch solution of linear equations");
6631 if (nr == 0 || b.
cols () == 0)
6636 int typ = mattype.
type ();
6641#if defined (HAVE_CHOLMOD)
6642 cholmod_common Common;
6643 cholmod_common *cm = &Common;
6647 cm->prefer_zomplex =
false;
6649 double spu = octave::sparse_params::get_key (
"spumoni");
6653 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6657 cm->print =
static_cast<int> (spu) + 2;
6658 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6662 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6663 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6665 cm->final_ll =
true;
6667 cholmod_sparse Astore;
6668 cholmod_sparse *
A = &Astore;
6678#if defined (OCTAVE_ENABLE_64)
6679 A->itype = CHOLMOD_LONG;
6681 A->itype = CHOLMOD_INT;
6683 A->dtype = CHOLMOD_DOUBLE;
6685 A->xtype = CHOLMOD_REAL;
6689 cholmod_sparse Bstore;
6690 cholmod_sparse *
B = &Bstore;
6691 B->nrow = b.
rows ();
6692 B->ncol = b.
cols ();
6695 B->nzmax = b.
nnz ();
6699#if defined (OCTAVE_ENABLE_64)
6700 B->itype = CHOLMOD_LONG;
6702 B->itype = CHOLMOD_INT;
6704 B->dtype = CHOLMOD_DOUBLE;
6706 B->xtype = CHOLMOD_COMPLEX;
6725 if (is_singular (rcond) || octave::math::isnan (rcond))
6731 sing_handler (rcond);
6735 octave::warn_singular_matrix (rcond);
6740 cholmod_sparse *X =
CHOLMOD_NAME(spsolve) (CHOLMOD_A, L,
B, cm);
6743 (octave::from_size_t (X->nrow),
6744 octave::from_size_t (X->ncol),
6748 j <= static_cast<octave_idx_type> (X->ncol); j++)
6761 static char blank_name[] =
" ";
6765 (*current_liboctave_warning_with_id_handler)
6766 (
"Octave:missing-dependency",
6767 "support for CHOLMOD was unavailable or disabled "
6768 "when liboctave was built");
6777#if defined (HAVE_UMFPACK)
6779 void *Numeric = factorize (err, rcond, Control, Info,
6780 sing_handler, calc_cond);
6785 Control (UMFPACK_IRSTEP) = 1;
6789 double *control = Control.
rwdata ();
6790 double *info = Info.
rwdata ();
6793 const double *Ax =
data ();
6807 retval.
xcidx (0) = 0;
6818 octave::to_suitesparse_intptr (Ap),
6819 octave::to_suitesparse_intptr (Ai),
6820 Ax, Xx, Bx, Numeric,
6823 octave::to_suitesparse_intptr (Ap),
6824 octave::to_suitesparse_intptr (Ai),
6826 Numeric, control, info);
6828 if (status < 0 || status2 < 0)
6833 (*current_liboctave_error_handler)
6834 (
"SparseMatrix::solve solve failed");
6849 sz = (
static_cast<double> (b_nc) - j) / b_nc
6851 sz = x_nz + (sz > 100 ? sz : 100);
6855 retval.
xdata (ii) = tmp;
6856 retval.
xridx (ii++) = i;
6859 retval.
xcidx (j+1) = ii;
6871 octave_unused_parameter (rcond);
6872 octave_unused_parameter (sing_handler);
6873 octave_unused_parameter (calc_cond);
6875 (*current_liboctave_error_handler)
6876 (
"support for UMFPACK was unavailable or disabled "
6877 "when liboctave was built");
6881 (*current_liboctave_error_handler) (
"incorrect matrix type");
6892 return solve (mattype, b, info, rcond,
nullptr);
6900 return solve (mattype, b, info, rcond,
nullptr);
6907 return solve (mattype, b, info, rcond,
nullptr);
6912 double& rcond, solve_singularity_handler sing_handler,
6913 bool singular_fallback)
const
6916 int typ = mattype.
type (
false);
6919 typ = mattype.
type (*
this);
6923 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
6925 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
6927 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
6929 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
6932 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
6934 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
6936 (*current_liboctave_error_handler) (
"unknown matrix type");
6953 return solve (mattype, b, info, rcond,
nullptr);
6961 return solve (mattype, b, info, rcond,
nullptr);
6968 return solve (mattype, b, info, rcond,
nullptr);
6974 solve_singularity_handler sing_handler,
6975 bool singular_fallback)
const
6978 int typ = mattype.
type (
false);
6981 typ = mattype.
type (*
this);
6984 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
6986 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
6988 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
6990 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
6993 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
6995 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
6997 (*current_liboctave_error_handler) (
"unknown matrix type");
7014 return solve (mattype, b, info, rcond,
nullptr);
7022 return solve (mattype, b, info, rcond,
nullptr);
7029 return solve (mattype, b, info, rcond,
nullptr);
7035 solve_singularity_handler sing_handler,
7036 bool singular_fallback)
const
7039 int typ = mattype.
type (
false);
7042 typ = mattype.
type (*
this);
7045 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
7047 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
7049 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
7051 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
7054 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
7056 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
7058 (*current_liboctave_error_handler) (
"unknown matrix type");
7075 return solve (mattype, b, info, rcond,
nullptr);
7083 return solve (mattype, b, info, rcond,
nullptr);
7090 return solve (mattype, b, info, rcond,
nullptr);
7096 solve_singularity_handler sing_handler,
7097 bool singular_fallback)
const
7100 int typ = mattype.
type (
false);
7103 typ = mattype.
type (*
this);
7106 retval = dsolve (mattype, b, err, rcond, sing_handler,
false);
7108 retval = utsolve (mattype, b, err, rcond, sing_handler,
false);
7110 retval = ltsolve (mattype, b, err, rcond, sing_handler,
false);
7112 retval = bsolve (mattype, b, err, rcond, sing_handler,
false);
7115 retval = trisolve (mattype, b, err, rcond, sing_handler,
false);
7117 retval = fsolve (mattype, b, err, rcond, sing_handler,
true);
7119 (*current_liboctave_error_handler) (
"unknown matrix type");
7135 return solve (mattype, b, info, rcond);
7143 return solve (mattype, b, info, rcond);
7150 return solve (mattype, b, info, rcond,
nullptr);
7156 solve_singularity_handler sing_handler)
const
7159 return solve (mattype, tmp, info, rcond,
7168 return solve (mattype, b, info, rcond,
nullptr);
7176 return solve (mattype, b, info, rcond,
nullptr);
7182 double& rcond)
const
7184 return solve (mattype, b, info, rcond,
nullptr);
7190 solve_singularity_handler sing_handler)
const
7193 return solve (mattype, tmp, info, rcond,
7202 return solve (b, info, rcond,
nullptr);
7209 return solve (b, info, rcond,
nullptr);
7214 double& rcond)
const
7216 return solve (b, info, rcond,
nullptr);
7221 solve_singularity_handler sing_handler)
const
7224 return solve (mattype, b, err, rcond, sing_handler);
7232 return solve (b, info, rcond,
nullptr);
7240 return solve (b, info, rcond,
nullptr);
7247 return solve (b, info, rcond,
nullptr);
7252 solve_singularity_handler sing_handler)
const
7255 return solve (mattype, b, err, rcond, sing_handler);
7262 return solve (b, info, rcond,
nullptr);
7267 double& rcond)
const
7269 return solve (b, info, rcond,
nullptr);
7275 solve_singularity_handler sing_handler)
const
7278 return solve (mattype, b, err, rcond, sing_handler);
7286 return solve (b, info, rcond,
nullptr);
7293 return solve (b, info, rcond,
nullptr);
7298 double& rcond)
const
7300 return solve (b, info, rcond,
nullptr);
7306 solve_singularity_handler sing_handler)
const
7309 return solve (mattype, b, err, rcond, sing_handler);
7316 return solve (b, info, rcond);
7323 return solve (b, info, rcond);
7328 double& rcond)
const
7330 return solve (b, info, rcond,
nullptr);
7336 solve_singularity_handler sing_handler)
const
7339 return solve (tmp, info, rcond,
7348 return solve (b, info, rcond,
nullptr);
7355 return solve (b, info, rcond,
nullptr);
7360 double& rcond)
const
7362 return solve (b, info, rcond,
nullptr);
7368 solve_singularity_handler sing_handler)
const
7371 return solve (tmp, info, rcond,
7385 if (octave::math::signbit (
data (i)))
7405 double val =
data (i);
7406 if (octave::math::isnan (val))
7420 double val =
data (i);
7421 if (octave::math::isinf (val) || octave::math::isnan (val))
7435 double val =
data (i);
7436 if (val != 0.0 && val != 1.0)
7462 double val =
data (i);
7463 if (octave::math::isnan (val) || octave::math::is_integer (val))
7488 double val =
data (i);
7496 if (! octave::math::is_integer (val))
7506 return test_any (octave::too_large_for_float);
7513 octave::err_nan_to_logical_conversion ();
7529 if (jj <
cidx (i+1) &&
ridx (jj) == j)
7571 if (! octave::math::isnan (data (j))) \
7592 else if ((
rows () == 1 && dim == -1) || dim == 1)
7597 double d = data (i); \
7598 if (nanflag && octave::math::isnan (d)) \
7599 tmp[ridx (i)] *= 1.0; \
7604 double d = data (i); \
7605 if (nanflag && octave::math::isnan (d)) \
7611 (
cidx (j+1) -
cidx (j) < nr ? 0.0 : 1.0), 1.0);
7622 double d = data (i); \
7623 if (nanflag && octave::math::isnan (d)) \
7624 tmp[ridx (i)] += 0.0; \
7629 double d = data (i); \
7630 if (nanflag && octave::math::isnan (d)) \
7655#define EXPR r.data (nel++) = data (i) * data (i);
7658 double d = data (i); \
7659 if (nanflag && octave::math::isnan (d)) \
7660 tmp[ridx (i)] += 0.0; \
7662 tmp[ridx (i)] += d * d
7665 double d = data (i); \
7666 if (nanflag && octave::math::isnan (d)) \
7688 retval.
data (i) = fabs (retval.
data (i));
7717 os << a.
ridx (i) + 1 <<
' ' << j + 1 <<
' ';
7718 octave::write_value<double> (os, a.
data (i));
7731 return read_sparse_matrix<elt_type> (is, a, octave::read_value<double>);
7795 return do_mul_dm_sm<SparseMatrix> (
d, a);
7801 return do_mul_sm_dm<SparseMatrix> (a,
d);
7807 return do_add_dm_sm<SparseMatrix> (
d, a);
7813 return do_sub_dm_sm<SparseMatrix> (
d, a);
7819 return do_add_sm_dm<SparseMatrix> (a,
d);
7825 return do_sub_sm_dm<SparseMatrix> (a,
d);
7844#define EMPTY_RETURN_CHECK(T) \
7845 if (nr == 0 || nc == 0) \
7851 const bool nanflag =
true;
7852 const bool realabs =
true;
7853 return min (
d, m, nanflag, realabs);
7859 const bool realabs =
true;
7860 return min (
d, m, nanflag, realabs);
7880 double tmp = octave::math::min (
d, m.
data (i), nanflag, realabs);
7884 result.
xdata (idx) = tmp;
7894 if (octave::math::min (
d, m.
data (i), nanflag, realabs) != 0.)
7900 result.
xcidx (0) = 0;
7905 double tmp = octave::math::min (
d, m.
data (i), nanflag, realabs);
7909 result.
xdata (ii) = tmp;
7913 result.
xcidx (j+1) = ii;
7923 const bool nanflag =
true;
7924 const bool realabs =
true;
7925 return min (
d, m, nanflag, realabs);
7931 const bool realabs =
true;
7932 return min (
d, m, nanflag, realabs);
7938 return min (
d, m, nanflag, realabs);
7944 const bool nanflag =
true;
7945 const bool realabs =
true;
7946 return min (a, b, nanflag, realabs);
7952 const bool realabs =
true;
7953 return min (a, b, nanflag, realabs);
7967 if (a_nr == b_nr && a_nc == b_nc)
7977 bool ja_lt_max = ja < ja_max;
7981 bool jb_lt_max = jb < jb_max;
7983 while (ja_lt_max || jb_lt_max)
7986 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
7988 double tmp = octave::math::min (a.
data (ja), 0.,
7997 ja_lt_max= ja < ja_max;
7999 else if ((! ja_lt_max)
8000 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
8002 double tmp = octave::math::min (0., b.
data (jb),
8011 jb_lt_max= jb < jb_max;
8015 double tmp = octave::math::min (a.
data (ja), b.
data (jb),
8024 ja_lt_max= ja < ja_max;
8026 jb_lt_max= jb < jb_max;
8036 if (a_nr == 0 && (b_nr == 0 || b_nr == 1))
8038 if (a_nc == 1 || b_nc == 1 || a_nc == b_nc)
8039 r.
resize (a_nr, std::max (a_nc, b_nc));
8041 octave::err_nonconformant (
"min", a_nr, a_nc, b_nr, b_nc);
8043 else if (a_nc == 0 && (b_nc == 0 || b_nc == 1))
8045 if (a_nr == 1 || b_nr == 1 || a_nr == b_nr)
8046 r.
resize (std::max (a_nr, b_nr), a_nc);
8048 octave::err_nonconformant (
"min", a_nr, a_nc, b_nr, b_nc);
8050 else if (b_nr == 0 && (a_nr == 0 || a_nr == 1))
8052 if (b_nc == 1 || a_nc == 1 || b_nc == a_nc)
8053 r.
resize (b_nr, std::max (a_nc, b_nc));
8055 octave::err_nonconformant (
"min", a_nr, a_nc, b_nr, b_nc);
8057 else if (b_nc == 0 && (a_nc == 0 || a_nc == 1))
8059 if (b_nr == 1 || a_nr == 1 || b_nr == a_nr)
8060 r.
resize (std::max (a_nr, b_nr), b_nc);
8062 octave::err_nonconformant (
"min", a_nr, a_nc, b_nr, b_nc);
8065 octave::err_nonconformant (
"min", a_nr, a_nc, b_nr, b_nc);
8074 const bool nanflag =
true;
8075 const bool realabs =
true;
8076 return max (
d, m, nanflag, realabs);
8082 const bool realabs =
true;
8083 return max (
d, m, nanflag, realabs);
8103 double tmp = octave::math::max (
d, m.
data (i), nanflag, realabs);
8108 result.
xdata (idx) = tmp;
8118 if (octave::math::max (
d, m.
data (i), nanflag, realabs) != 0.)
8124 result.
xcidx (0) = 0;
8129 double tmp = octave::math::max (
d, m.
data (i), nanflag, realabs);
8132 result.
xdata (ii) = tmp;
8136 result.
xcidx (j+1) = ii;
8146 const bool nanflag =
true;
8147 const bool realabs =
true;
8148 return max (
d, m, nanflag, realabs);
8154 const bool realabs =
true;
8155 return max (
d, m, nanflag, realabs);
8161 return max (
d, m, nanflag, realabs);
8167 const bool nanflag =
true;
8168 const bool realabs =
true;
8169 return max (a, b, nanflag, realabs);
8175 const bool realabs =
true;
8176 return max (a, b, nanflag, realabs);
8190 if (a_nr == b_nr && a_nc == b_nc)
8200 bool ja_lt_max = ja < ja_max;
8204 bool jb_lt_max = jb < jb_max;
8206 while (ja_lt_max || jb_lt_max)
8209 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
8211 double tmp = octave::math::max (a.
data (ja), 0.,
8220 ja_lt_max= ja < ja_max;
8222 else if ((! ja_lt_max)
8223 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
8225 double tmp = octave::math::max (0., b.
data (jb),
8234 jb_lt_max= jb < jb_max;
8238 double tmp = octave::math::max (a.
data (ja), b.
data (jb),
8247 ja_lt_max= ja < ja_max;
8249 jb_lt_max= jb < jb_max;
8259 if (a_nr == 0 && (b_nr == 0 || b_nr == 1))
8261 if (a_nc == 1 || b_nc == 1 || a_nc == b_nc)
8262 r.
resize (a_nr, std::max (a_nc, b_nc));
8264 octave::err_nonconformant (
"max", a_nr, a_nc, b_nr, b_nc);
8266 else if (a_nc == 0 && (b_nc == 0 || b_nc == 1))
8268 if (a_nr == 1 || b_nr == 1 || a_nr == b_nr)
8269 r.
resize (std::max (a_nr, b_nr), a_nc);
8271 octave::err_nonconformant (
"max", a_nr, a_nc, b_nr, b_nc);
8273 else if (b_nr == 0 && (a_nr == 0 || a_nr == 1))
8275 if (b_nc == 1 || a_nc == 1 || b_nc == a_nc)
8276 r.
resize (b_nr, std::max (a_nc, b_nc));
8278 octave::err_nonconformant (
"max", a_nr, a_nc, b_nr, b_nc);
8280 else if (b_nc == 0 && (a_nc == 0 || a_nc == 1))
8282 if (b_nr == 1 || a_nr == 1 || b_nr == a_nr)
8283 r.
resize (std::max (a_nr, b_nr), b_nc);
8285 octave::err_nonconformant (
"max", a_nr, a_nc, b_nr, b_nc);
8288 octave::err_nonconformant (
"max", a_nr, a_nc, b_nr, b_nc);