26 #if defined (HAVE_CONFIG_H)
43 #include "mx-fcm-fs.h"
45 #include "mx-fs-fcm.h"
69 #if ! defined (USE_QRSOLVE)
121 if (nr != nr_a || nc != nc_a || nz != nz_a)
138 return !(*
this == a);
147 if (nr == nc && nr > 0)
189 return max (dummy_idx, dim);
200 if (dim >= dv.
ndims ())
213 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
224 if (
ridx (i) != idx_j)
261 result.
xcidx (0) = 0;
267 result.
xdata (ii) = tmp;
268 result.
xridx (ii++) = 0;
270 result.
xcidx (j+1) = ii;
277 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
287 if (found[
ridx (i)] == -j)
288 found[
ridx (i)] = -j - 1;
291 if (found[i] > -nc && found[i] < 0)
292 idx_arg.
elem (i) = -found[i];
305 idx_arg.
elem (ir) = j;
311 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
317 result.
xcidx (0) = 0;
318 result.
xcidx (1) = nel;
321 if (idx_arg(j) == -1)
325 result.
xridx (ii++) = j;
332 result.
xdata (ii) = tmp;
333 result.
xridx (ii++) = j;
346 return min (dummy_idx, dim);
357 if (dim >= dv.
ndims ())
370 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
381 if (
ridx (i) != idx_j)
418 result.
xcidx (0) = 0;
424 result.
xdata (ii) = tmp;
425 result.
xridx (ii++) = 0;
427 result.
xcidx (j+1) = ii;
434 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
444 if (found[
ridx (i)] == -j)
445 found[
ridx (i)] = -j - 1;
448 if (found[i] > -nc && found[i] < 0)
449 idx_arg.
elem (i) = -found[i];
462 idx_arg.
elem (ir) = j;
468 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
474 result.
xcidx (0) = 0;
475 result.
xcidx (1) = nel;
478 if (idx_arg(j) == -1)
482 result.
xridx (ii++) = j;
489 result.
xdata (ii) = tmp;
490 result.
xridx (ii++) = j;
568 return insert (tmp , indx);
584 if (rb.
rows () > 0 && rb.
cols () > 0)
594 if (rb.
rows () > 0 && rb.
cols () > 0)
665 return inverse (mattype, info, rcond, 0, 0);
673 return inverse (mattype, info, rcond, 0, 0);
680 return inverse (mattype, info, rcond, 0, 0);
685 double& rcond,
const bool,
686 const bool calccond)
const
694 if (nr == 0 || nc == 0 || nr != nc)
695 (*current_liboctave_error_handler) (
"inverse requires square matrix");
698 int typ = mattype.
type ();
702 (*current_liboctave_error_handler) (
"incorrect matrix type");
735 double& rcond,
const bool,
736 const bool calccond)
const
744 if (nr == 0 || nc == 0 || nr != nc)
745 (*current_liboctave_error_handler) (
"inverse requires square matrix");
748 int typ = mattype.
type ();
753 (*current_liboctave_error_handler) (
"incorrect matrix type");
756 double ainvnorm = 0.;
787 retval.change_capacity (nz2);
805 (*current_liboctave_error_handler) (
"division by zero");
810 rpX =
retval.xridx (colXp);
824 while (rpX < j && rpU < j && colXp < cx && colUp < nz);
828 colUp =
cidx (j+1) - 1;
832 if (pivot == 0. ||
ridx (colUp) != j)
833 (*current_liboctave_error_handler) (
"division by zero");
840 retval.change_capacity (nz2);
844 retval.xdata (cx) = v / pivot;
852 colUp =
cidx (i+1) - 1;
856 if (pivot == 0. ||
ridx (colUp) != i)
857 (*current_liboctave_error_handler) (
"division by zero");
861 retval.xdata (j) /= pivot;
908 k <
cidx (jidx+1); k++)
921 (*current_liboctave_error_handler) (
"division by zero");
929 colUp =
cidx (perm[iidx]+1) - 1;
931 colUp =
cidx (perm[iidx]);
935 (*current_liboctave_error_handler) (
"division by zero");
948 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
949 retval.change_capacity (nz2);
957 retval.xdata (cx++) = work[j];
972 i <
retval.cidx (j+1); i++)
978 rcond = 1. / ainvnorm / anorm;
986 double& rcond,
bool,
bool calc_cond)
const
990 (*current_liboctave_error_handler)
991 (
"inverse of the null matrix not defined");
994 int typ = mattype.
type (
false);
998 typ = mattype.
type (*
this);
1001 ret =
dinverse (mattype, info, rcond,
true, calc_cond);
1015 rcond = fact.
rcond ();
1023 ret =
Q * InvL.
hermitian () * InvL *
Q.transpose ();
1043 rcond = fact.
rcond ();
1051 std::copy_n (
ridx (), nz, ret.
xridx ());
1091 #if defined (HAVE_UMFPACK)
1096 if (nr == 0 || nc == 0 || nr != nc)
1105 Matrix Control (UMFPACK_CONTROL, 1);
1111 Control (UMFPACK_PRL) = tmp;
1116 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
1117 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
1123 Control (UMFPACK_FIXQ) = tmp;
1126 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1137 reinterpret_cast<const double *
> (Ax),
1138 nullptr, 1, control);
1141 Matrix Info (1, UMFPACK_INFO);
1146 reinterpret_cast<const double *
> (Ax),
1147 nullptr,
nullptr, &Symbolic, control, info);
1156 (*current_liboctave_error_handler)
1157 (
"SparseComplexMatrix::determinant symbolic factorization failed");
1167 reinterpret_cast<const double *
> (Ax),
1168 nullptr, Symbolic, &Numeric, control, info);
1171 rcond = Info (UMFPACK_RCOND);
1180 (*current_liboctave_error_handler)
1181 (
"SparseComplexMatrix::determinant numeric factorization failed");
1189 status =
UMFPACK_ZNAME (get_determinant) (c10,
nullptr, &e10,
1197 (*current_liboctave_error_handler)
1198 (
"SparseComplexMatrix::determinant error calculating determinant");
1210 octave_unused_parameter (err);
1211 octave_unused_parameter (rcond);
1213 (*current_liboctave_error_handler)
1214 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
1224 solve_singularity_handler,
bool calc_cond)
const
1233 if (nr != b.
rows ())
1235 (
"matrix dimension mismatch solution of linear equations");
1237 if (nr == 0 || nc == 0 || b.
cols () == 0)
1242 int typ = mattype.
type ();
1246 (*current_liboctave_error_handler) (
"incorrect matrix type");
1271 rcond = dmin / dmax;
1283 solve_singularity_handler,
1284 bool calc_cond)
const
1293 if (nr != b.
rows ())
1295 (
"matrix dimension mismatch solution of linear equations");
1297 if (nr == 0 || nc == 0 || b.
cols () == 0)
1302 int typ = mattype.
type ();
1306 (*current_liboctave_error_handler) (
"incorrect matrix type");
1319 if (b.
ridx (i) >= nm)
1334 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1361 rcond = dmin / dmax;
1373 solve_singularity_handler,
1374 bool calc_cond)
const
1383 if (nr != b.
rows ())
1385 (
"matrix dimension mismatch solution of linear equations");
1387 if (nr == 0 || nc == 0 || b.
cols () == 0)
1392 int typ = mattype.
type ();
1396 (*current_liboctave_error_handler) (
"incorrect matrix type");
1421 rcond = dmin / dmax;
1433 solve_singularity_handler,
1434 bool calc_cond)
const
1443 if (nr != b.
rows ())
1445 (
"matrix dimension mismatch solution of linear equations");
1447 if (nr == 0 || nc == 0 || b.
cols () == 0)
1452 int typ = mattype.
type ();
1456 (*current_liboctave_error_handler) (
"incorrect matrix type");
1469 if (b.
ridx (i) >= nm)
1484 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1511 rcond = dmin / dmax;
1523 solve_singularity_handler sing_handler,
1524 bool calc_cond)
const
1533 if (nr != b.
rows ())
1535 (
"matrix dimension mismatch solution of linear equations");
1537 if (nr == 0 || nc == 0 || b.
cols () == 0)
1542 int typ = mattype.
type ();
1546 (*current_liboctave_error_handler) (
"incorrect matrix type");
1549 double ainvnorm = 0.;
1589 goto triangular_error;
1595 i <
cidx (kidx+1)-1; i++)
1598 work[iidx] = work[iidx] - tmp *
data (i);
1604 retval(perm[i], j) = work[i];
1626 i <
cidx (iidx+1)-1; i++)
1629 work[idx2] = work[idx2] - tmp *
data (i);
1639 if (atmp > ainvnorm)
1642 rcond = 1. / ainvnorm / anorm;
1665 goto triangular_error;
1673 work[iidx] = work[iidx] - tmp *
data (i);
1699 i <
cidx (k+1)-1; i++)
1702 work[iidx] = work[iidx] - tmp *
data (i);
1712 if (atmp > ainvnorm)
1715 rcond = 1. / ainvnorm / anorm;
1724 sing_handler (rcond);
1731 volatile double rcond_plus_one = rcond + 1.0;
1739 sing_handler (rcond);
1753 solve_singularity_handler sing_handler,
1754 bool calc_cond)
const
1763 if (nr != b.
rows ())
1765 (
"matrix dimension mismatch solution of linear equations");
1767 if (nr == 0 || nc == 0 || b.
cols () == 0)
1772 int typ = mattype.
type ();
1776 (*current_liboctave_error_handler) (
"incorrect matrix type");
1779 double ainvnorm = 0.;
1828 goto triangular_error;
1834 i <
cidx (kidx+1)-1; i++)
1837 work[iidx] = work[iidx] - tmp *
data (i);
1849 if (ii + new_nnz > x_nz)
1853 retval.change_capacity (sz);
1858 if (work[rperm[i]] != 0.)
1861 retval.xdata (ii++) = work[rperm[i]];
1866 retval.maybe_compress ();
1887 i <
cidx (iidx+1)-1; i++)
1890 work[idx2] = work[idx2] - tmp *
data (i);
1900 if (atmp > ainvnorm)
1903 rcond = 1. / ainvnorm / anorm;
1925 goto triangular_error;
1933 work[iidx] = work[iidx] - tmp *
data (i);
1945 if (ii + new_nnz > x_nz)
1949 retval.change_capacity (sz);
1957 retval.xdata (ii++) = work[i];
1962 retval.maybe_compress ();
1981 i <
cidx (k+1)-1; i++)
1984 work[iidx] = work[iidx] - tmp *
data (i);
1994 if (atmp > ainvnorm)
1997 rcond = 1. / ainvnorm / anorm;
2006 sing_handler (rcond);
2013 volatile double rcond_plus_one = rcond + 1.0;
2021 sing_handler (rcond);
2034 solve_singularity_handler sing_handler,
2035 bool calc_cond)
const
2044 if (nr != b.
rows ())
2046 (
"matrix dimension mismatch solution of linear equations");
2048 if (nr == 0 || nc == 0 || b.
cols () == 0)
2053 int typ = mattype.
type ();
2057 (*current_liboctave_error_handler) (
"incorrect matrix type");
2060 double ainvnorm = 0.;
2100 goto triangular_error;
2106 i <
cidx (kidx+1)-1; i++)
2109 work[iidx] = work[iidx] - tmp *
data (i);
2115 retval(perm[i], j) = work[i];
2137 i <
cidx (iidx+1)-1; i++)
2140 work[idx2] = work[idx2] - tmp *
data (i);
2150 if (atmp > ainvnorm)
2153 rcond = 1. / ainvnorm / anorm;
2176 goto triangular_error;
2184 work[iidx] = work[iidx] - tmp *
data (i);
2210 i <
cidx (k+1)-1; i++)
2213 work[iidx] = work[iidx] - tmp *
data (i);
2223 if (atmp > ainvnorm)
2226 rcond = 1. / ainvnorm / anorm;
2235 sing_handler (rcond);
2242 volatile double rcond_plus_one = rcond + 1.0;
2250 sing_handler (rcond);
2264 solve_singularity_handler sing_handler,
2265 bool calc_cond)
const
2274 if (nr != b.
rows ())
2276 (
"matrix dimension mismatch solution of linear equations");
2278 if (nr == 0 || nc == 0 || b.
cols () == 0)
2283 int typ = mattype.
type ();
2287 (*current_liboctave_error_handler) (
"incorrect matrix type");
2290 double ainvnorm = 0.;
2339 goto triangular_error;
2345 i <
cidx (kidx+1)-1; i++)
2348 work[iidx] = work[iidx] - tmp *
data (i);
2360 if (ii + new_nnz > x_nz)
2364 retval.change_capacity (sz);
2369 if (work[rperm[i]] != 0.)
2372 retval.xdata (ii++) = work[rperm[i]];
2377 retval.maybe_compress ();
2398 i <
cidx (iidx+1)-1; i++)
2401 work[idx2] = work[idx2] - tmp *
data (i);
2411 if (atmp > ainvnorm)
2414 rcond = 1. / ainvnorm / anorm;
2436 goto triangular_error;
2444 work[iidx] = work[iidx] - tmp *
data (i);
2456 if (ii + new_nnz > x_nz)
2460 retval.change_capacity (sz);
2468 retval.xdata (ii++) = work[i];
2473 retval.maybe_compress ();
2492 i <
cidx (k+1)-1; i++)
2495 work[iidx] = work[iidx] - tmp *
data (i);
2505 if (atmp > ainvnorm)
2508 rcond = 1. / ainvnorm / anorm;
2517 sing_handler (rcond);
2524 volatile double rcond_plus_one = rcond + 1.0;
2532 sing_handler (rcond);
2546 solve_singularity_handler sing_handler,
2547 bool calc_cond)
const
2556 if (nr != b.
rows ())
2558 (
"matrix dimension mismatch solution of linear equations");
2560 if (nr == 0 || nc == 0 || b.
cols () == 0)
2565 int typ = mattype.
type ();
2569 (*current_liboctave_error_handler) (
"incorrect matrix type");
2572 double ainvnorm = 0.;
2600 work[perm[i]] = b(i,j);
2610 if (perm[
ridx (i)] < minr)
2612 minr = perm[
ridx (i)];
2616 if (minr != k ||
data (mini) == 0.)
2619 goto triangular_error;
2630 work[iidx] = work[iidx] - tmp *
data (i);
2657 i <
cidx (k+1); i++)
2658 if (perm[
ridx (i)] < minr)
2660 minr = perm[
ridx (i)];
2667 i <
cidx (k+1); i++)
2673 work[iidx] = work[iidx] - tmp *
data (i);
2684 if (atmp > ainvnorm)
2687 rcond = 1. / ainvnorm / anorm;
2708 goto triangular_error;
2716 work[iidx] = work[iidx] - tmp *
data (i);
2742 i <
cidx (k+1); i++)
2745 work[iidx] = work[iidx] - tmp *
data (i);
2755 if (atmp > ainvnorm)
2758 rcond = 1. / ainvnorm / anorm;
2766 sing_handler (rcond);
2773 volatile double rcond_plus_one = rcond + 1.0;
2781 sing_handler (rcond);
2795 solve_singularity_handler sing_handler,
2796 bool calc_cond)
const
2806 if (nr != b.
rows ())
2808 (
"matrix dimension mismatch solution of linear equations");
2810 if (nr == 0 || nc == 0 || b.
cols () == 0)
2815 int typ = mattype.
type ();
2819 (*current_liboctave_error_handler) (
"incorrect matrix type");
2822 double ainvnorm = 0.;
2855 work[perm[b.
ridx (i)]] = b.
data (i);
2865 if (perm[
ridx (i)] < minr)
2867 minr = perm[
ridx (i)];
2871 if (minr != k ||
data (mini) == 0.)
2874 goto triangular_error;
2885 work[iidx] = work[iidx] - tmp *
data (i);
2897 if (ii + new_nnz > x_nz)
2901 retval.change_capacity (sz);
2909 retval.xdata (ii++) = work[i];
2914 retval.maybe_compress ();
2934 i <
cidx (k+1); i++)
2935 if (perm[
ridx (i)] < minr)
2937 minr = perm[
ridx (i)];
2944 i <
cidx (k+1); i++)
2950 work[iidx] = work[iidx] - tmp *
data (i);
2961 if (atmp > ainvnorm)
2964 rcond = 1. / ainvnorm / anorm;
2985 goto triangular_error;
2993 work[iidx] = work[iidx] - tmp *
data (i);
3005 if (ii + new_nnz > x_nz)
3009 retval.change_capacity (sz);
3017 retval.xdata (ii++) = work[i];
3022 retval.maybe_compress ();
3042 i <
cidx (k+1); i++)
3045 work[iidx] = work[iidx] - tmp *
data (i);
3055 if (atmp > ainvnorm)
3058 rcond = 1. / ainvnorm / anorm;
3067 sing_handler (rcond);
3074 volatile double rcond_plus_one = rcond + 1.0;
3082 sing_handler (rcond);
3096 solve_singularity_handler sing_handler,
3097 bool calc_cond)
const
3106 if (nr != b.
rows ())
3108 (
"matrix dimension mismatch solution of linear equations");
3110 if (nr == 0 || nc == 0 || b.
cols () == 0)
3115 int typ = mattype.
type ();
3119 (*current_liboctave_error_handler) (
"incorrect matrix type");
3122 double ainvnorm = 0.;
3150 work[perm[i]] = b(i,j);
3160 if (perm[
ridx (i)] < minr)
3162 minr = perm[
ridx (i)];
3166 if (minr != k ||
data (mini) == 0.)
3169 goto triangular_error;
3180 work[iidx] = work[iidx] - tmp *
data (i);
3207 i <
cidx (k+1); i++)
3208 if (perm[
ridx (i)] < minr)
3210 minr = perm[
ridx (i)];
3217 i <
cidx (k+1); i++)
3223 work[iidx] = work[iidx] - tmp *
data (i);
3234 if (atmp > ainvnorm)
3237 rcond = 1. / ainvnorm / anorm;
3259 goto triangular_error;
3267 work[iidx] = work[iidx] - tmp *
data (i);
3294 i <
cidx (k+1); i++)
3297 work[iidx] = work[iidx] - tmp *
data (i);
3307 if (atmp > ainvnorm)
3310 rcond = 1. / ainvnorm / anorm;
3319 sing_handler (rcond);
3326 volatile double rcond_plus_one = rcond + 1.0;
3334 sing_handler (rcond);
3348 solve_singularity_handler sing_handler,
3349 bool calc_cond)
const
3358 if (nr != b.
rows ())
3360 (
"matrix dimension mismatch solution of linear equations");
3362 if (nr == 0 || nc == 0 || b.
cols () == 0)
3367 int typ = mattype.
type ();
3371 (*current_liboctave_error_handler) (
"incorrect matrix type");
3374 double ainvnorm = 0.;
3407 work[perm[b.
ridx (i)]] = b.
data (i);
3417 if (perm[
ridx (i)] < minr)
3419 minr = perm[
ridx (i)];
3423 if (minr != k ||
data (mini) == 0.)
3426 goto triangular_error;
3437 work[iidx] = work[iidx] - tmp *
data (i);
3449 if (ii + new_nnz > x_nz)
3453 retval.change_capacity (sz);
3461 retval.xdata (ii++) = work[i];
3466 retval.maybe_compress ();
3486 i <
cidx (k+1); i++)
3487 if (perm[
ridx (i)] < minr)
3489 minr = perm[
ridx (i)];
3496 i <
cidx (k+1); i++)
3502 work[iidx] = work[iidx] - tmp *
data (i);
3513 if (atmp > ainvnorm)
3516 rcond = 1. / ainvnorm / anorm;
3537 goto triangular_error;
3545 work[iidx] = work[iidx] - tmp *
data (i);
3557 if (ii + new_nnz > x_nz)
3561 retval.change_capacity (sz);
3569 retval.xdata (ii++) = work[i];
3574 retval.maybe_compress ();
3594 i <
cidx (k+1); i++)
3597 work[iidx] = work[iidx] - tmp *
data (i);
3607 if (atmp > ainvnorm)
3610 rcond = 1. / ainvnorm / anorm;
3619 sing_handler (rcond);
3626 volatile double rcond_plus_one = rcond + 1.0;
3634 sing_handler (rcond);
3648 solve_singularity_handler sing_handler,
3649 bool calc_cond)
const
3657 if (nr != nc || nr != b.
rows ())
3658 (*current_liboctave_error_handler)
3659 (
"matrix dimension mismatch solution of linear equations");
3661 if (nr == 0 || b.
cols () == 0)
3664 (*current_liboctave_error_handler)
3665 (
"calculation of condition number not implemented");
3669 volatile int typ = mattype.
type ();
3703 else if (
ridx (i) == j + 1)
3714 F77_INT tmp_nr = octave::to_f77_int (nr);
3747 DL[j] =
data (ii++);
3748 DU[j] =
data (ii++);
3750 D[nc-1] =
data (ii);
3767 else if (
ridx (i) == j + 1)
3769 else if (
ridx (i) == j - 1)
3780 F77_INT tmp_nr = octave::to_f77_int (nr);
3797 sing_handler (rcond);
3808 (*current_liboctave_error_handler) (
"incorrect matrix type");
3817 solve_singularity_handler sing_handler,
3818 bool calc_cond)
const
3826 if (nr != nc || nr != b.
rows ())
3827 (*current_liboctave_error_handler)
3828 (
"matrix dimension mismatch solution of linear equations");
3830 if (nr == 0 || b.
cols () == 0)
3833 (*current_liboctave_error_handler)
3834 (
"calculation of condition number not implemented");
3838 int typ = mattype.
type ();
3859 DL[j] =
data (ii++);
3860 DU[j] =
data (ii++);
3862 D[nc-1] =
data (ii);
3879 else if (
ridx (i) == j + 1)
3881 else if (
ridx (i) == j - 1)
3886 F77_INT tmp_nr = octave::to_f77_int (nr);
3905 sing_handler (rcond);
3932 (F77_CONST_CHAR_ARG2 (&job, 1),
3938 F77_CHAR_ARG_LEN (1)));
3949 if (ii + new_nnz > x_nz)
3953 retval.change_capacity (sz);
3961 retval.xdata (ii++) = work[i];
3966 retval.maybe_compress ();
3970 (*current_liboctave_error_handler) (
"incorrect matrix type");
3979 solve_singularity_handler sing_handler,
3980 bool calc_cond)
const
3988 if (nr != nc || nr != b.
rows ())
3989 (*current_liboctave_error_handler)
3990 (
"matrix dimension mismatch solution of linear equations");
3992 if (nr == 0 || b.
cols () == 0)
3995 (*current_liboctave_error_handler)
3996 (
"calculation of condition number not implemented");
4000 volatile int typ = mattype.
type ();
4034 else if (
ridx (i) == j + 1)
4047 F77_INT tmp_nr = octave::to_f77_int (nr);
4078 DL[j] =
data (ii++);
4079 DU[j] =
data (ii++);
4081 D[nc-1] =
data (ii);
4098 else if (
ridx (i) == j + 1)
4100 else if (
ridx (i) == j - 1)
4113 F77_INT tmp_nr = octave::to_f77_int (nr);
4130 sing_handler (rcond);
4138 (*current_liboctave_error_handler) (
"incorrect matrix type");
4148 solve_singularity_handler sing_handler,
4149 bool calc_cond)
const
4157 if (nr != nc || nr != b.
rows ())
4158 (*current_liboctave_error_handler)
4159 (
"matrix dimension mismatch solution of linear equations");
4161 if (nr == 0 || b.
cols () == 0)
4164 (*current_liboctave_error_handler)
4165 (
"calculation of condition number not implemented");
4169 int typ = mattype.
type ();
4190 DL[j] =
data (ii++);
4191 DU[j] =
data (ii++);
4193 D[nc-1] =
data (ii);
4210 else if (
ridx (i) == j + 1)
4212 else if (
ridx (i) == j - 1)
4217 F77_INT tmp_nr = octave::to_f77_int (nr);
4236 sing_handler (rcond);
4260 for (
F77_INT i = 0; i < b_nr; i++)
4264 (F77_CONST_CHAR_ARG2 (&job, 1),
4270 F77_CHAR_ARG_LEN (1)));
4278 (*current_liboctave_error_handler)
4279 (
"SparseComplexMatrix::solve solve failed");
4292 if (ii + new_nnz > x_nz)
4296 retval.change_capacity (sz);
4304 retval.xdata (ii++) = Bx[i];
4310 retval.maybe_compress ();
4314 (*current_liboctave_error_handler) (
"incorrect matrix type");
4323 solve_singularity_handler sing_handler,
4324 bool calc_cond)
const
4332 if (nr != nc || nr != b.
rows ())
4333 (*current_liboctave_error_handler)
4334 (
"matrix dimension mismatch solution of linear equations");
4336 if (nr == 0 || b.
cols () == 0)
4341 volatile int typ = mattype.
type ();
4355 for (
F77_INT j = 0; j < ldm; j++)
4357 tmp_data[ii++] = 0.;
4365 m_band(ri - j, j) =
data (i);
4373 F77_INT tmp_nr = octave::to_f77_int (nr);
4378 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4380 F77_CHAR_ARG_LEN (1)));
4403 (F77_CONST_CHAR_ARG2 (&job, 1),
4406 F77_CHAR_ARG_LEN (1)));
4413 volatile double rcond_plus_one = rcond + 1.0;
4421 sing_handler (rcond);
4440 (F77_CONST_CHAR_ARG2 (&job, 1),
4443 F77_CHAR_ARG_LEN (1)));
4450 (*current_liboctave_error_handler)
4451 (
"SparseMatrix::solve solve failed");
4463 F77_INT ldm = n_upper + 2 * n_lower + 1;
4472 for (
F77_INT j = 0; j < ldm; j++)
4474 tmp_data[ii++] = 0.;
4479 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
4498 F77_INT tmp_nr = octave::to_f77_int (nr);
4502 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4504 ldm, pipvt, tmp_err));
4517 sing_handler (rcond);
4533 F77_INT tmp_nc = octave::to_f77_int (nc);
4536 (F77_CONST_CHAR_ARG2 (&job, 1),
4539 F77_CHAR_ARG_LEN (1)));
4546 volatile double rcond_plus_one = rcond + 1.0;
4554 sing_handler (rcond);
4574 (F77_CONST_CHAR_ARG2 (&job, 1),
4577 F77_CHAR_ARG_LEN (1)));
4584 (*current_liboctave_error_handler) (
"incorrect matrix type");
4593 solve_singularity_handler sing_handler,
4594 bool calc_cond)
const
4602 if (nr != nc || nr != b.
rows ())
4603 (*current_liboctave_error_handler)
4604 (
"matrix dimension mismatch solution of linear equations");
4606 if (nr == 0 || b.
cols () == 0)
4611 volatile int typ = mattype.
type ();
4626 for (
F77_INT j = 0; j < ldm; j++)
4628 tmp_data[ii++] = 0.;
4636 m_band(ri - j, j) =
data (i);
4644 F77_INT tmp_nr = octave::to_f77_int (nr);
4649 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4651 F77_CHAR_ARG_LEN (1)));
4672 (F77_CONST_CHAR_ARG2 (&job, 1),
4675 F77_CHAR_ARG_LEN (1)));
4682 volatile double rcond_plus_one = rcond + 1.0;
4690 sing_handler (rcond);
4715 for (
F77_INT i = 0; i < b_nr; i++)
4716 Bx[i] = b.
elem (i, j);
4719 (F77_CONST_CHAR_ARG2 (&job, 1),
4722 F77_CHAR_ARG_LEN (1)));
4729 (*current_liboctave_error_handler)
4730 (
"SparseComplexMatrix::solve solve failed");
4744 sz = (
static_cast<double> (b_nc) - j) / b_nc
4746 sz = x_nz + (sz > 100 ? sz : 100);
4747 retval.change_capacity (sz);
4757 retval.maybe_compress ();
4767 F77_INT ldm = n_upper + 2 * n_lower + 1;
4776 for (
F77_INT j = 0; j < ldm; j++)
4778 tmp_data[ii++] = 0.;
4783 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
4802 F77_INT tmp_nr = octave::to_f77_int (nr);
4806 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4808 ldm, pipvt, tmp_err));
4819 sing_handler (rcond);
4835 F77_INT tmp_nc = octave::to_f77_int (nc);
4838 (F77_CONST_CHAR_ARG2 (&job, 1),
4841 F77_CHAR_ARG_LEN (1)));
4848 volatile double rcond_plus_one = rcond + 1.0;
4856 sing_handler (rcond);
4883 i < b.
cidx (j+1); i++)
4887 (F77_CONST_CHAR_ARG2 (&job, 1),
4890 F77_CHAR_ARG_LEN (1)));
4901 if (ii + new_nnz > x_nz)
4905 retval.change_capacity (sz);
4913 retval.xdata (ii++) = work[i];
4918 retval.maybe_compress ();
4923 (*current_liboctave_error_handler) (
"incorrect matrix type");
4932 solve_singularity_handler sing_handler,
4933 bool calc_cond)
const
4941 if (nr != nc || nr != b.
rows ())
4942 (*current_liboctave_error_handler)
4943 (
"matrix dimension mismatch solution of linear equations");
4945 if (nr == 0 || b.
cols () == 0)
4950 volatile int typ = mattype.
type ();
4965 for (
F77_INT j = 0; j < ldm; j++)
4967 tmp_data[ii++] = 0.;
4975 m_band(ri - j, j) =
data (i);
4983 F77_INT tmp_nr = octave::to_f77_int (nr);
4988 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4990 F77_CHAR_ARG_LEN (1)));
5013 (F77_CONST_CHAR_ARG2 (&job, 1),
5016 F77_CHAR_ARG_LEN (1)));
5023 volatile double rcond_plus_one = rcond + 1.0;
5031 sing_handler (rcond);
5049 (F77_CONST_CHAR_ARG2 (&job, 1),
5052 F77_CHAR_ARG_LEN (1)));
5059 (*current_liboctave_error_handler)
5060 (
"SparseComplexMatrix::solve solve failed");
5072 F77_INT ldm = n_upper + 2 * n_lower + 1;
5081 for (
F77_INT j = 0; j < ldm; j++)
5083 tmp_data[ii++] = 0.;
5088 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
5107 F77_INT tmp_nr = octave::to_f77_int (nr);
5111 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5113 ldm, pipvt, tmp_err));
5124 sing_handler (rcond);
5140 F77_INT tmp_nc = octave::to_f77_int (nc);
5143 (F77_CONST_CHAR_ARG2 (&job, 1),
5146 F77_CHAR_ARG_LEN (1)));
5153 volatile double rcond_plus_one = rcond + 1.0;
5161 sing_handler (rcond);
5180 (F77_CONST_CHAR_ARG2 (&job, 1),
5183 F77_CHAR_ARG_LEN (1)));
5190 (*current_liboctave_error_handler) (
"incorrect matrix type");
5199 solve_singularity_handler sing_handler,
5200 bool calc_cond)
const
5208 if (nr != nc || nr != b.
rows ())
5209 (*current_liboctave_error_handler)
5210 (
"matrix dimension mismatch solution of linear equations");
5212 if (nr == 0 || b.
cols () == 0)
5217 volatile int typ = mattype.
type ();
5232 for (
F77_INT j = 0; j < ldm; j++)
5234 tmp_data[ii++] = 0.;
5242 m_band(ri - j, j) =
data (i);
5250 F77_INT tmp_nr = octave::to_f77_int (nr);
5255 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5257 F77_CHAR_ARG_LEN (1)));
5281 (F77_CONST_CHAR_ARG2 (&job, 1),
5284 F77_CHAR_ARG_LEN (1)));
5291 volatile double rcond_plus_one = rcond + 1.0;
5299 sing_handler (rcond);
5325 for (
F77_INT i = 0; i < b_nr; i++)
5329 (F77_CONST_CHAR_ARG2 (&job, 1),
5332 F77_CHAR_ARG_LEN (1)));
5339 (*current_liboctave_error_handler)
5340 (
"SparseMatrix::solve solve failed");
5352 if (ii + new_nnz > x_nz)
5356 retval.change_capacity (sz);
5364 retval.xdata (ii++) = Bx[i];
5370 retval.maybe_compress ();
5380 F77_INT ldm = n_upper + 2 * n_lower + 1;
5389 for (
F77_INT j = 0; j < ldm; j++)
5391 tmp_data[ii++] = 0.;
5396 m_band(
ridx (i) - j + n_lower + n_upper, j) =
data (i);
5415 F77_INT tmp_nr = octave::to_f77_int (nr);
5419 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5421 ldm, pipvt, tmp_err));
5432 sing_handler (rcond);
5448 F77_INT tmp_nc = octave::to_f77_int (nc);
5451 (F77_CONST_CHAR_ARG2 (&job, 1),
5454 F77_CHAR_ARG_LEN (1)));
5461 volatile double rcond_plus_one = rcond + 1.0;
5469 sing_handler (rcond);
5497 i < b.
cidx (j+1); i++)
5501 (F77_CONST_CHAR_ARG2 (&job, 1),
5504 F77_CHAR_ARG_LEN (1)));
5515 if (ii + new_nnz > x_nz)
5519 retval.change_capacity (sz);
5527 retval.xdata (ii++) = Bx[i];
5532 retval.maybe_compress ();
5537 (*current_liboctave_error_handler) (
"incorrect matrix type");
5546 solve_singularity_handler sing_handler,
5547 bool calc_cond)
const
5550 void *Numeric =
nullptr;
5553 #if defined (HAVE_UMFPACK)
5556 Control =
Matrix (UMFPACK_CONTROL, 1);
5562 Control (UMFPACK_PRL) = tmp;
5566 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
5567 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
5573 Control (UMFPACK_FIXQ) = tmp;
5586 reinterpret_cast<const double *
> (Ax),
5587 nullptr, 1, control);
5590 Info =
Matrix (1, UMFPACK_INFO);
5595 reinterpret_cast<const double *
> (Ax),
5596 nullptr,
nullptr, &Symbolic, control, info);
5606 (*current_liboctave_error_handler)
5607 (
"SparseComplexMatrix::solve symbolic factorization failed");
5616 reinterpret_cast<const double *
> (Ax),
5617 nullptr, Symbolic, &Numeric, control, info);
5621 rcond = Info (UMFPACK_RCOND);
5624 volatile double rcond_plus_one = rcond + 1.0;
5626 if (status == UMFPACK_WARNING_singular_matrix
5634 sing_handler (rcond);
5638 else if (status < 0)
5644 (*current_liboctave_error_handler)
5645 (
"SparseComplexMatrix::solve numeric factorization failed");
5660 octave_unused_parameter (rcond);
5661 octave_unused_parameter (Control);
5662 octave_unused_parameter (Info);
5663 octave_unused_parameter (sing_handler);
5664 octave_unused_parameter (calc_cond);
5666 (*current_liboctave_error_handler)
5667 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
5677 solve_singularity_handler sing_handler,
5678 bool calc_cond)
const
5686 if (nr != nc || nr != b.
rows ())
5687 (*current_liboctave_error_handler)
5688 (
"matrix dimension mismatch solution of linear equations");
5690 if (nr == 0 || b.
cols () == 0)
5695 volatile int typ = mattype.
type ();
5700 #if defined (HAVE_CHOLMOD)
5701 cholmod_common Common;
5702 cholmod_common *cm = &Common;
5706 cm->prefer_zomplex =
false;
5712 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
5716 cm->print =
static_cast<int> (spu) + 2;
5717 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5721 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5722 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5724 cm->final_ll =
true;
5726 cholmod_sparse Astore;
5727 cholmod_sparse *
A = &Astore;
5737 #if defined (OCTAVE_ENABLE_64)
5738 A->itype = CHOLMOD_LONG;
5740 A->itype = CHOLMOD_INT;
5742 A->dtype = CHOLMOD_DOUBLE;
5744 A->xtype = CHOLMOD_COMPLEX;
5748 cholmod_dense Bstore;
5749 cholmod_dense *
B = &Bstore;
5750 B->nrow = b.
rows ();
5751 B->ncol = b.
cols ();
5753 B->nzmax =
B->nrow *
B->ncol;
5754 B->dtype = CHOLMOD_DOUBLE;
5755 B->xtype = CHOLMOD_REAL;
5777 volatile double rcond_plus_one = rcond + 1.0;
5785 sing_handler (rcond);
5811 static char blank_name[] =
" ";
5816 (*current_liboctave_warning_with_id_handler)
5817 (
"Octave:missing-dependency",
5818 "support for CHOLMOD was unavailable or disabled "
5819 "when liboctave was built");
5828 #if defined (HAVE_UMFPACK)
5830 void *Numeric =
factorize (err, rcond, Control, Info,
5831 sing_handler, calc_cond);
5836 Control (UMFPACK_IRSTEP) = 1;
5845 #if defined (UMFPACK_SEPARATE_SPLIT)
5858 #if defined (UMFPACK_SEPARATE_SPLIT)
5862 reinterpret_cast<const double *
> (Ax),
5864 reinterpret_cast<double *
> (&Xx[iidx]),
5866 &Bx[iidx], Bz, Numeric,
5870 Bz[i] = b.
elem (i, j);
5875 reinterpret_cast<const double *
> (Ax),
5877 reinterpret_cast<double *
> (&Xx[iidx]),
5879 reinterpret_cast<const double *
> (Bz),
5889 (*current_liboctave_error_handler)
5890 (
"SparseComplexMatrix::solve solve failed");
5905 octave_unused_parameter (rcond);
5906 octave_unused_parameter (sing_handler);
5907 octave_unused_parameter (calc_cond);
5909 (*current_liboctave_error_handler)
5910 (
"support for UMFPACK was unavailable or disabled "
5911 "when liboctave was built");
5915 (*current_liboctave_error_handler) (
"incorrect matrix type");
5924 solve_singularity_handler sing_handler,
5925 bool calc_cond)
const
5933 if (nr != nc || nr != b.
rows ())
5934 (*current_liboctave_error_handler)
5935 (
"matrix dimension mismatch solution of linear equations");
5937 if (nr == 0 || b.
cols () == 0)
5942 volatile int typ = mattype.
type ();
5947 #if defined (HAVE_CHOLMOD)
5948 cholmod_common Common;
5949 cholmod_common *cm = &Common;
5953 cm->prefer_zomplex =
false;
5959 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
5963 cm->print =
static_cast<int> (spu) + 2;
5964 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5968 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5969 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5971 cm->final_ll =
true;
5973 cholmod_sparse Astore;
5974 cholmod_sparse *
A = &Astore;
5984 #if defined (OCTAVE_ENABLE_64)
5985 A->itype = CHOLMOD_LONG;
5987 A->itype = CHOLMOD_INT;
5989 A->dtype = CHOLMOD_DOUBLE;
5991 A->xtype = CHOLMOD_COMPLEX;
5995 cholmod_sparse Bstore;
5996 cholmod_sparse *
B = &Bstore;
5997 B->nrow = b.
rows ();
5998 B->ncol = b.
cols ();
6001 B->nzmax = b.
nnz ();
6005 #if defined (OCTAVE_ENABLE_64)
6006 B->itype = CHOLMOD_LONG;
6008 B->itype = CHOLMOD_INT;
6010 B->dtype = CHOLMOD_DOUBLE;
6012 B->xtype = CHOLMOD_REAL;
6034 volatile double rcond_plus_one = rcond + 1.0;
6042 sing_handler (rcond);
6061 j <= static_cast<octave_idx_type> (X->ncol); j++)
6064 j < static_cast<octave_idx_type> (X->nzmax); j++)
6074 static char blank_name[] =
" ";
6079 (*current_liboctave_warning_with_id_handler)
6080 (
"Octave:missing-dependency",
6081 "support for CHOLMOD was unavailable or disabled "
6082 "when liboctave was built");
6091 #if defined (HAVE_UMFPACK)
6093 void *Numeric =
factorize (err, rcond, Control, Info,
6094 sing_handler, calc_cond);
6099 Control (UMFPACK_IRSTEP) = 1;
6109 #if defined (UMFPACK_SEPARATE_SPLIT)
6130 #if defined (UMFPACK_SEPARATE_SPLIT)
6132 Bx[i] = b.
elem (i, j);
6137 reinterpret_cast<const double *
> (Ax),
6139 reinterpret_cast<double *
> (Xx),
6141 Bx, Bz, Numeric, control,
6145 Bz[i] = b.
elem (i, j);
6150 reinterpret_cast<const double *
> (Ax),
6152 reinterpret_cast<double *
> (Xx),
6154 reinterpret_cast<double *
> (Bz),
6164 (*current_liboctave_error_handler)
6165 (
"SparseComplexMatrix::solve solve failed");
6180 sz = (
static_cast<double> (b_nc) - j) / b_nc
6182 sz = x_nz + (sz > 100 ? sz : 100);
6183 retval.change_capacity (sz);
6193 retval.maybe_compress ();
6203 octave_unused_parameter (rcond);
6204 octave_unused_parameter (sing_handler);
6205 octave_unused_parameter (calc_cond);
6207 (*current_liboctave_error_handler)
6208 (
"support for UMFPACK was unavailable or disabled "
6209 "when liboctave was built");
6213 (*current_liboctave_error_handler) (
"incorrect matrix type");
6222 solve_singularity_handler sing_handler,
6223 bool calc_cond)
const
6231 if (nr != nc || nr != b.
rows ())
6232 (*current_liboctave_error_handler)
6233 (
"matrix dimension mismatch solution of linear equations");
6235 if (nr == 0 || b.
cols () == 0)
6240 volatile int typ = mattype.
type ();
6245 #if defined (HAVE_CHOLMOD)
6246 cholmod_common Common;
6247 cholmod_common *cm = &Common;
6251 cm->prefer_zomplex =
false;
6257 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6261 cm->print =
static_cast<int> (spu) + 2;
6262 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6266 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6267 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6269 cm->final_ll =
true;
6271 cholmod_sparse Astore;
6272 cholmod_sparse *
A = &Astore;
6282 #if defined (OCTAVE_ENABLE_64)
6283 A->itype = CHOLMOD_LONG;
6285 A->itype = CHOLMOD_INT;
6287 A->dtype = CHOLMOD_DOUBLE;
6289 A->xtype = CHOLMOD_COMPLEX;
6293 cholmod_dense Bstore;
6294 cholmod_dense *
B = &Bstore;
6295 B->nrow = b.
rows ();
6296 B->ncol = b.
cols ();
6298 B->nzmax =
B->nrow *
B->ncol;
6299 B->dtype = CHOLMOD_DOUBLE;
6300 B->xtype = CHOLMOD_COMPLEX;
6322 volatile double rcond_plus_one = rcond + 1.0;
6330 sing_handler (rcond);
6356 static char blank_name[] =
" ";
6361 (*current_liboctave_warning_with_id_handler)
6362 (
"Octave:missing-dependency",
6363 "support for CHOLMOD was unavailable or disabled "
6364 "when liboctave was built");
6373 #if defined (HAVE_UMFPACK)
6375 void *Numeric =
factorize (err, rcond, Control, Info,
6376 sing_handler, calc_cond);
6381 Control (UMFPACK_IRSTEP) = 1;
6401 reinterpret_cast<const double *
> (Ax),
6403 reinterpret_cast<double *
> (&Xx[iidx]),
6405 reinterpret_cast<const double *
> (&Bx[iidx]),
6406 nullptr, Numeric, control, info);
6413 (*current_liboctave_error_handler)
6414 (
"SparseComplexMatrix::solve solve failed");
6429 octave_unused_parameter (rcond);
6430 octave_unused_parameter (sing_handler);
6431 octave_unused_parameter (calc_cond);
6433 (*current_liboctave_error_handler)
6434 (
"support for UMFPACK was unavailable or disabled "
6435 "when liboctave was built");
6439 (*current_liboctave_error_handler) (
"incorrect matrix type");
6448 solve_singularity_handler sing_handler,
6449 bool calc_cond)
const
6457 if (nr != nc || nr != b.
rows ())
6458 (*current_liboctave_error_handler)
6459 (
"matrix dimension mismatch solution of linear equations");
6461 if (nr == 0 || b.
cols () == 0)
6466 volatile int typ = mattype.
type ();
6471 #if defined (HAVE_CHOLMOD)
6472 cholmod_common Common;
6473 cholmod_common *cm = &Common;
6477 cm->prefer_zomplex =
false;
6483 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6487 cm->print =
static_cast<int> (spu) + 2;
6488 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6492 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6493 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6495 cm->final_ll =
true;
6497 cholmod_sparse Astore;
6498 cholmod_sparse *
A = &Astore;
6508 #if defined (OCTAVE_ENABLE_64)
6509 A->itype = CHOLMOD_LONG;
6511 A->itype = CHOLMOD_INT;
6513 A->dtype = CHOLMOD_DOUBLE;
6515 A->xtype = CHOLMOD_COMPLEX;
6519 cholmod_sparse Bstore;
6520 cholmod_sparse *
B = &Bstore;
6521 B->nrow = b.
rows ();
6522 B->ncol = b.
cols ();
6525 B->nzmax = b.
nnz ();
6529 #if defined (OCTAVE_ENABLE_64)
6530 B->itype = CHOLMOD_LONG;
6532 B->itype = CHOLMOD_INT;
6534 B->dtype = CHOLMOD_DOUBLE;
6536 B->xtype = CHOLMOD_COMPLEX;
6558 volatile double rcond_plus_one = rcond + 1.0;
6566 sing_handler (rcond);
6585 j <= static_cast<octave_idx_type> (X->ncol); j++)
6588 j < static_cast<octave_idx_type> (X->nzmax); j++)
6598 static char blank_name[] =
" ";
6603 (*current_liboctave_warning_with_id_handler)
6604 (
"Octave:missing-dependency",
6605 "support for CHOLMOD was unavailable or disabled "
6606 "when liboctave was built");
6615 #if defined (HAVE_UMFPACK)
6617 void *Numeric =
factorize (err, rcond, Control, Info,
6618 sing_handler, calc_cond);
6623 Control (UMFPACK_IRSTEP) = 1;
6652 reinterpret_cast<const double *
> (Ax),
6654 reinterpret_cast<double *
> (Xx),
6656 reinterpret_cast<double *
> (Bx),
6657 nullptr, Numeric, control, info);
6664 (*current_liboctave_error_handler)
6665 (
"SparseComplexMatrix::solve solve failed");
6680 sz = (
static_cast<double> (b_nc) - j) / b_nc
6682 sz = x_nz + (sz > 100 ? sz : 100);
6683 retval.change_capacity (sz);
6693 retval.maybe_compress ();
6695 rcond = Info (UMFPACK_RCOND);
6696 volatile double rcond_plus_one = rcond + 1.0;
6698 if (status == UMFPACK_WARNING_singular_matrix
6704 sing_handler (rcond);
6717 octave_unused_parameter (rcond);
6718 octave_unused_parameter (sing_handler);
6719 octave_unused_parameter (calc_cond);
6721 (*current_liboctave_error_handler)
6722 (
"support for UMFPACK was unavailable or disabled "
6723 "when liboctave was built");
6727 (*current_liboctave_error_handler) (
"incorrect matrix type");
6738 return solve (mattype, b, info, rcond,
nullptr);
6746 return solve (mattype, b, info, rcond,
nullptr);
6753 return solve (mattype, b, info, rcond,
nullptr);
6759 solve_singularity_handler sing_handler,
6760 bool singular_fallback)
const
6763 int typ = mattype.
type (
false);
6766 typ = mattype.
type (*
this);
6769 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6771 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6773 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6775 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6780 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6782 (*current_liboctave_error_handler) (
"unknown matrix type");
6787 #if defined (USE_QRSOLVE)
6803 return solve (mattype, b, info, rcond,
nullptr);
6811 return solve (mattype, b, info, rcond,
nullptr);
6818 return solve (mattype, b, info, rcond,
nullptr);
6824 solve_singularity_handler sing_handler,
6825 bool singular_fallback)
const
6828 int typ = mattype.
type (
false);
6831 typ = mattype.
type (*
this);
6834 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6836 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6838 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6840 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6845 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6847 (*current_liboctave_error_handler) (
"unknown matrix type");
6852 #if defined (USE_QRSOLVE)
6868 return solve (mattype, b, info, rcond,
nullptr);
6876 return solve (mattype, b, info, rcond,
nullptr);
6883 return solve (mattype, b, info, rcond,
nullptr);
6889 solve_singularity_handler sing_handler,
6890 bool singular_fallback)
const
6893 int typ = mattype.
type (
false);
6896 typ = mattype.
type (*
this);
6899 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6901 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6903 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6905 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6910 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6912 (*current_liboctave_error_handler) (
"unknown matrix type");
6917 #if defined (USE_QRSOLVE)
6934 return solve (mattype, b, info, rcond,
nullptr);
6942 return solve (mattype, b, info, rcond,
nullptr);
6949 return solve (mattype, b, info, rcond,
nullptr);
6955 solve_singularity_handler sing_handler,
6956 bool singular_fallback)
const
6959 int typ = mattype.
type (
false);
6962 typ = mattype.
type (*
this);
6965 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6967 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6969 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6971 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6976 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6978 (*current_liboctave_error_handler) (
"unknown matrix type");
6983 #if defined (USE_QRSOLVE)
6998 return solve (mattype, b, info, rcond);
7006 return solve (mattype, b, info, rcond);
7013 return solve (mattype, b, info, rcond,
nullptr);
7019 solve_singularity_handler sing_handler)
const
7022 return solve (mattype, tmp, info, rcond,
7032 return solve (mattype, b, info, rcond,
nullptr);
7040 return solve (mattype, b, info, rcond,
nullptr);
7047 return solve (mattype, b, info, rcond,
nullptr);
7053 solve_singularity_handler sing_handler)
const
7056 return solve (mattype, tmp, info, rcond,
7065 return solve (b, info, rcond,
nullptr);
7072 return solve (b, info, rcond,
nullptr);
7077 double& rcond)
const
7079 return solve (b, info, rcond,
nullptr);
7085 solve_singularity_handler sing_handler)
const
7088 return solve (mattype, b, err, rcond, sing_handler);
7096 return solve (b, info, rcond,
nullptr);
7104 return solve (b, info, rcond,
nullptr);
7111 return solve (b, info, rcond,
nullptr);
7117 solve_singularity_handler sing_handler)
const
7120 return solve (mattype, b, err, rcond, sing_handler);
7128 return solve (b, info, rcond,
nullptr);
7135 return solve (b, info, rcond,
nullptr);
7141 solve_singularity_handler sing_handler)
const
7144 return solve (mattype, b, err, rcond, sing_handler);
7152 return solve (b, info, rcond,
nullptr);
7160 return solve (b, info, rcond,
nullptr);
7167 return solve (b, info, rcond,
nullptr);
7173 solve_singularity_handler sing_handler)
const
7176 return solve (mattype, b, err, rcond, sing_handler);
7183 return solve (b, info, rcond);
7190 return solve (b, info, rcond);
7195 double& rcond)
const
7197 return solve (b, info, rcond,
nullptr);
7203 solve_singularity_handler sing_handler)
const
7206 return solve (tmp, info, rcond,
7215 return solve (b, info, rcond,
nullptr);
7223 return solve (b, info, rcond,
nullptr);
7228 double& rcond)
const
7230 return solve (b, info, rcond,
nullptr);
7236 solve_singularity_handler sing_handler)
const
7239 return solve (tmp, info, rcond,
7264 if (jj <
cidx (i+1) &&
ridx (jj) == j)
7361 double r_val = val.real ();
7362 double i_val = val.imag ();
7364 if (r_val > max_val)
7367 if (i_val > max_val)
7370 if (r_val < min_val)
7373 if (i_val < min_val)
7419 if ((
rows () == 1 && dim == -1) || dim == 1)
7424 (
cidx (j+1) -
cidx (j) < nr ? 0.0 : 1.0), 1.0);
7438 Complex d = data (i); \
7439 tmp[ridx (i)] += d * conj (d)
7442 Complex d = data (i); \
7443 tmp[j] += d * conj (d)
7489 os << a.
ridx (i) + 1 <<
' ' << j + 1 <<
' ';
7503 return read_sparse_matrix<elt_type> (is, a, octave_read_value<Complex>);
7588 return do_mul_dm_sm<SparseComplexMatrix> (
d, a);
7593 return do_mul_sm_dm<SparseComplexMatrix> (a,
d);
7599 return do_mul_dm_sm<SparseComplexMatrix> (
d, a);
7604 return do_mul_sm_dm<SparseComplexMatrix> (a,
d);
7610 return do_mul_dm_sm<SparseComplexMatrix> (
d, a);
7615 return do_mul_sm_dm<SparseComplexMatrix> (a,
d);
7621 return do_add_dm_sm<SparseComplexMatrix> (
d, a);
7626 return do_add_dm_sm<SparseComplexMatrix> (
d, a);
7631 return do_add_dm_sm<SparseComplexMatrix> (
d, a);
7636 return do_add_sm_dm<SparseComplexMatrix> (a,
d);
7641 return do_add_sm_dm<SparseComplexMatrix> (a,
d);
7646 return do_add_sm_dm<SparseComplexMatrix> (a,
d);
7652 return do_sub_dm_sm<SparseComplexMatrix> (
d, a);
7657 return do_sub_dm_sm<SparseComplexMatrix> (
d, a);
7662 return do_sub_dm_sm<SparseComplexMatrix> (
d, a);
7667 return do_sub_sm_dm<SparseComplexMatrix> (a,
d);
7672 return do_sub_sm_dm<SparseComplexMatrix> (a,
d);
7677 return do_sub_sm_dm<SparseComplexMatrix> (a,
d);
7696 #define EMPTY_RETURN_CHECK(T) \
7697 if (nr == 0 || nc == 0) \
7740 if (a_nr == b_nr && a_nc == b_nc)
7750 bool ja_lt_max = ja < ja_max;
7754 bool jb_lt_max = jb < jb_max;
7756 while (ja_lt_max || jb_lt_max)
7759 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
7764 r.ridx (jx) = a.
ridx (ja);
7769 ja_lt_max= ja < ja_max;
7771 else if ((! ja_lt_max)
7772 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
7777 r.ridx (jx) = b.
ridx (jb);
7782 jb_lt_max= jb < jb_max;
7790 r.ridx (jx) = a.
ridx (ja);
7794 ja_lt_max= ja < ja_max;
7796 jb_lt_max= jb < jb_max;
7802 r.maybe_compress ();
7806 if (a_nr == 0 || a_nc == 0)
7807 r.resize (a_nr, a_nc);
7808 else if (b_nr == 0 || b_nc == 0)
7809 r.resize (b_nr, b_nc);
7857 if (a_nr == b_nr && a_nc == b_nc)
7867 bool ja_lt_max = ja < ja_max;
7871 bool jb_lt_max = jb < jb_max;
7873 while (ja_lt_max || jb_lt_max)
7876 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
7881 r.ridx (jx) = a.
ridx (ja);
7886 ja_lt_max= ja < ja_max;
7888 else if ((! ja_lt_max)
7889 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
7894 r.ridx (jx) = b.
ridx (jb);
7899 jb_lt_max= jb < jb_max;
7907 r.ridx (jx) = a.
ridx (ja);
7911 ja_lt_max= ja < ja_max;
7913 jb_lt_max= jb < jb_max;
7919 r.maybe_compress ();
7923 if (a_nr == 0 || a_nc == 0)
7924 r.resize (a_nr, a_nc);
7925 else if (b_nr == 0 || b_nc == 0)
7926 r.resize (b_nr, b_nc);
SparseComplexMatrix max(const Complex &c, const SparseComplexMatrix &m)
std::istream & operator>>(std::istream &is, SparseComplexMatrix &a)
std::ostream & operator<<(std::ostream &os, const SparseComplexMatrix &a)
SparseComplexMatrix operator-(const ComplexDiagMatrix &d, const SparseMatrix &a)
SparseComplexMatrix conj(const SparseComplexMatrix &a)
static const Complex Complex_NaN_result(octave::numeric_limits< double >::NaN(), octave::numeric_limits< double >::NaN())
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)
SparseComplexMatrix operator+(const ComplexDiagMatrix &d, const SparseMatrix &a)
ComplexMatrix trans_mul(const SparseComplexMatrix &m, const ComplexMatrix &a)
ComplexMatrix herm_mul(const SparseComplexMatrix &m, const ComplexMatrix &a)
base_det< Complex > ComplexDET
#define SPARSE_SSM_BOOL_OPS(S, M)
#define SPARSE_ANY_OP(DIM)
#define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, INIT_VAL, MT_RESULT)
#define SPARSE_FULL_MUL(RET_TYPE, EL_TYPE)
#define FULL_SPARSE_MUL_TRANS(RET_TYPE, EL_TYPE, CONJ_OP)
#define SPARSE_SSM_CMP_OPS(S, M)
#define SPARSE_SMSM_CMP_OPS(M1, M2)
#define SPARSE_REDUCTION_OP(RET_TYPE, EL_TYPE, OP, INIT_VAL, MT_RESULT)
#define SPARSE_SMS_CMP_OPS(M, S)
#define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE)
#define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN)
#define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN)
#define SPARSE_FULL_TRANS_MUL(RET_TYPE, EL_TYPE, CONJ_OP)
#define FULL_SPARSE_MUL(RET_TYPE, EL_TYPE)
#define SPARSE_SMSM_BOOL_OPS(M1, M2)
#define SPARSE_SMS_BOOL_OPS(M, S)
#define SPARSE_ALL_OP(DIM)
SM octinternal_do_mul_pm_sm(const PermMatrix &p, const SM &a)
SM octinternal_do_mul_sm_pm(const SM &a, const PermMatrix &p)
void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
T & xelem(octave_idx_type n)
Size of the specified dimension.
T & elem(octave_idx_type n)
Size of the specified dimension.
const T * data(void) const
Size of the specified dimension.
octave_idx_type cols(void) const
octave_idx_type rows(void) const
const T * fortran_vec(void) const
Size of the specified dimension.
ComplexColumnVector column(octave_idx_type i) const
octave_idx_type length(void) const
octave_idx_type cols(void) const
MSparse< T > permute(const Array< octave_idx_type > &vec, bool inv=false) const
MSparse< T > squeeze(void) const
MSparse< T > ipermute(const Array< octave_idx_type > &vec) const
MSparse< T > reshape(const dim_vector &new_dims) const
MSparse< T > & insert(const Sparse< T > &a, octave_idx_type r, octave_idx_type c)
MSparse< T > diag(octave_idx_type k=0) const
bool ishermitian(void) const
MatrixType transpose(void) const
void mark_as_rectangular(void)
bool is_dense(void) const
void mark_as_unsymmetric(void)
octave_idx_type * triangular_perm(void) const
int type(bool quiet=true)
Matrix sum(int dim=-1) const
RowVector row(octave_idx_type i) const
SparseComplexMatrix min(int dim=-1) const
ComplexRowVector row(octave_idx_type i) const
bool any_element_is_inf_or_nan(void) const
void * factorize(octave_idx_type &err, double &rcond, Matrix &Control, Matrix &Info, solve_singularity_handler sing_handler, bool calc_cond) const
bool any_element_is_nan(void) const
friend SparseComplexMatrix conj(const SparseComplexMatrix &a)
ComplexMatrix dsolve(MatrixType &mattype, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
SparseComplexMatrix sumsq(int dim=-1) const
SparseComplexMatrix diag(octave_idx_type k=0) const
SparseComplexMatrix dinverse(MatrixType &mattype, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
ComplexColumnVector column(octave_idx_type i) const
SparseBoolMatrix operator!(void) const
bool ishermitian(void) const
SparseComplexMatrix transpose(void) const
ComplexMatrix bsolve(MatrixType &mattype, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
SparseComplexMatrix(void)
SparseComplexMatrix & insert(const SparseComplexMatrix &a, octave_idx_type r, octave_idx_type c)
SparseComplexMatrix max(int dim=-1) const
bool all_elements_are_real(void) const
SparseComplexMatrix squeeze(void) const
bool all_integers(double &max_val, double &min_val) const
ComplexMatrix utsolve(MatrixType &mattype, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) 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
ComplexMatrix fsolve(MatrixType &mattype, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
SparseComplexMatrix inverse(void) const
SparseComplexMatrix hermitian(void) const
SparseComplexMatrix tinverse(MatrixType &mattype, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
SparseBoolMatrix any(int dim=-1) const
SparseComplexMatrix prod(int dim=-1) const
ComplexMatrix ltsolve(MatrixType &mattype, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
ComplexMatrix trisolve(MatrixType &mattype, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
SparseComplexMatrix sum(int dim=-1) const
ComplexDET determinant(void) const
SparseMatrix abs(void) 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)
ComplexMatrix solve(MatrixType &mattype, const Matrix &b) const
SparseComplexMatrix cumsum(int dim=-1) const
SparseComplexMatrix cumprod(int dim=-1) const
bool too_large_for_float(void) const
SparseBoolMatrix all(int dim=-1) const
ComplexMatrix matrix_value(void) const
SparseMatrix transpose(void) const
Array< T > array_value(void) const
octave_idx_type cols(void) const
octave_idx_type * xridx(void)
bool test_any(F fcn) const
octave_idx_type * cidx(void)
octave_idx_type nnz(void) const
Actual number of nonzero terms.
octave_idx_type rows(void) const
octave_idx_type columns(void) const
dim_vector dims(void) const
octave_idx_type * xcidx(void)
T & elem(octave_idx_type n)
octave_idx_type * ridx(void)
Vector representing the dimensions (size) of an Array.
int first_non_singleton(int def=0) const
octave_idx_type ndims(void) const
Number of dimensions.
SparseMatrix Q(void) const
SparseMatrix Pr(void) const
SparseMatrix Pc(void) const
static double get_key(const std::string &key)
ColumnVector real(const ComplexColumnVector &a)
#define F77_DBLE_CMPLX_ARG(x)
#define F77_XFCN(f, F, args)
octave_f77_int_type F77_INT
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Q
F77_RET_T const F77_INT F77_CMPLX * A
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
bool xtoo_large_for_float(double x)
void octave_write_complex(std::ostream &os, const Complex &c)
bool mx_inline_all_real(size_t n, const std::complex< T > *x)
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
void err_nan_to_logical_conversion(void)
void warn_singular_matrix(double rcond)
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
std::complex< double > Complex
#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
octave_value::octave_value(const Array< char > &chm, char type) return retval
#define BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
#define END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
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,...)