25 #if defined (HAVE_CONFIG_H)
70 #if ! defined (USE_QRSOLVE)
122 if (nr != nr_a || nc != nc_a || nz != nz_a)
139 return !(*
this ==
a);
148 if (nr == nc && nr > 0)
190 return max (dummy_idx, dim);
201 if (dim >= dv.
ndims ())
214 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
225 if (
ridx (
i) != idx_j)
262 result.
xcidx (0) = 0;
269 result.
xridx (ii++) = 0;
271 result.
xcidx (j+1) = ii;
278 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
306 idx_arg.
elem (ir) = j;
312 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
318 result.
xcidx (0) = 0;
319 result.
xcidx (1) = nel;
322 if (idx_arg(j) == -1)
326 result.
xridx (ii++) = j;
334 result.
xridx (ii++) = j;
347 return min (dummy_idx, dim);
358 if (dim >= dv.
ndims ())
371 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
382 if (
ridx (
i) != idx_j)
419 result.
xcidx (0) = 0;
426 result.
xridx (ii++) = 0;
428 result.
xcidx (j+1) = ii;
435 if (nr == 0 || nc == 0 || dim >= dv.
ndims ())
463 idx_arg.
elem (ir) = j;
469 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
475 result.
xcidx (0) = 0;
476 result.
xcidx (1) = nel;
479 if (idx_arg(j) == -1)
483 result.
xridx (ii++) = j;
491 result.
xridx (ii++) = j;
553 return insert (tmp , r, c);
569 return insert (tmp , indx);
585 if (rb.
rows () > 0 && rb.
cols () > 0)
595 if (rb.
rows () > 0 && rb.
cols () > 0)
630 retval.
xridx (q) = j;
633 assert (
nnz () == retval.
xcidx (nr));
666 return inverse (mattype, info, rcond, 0, 0);
674 return inverse (mattype, info, rcond, 0, 0);
681 return inverse (mattype, info, rcond, 0, 0);
686 double& rcond,
const bool,
687 const bool calccond)
const
695 if (nr == 0 || nc == 0 || nr != nc)
696 (*current_liboctave_error_handler) (
"inverse requires square matrix");
699 int typ = mattyp.
type ();
703 (*current_liboctave_error_handler) (
"incorrect matrix type");
736 double& rcond,
const bool,
737 const bool calccond)
const
745 if (nr == 0 || nc == 0 || nr != nc)
746 (*current_liboctave_error_handler) (
"inverse requires square matrix");
749 int typ = mattyp.
type ();
754 (*current_liboctave_error_handler) (
"incorrect matrix type");
757 double ainvnorm = 0.;
793 retval.
xdata (cx) = 1.0;
806 (*current_liboctave_error_handler) (
"division by zero");
811 rpX = retval.
xridx (colXp);
820 v -= retval.
xdata (colXp) *
data (colUp);
825 while (rpX < j && rpU < j && colXp < cx && colUp < nz);
829 colUp =
cidx (j+1) - 1;
833 if (pivot == 0. ||
ridx (colUp) != j)
834 (*current_liboctave_error_handler) (
"division by zero");
844 retval.
xridx (cx) = j;
845 retval.
xdata (cx) = v / pivot;
853 colUp =
cidx (
i+1) - 1;
857 if (pivot == 0. ||
ridx (colUp) !=
i)
858 (*current_liboctave_error_handler) (
"division by zero");
862 retval.
xdata (j) /= pivot;
864 retval.
xcidx (nr) = cx;
918 pivot =
data (cidx (jidx+1) - 1);
920 pivot =
data (cidx (jidx));
922 (*current_liboctave_error_handler) (
"division by zero");
930 colUp =
cidx (perm[iidx]+1) - 1;
932 colUp =
cidx (perm[iidx]);
936 (*current_liboctave_error_handler) (
"division by zero");
949 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
957 retval.
xridx (cx) = j;
958 retval.
xdata (cx++) = work[j];
962 retval.
xcidx (nr) = cx;
973 i < retval.
cidx (j+1);
i++)
979 rcond = 1. / ainvnorm / anorm;
987 double& rcond,
bool,
bool calc_cond)
const
989 int typ = mattype.
type (
false);
993 typ = mattype.
type (*
this);
996 ret =
dinverse (mattype, info, rcond,
true, calc_cond);
1010 rcond = fact.
rcond ();
1038 rcond = fact.
rcond ();
1074 #if defined (HAVE_UMFPACK)
1079 if (nr == 0 || nc == 0 || nr != nc)
1088 Matrix Control (UMFPACK_CONTROL, 1);
1094 Control (UMFPACK_PRL) =
tmp;
1099 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
1100 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
1106 Control (UMFPACK_FIXQ) =
tmp;
1109 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1118 reinterpret_cast<const double *
> (Ax),
1122 Matrix Info (1, UMFPACK_INFO);
1125 (nr, nc, Ap, Ai,
reinterpret_cast<const double *
> (Ax), 0,
1126 0, &Symbolic, control, info);
1135 (*current_liboctave_error_handler)
1136 (
"SparseComplexMatrix::determinant symbolic factorization failed");
1145 reinterpret_cast<const double *
> (Ax),
1146 0, Symbolic, &Numeric, control, info);
1149 rcond = Info (UMFPACK_RCOND);
1158 (*current_liboctave_error_handler)
1159 (
"SparseComplexMatrix::determinant numeric factorization failed");
1175 (*current_liboctave_error_handler)
1176 (
"SparseComplexMatrix::determinant error calculating determinant");
1188 octave_unused_parameter (err);
1189 octave_unused_parameter (rcond);
1191 (*current_liboctave_error_handler)
1192 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
1202 solve_singularity_handler,
bool calc_cond)
const
1211 if (nr != b.
rows ())
1213 (
"matrix dimension mismatch solution of linear equations");
1215 if (nr == 0 || nc == 0 || b.
cols () == 0)
1220 int typ = mattype.
type ();
1224 (*current_liboctave_error_handler) (
"incorrect matrix type");
1249 rcond = dmin / dmax;
1261 solve_singularity_handler,
1262 bool calc_cond)
const
1271 if (nr != b.
rows ())
1273 (
"matrix dimension mismatch solution of linear equations");
1275 if (nr == 0 || nc == 0 || b.
cols () == 0)
1280 int typ = mattype.
type ();
1284 (*current_liboctave_error_handler) (
"incorrect matrix type");
1290 retval.
xcidx (0) = 0;
1302 retval.
xcidx (j+1) = ii;
1312 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1320 retval.
xridx (ii) = l;
1324 retval.
xcidx (j+1) = ii;
1339 rcond = dmin / dmax;
1351 solve_singularity_handler,
1352 bool calc_cond)
const
1361 if (nr != b.
rows ())
1363 (
"matrix dimension mismatch solution of linear equations");
1365 if (nr == 0 || nc == 0 || b.
cols () == 0)
1370 int typ = mattype.
type ();
1374 (*current_liboctave_error_handler) (
"incorrect matrix type");
1399 rcond = dmin / dmax;
1411 solve_singularity_handler,
1412 bool calc_cond)
const
1421 if (nr != b.
rows ())
1423 (
"matrix dimension mismatch solution of linear equations");
1425 if (nr == 0 || nc == 0 || b.
cols () == 0)
1430 int typ = mattype.
type ();
1434 (*current_liboctave_error_handler) (
"incorrect matrix type");
1440 retval.
xcidx (0) = 0;
1452 retval.
xcidx (j+1) = ii;
1462 for (k = b.
cidx (j); k < b.
cidx (j+1); k++)
1470 retval.
xridx (ii) = l;
1474 retval.
xcidx (j+1) = ii;
1489 rcond = dmin / dmax;
1501 solve_singularity_handler sing_handler,
1502 bool calc_cond)
const
1511 if (nr != b.
rows ())
1513 (
"matrix dimension mismatch solution of linear equations");
1515 if (nr == 0 || nc == 0 || b.
cols () == 0)
1520 int typ = mattype.
type ();
1524 (*current_liboctave_error_handler) (
"incorrect matrix type");
1527 double ainvnorm = 0.;
1546 retval.
resize (nc, b_nc);
1567 goto triangular_error;
1573 i <
cidx (kidx+1)-1;
i++)
1576 work[iidx] = work[iidx] - tmp *
data (
i);
1604 i <
cidx (iidx+1)-1;
i++)
1607 work[idx2] = work[idx2] - tmp *
data (
i);
1617 if (atmp > ainvnorm)
1620 rcond = 1. / ainvnorm / anorm;
1626 retval.
resize (nc, b_nc);
1643 goto triangular_error;
1651 work[iidx] = work[iidx] - tmp *
data (
i);
1657 retval.
xelem (
i, j) = work[
i];
1680 work[iidx] = work[iidx] - tmp *
data (
i);
1690 if (atmp > ainvnorm)
1693 rcond = 1. / ainvnorm / anorm;
1702 sing_handler (rcond);
1709 volatile double rcond_plus_one = rcond + 1.0;
1717 sing_handler (rcond);
1731 solve_singularity_handler sing_handler,
1732 bool calc_cond)
const
1741 if (nr != b.
rows ())
1743 (
"matrix dimension mismatch solution of linear equations");
1745 if (nr == 0 || nc == 0 || b.
cols () == 0)
1750 int typ = mattype.
type ();
1754 (*current_liboctave_error_handler) (
"incorrect matrix type");
1757 double ainvnorm = 0.;
1776 retval.
xcidx (0) = 0;
1806 goto triangular_error;
1812 i <
cidx (kidx+1)-1;
i++)
1815 work[iidx] = work[iidx] - tmp *
data (
i);
1827 if (ii + new_nnz > x_nz)
1836 if (work[rperm[
i]] != 0.)
1839 retval.
xdata (ii++) = work[rperm[
i]];
1841 retval.
xcidx (j+1) = ii;
1865 i <
cidx (iidx+1)-1;
i++)
1868 work[idx2] = work[idx2] - tmp *
data (
i);
1878 if (atmp > ainvnorm)
1881 rcond = 1. / ainvnorm / anorm;
1903 goto triangular_error;
1911 work[iidx] = work[iidx] - tmp *
data (
i);
1923 if (ii + new_nnz > x_nz)
1935 retval.
xdata (ii++) = work[
i];
1937 retval.
xcidx (j+1) = ii;
1962 work[iidx] = work[iidx] - tmp *
data (
i);
1972 if (atmp > ainvnorm)
1975 rcond = 1. / ainvnorm / anorm;
1984 sing_handler (rcond);
1991 volatile double rcond_plus_one = rcond + 1.0;
1999 sing_handler (rcond);
2012 solve_singularity_handler sing_handler,
2013 bool calc_cond)
const
2022 if (nr != b.
rows ())
2024 (
"matrix dimension mismatch solution of linear equations");
2026 if (nr == 0 || nc == 0 || b.
cols () == 0)
2031 int typ = mattype.
type ();
2035 (*current_liboctave_error_handler) (
"incorrect matrix type");
2038 double ainvnorm = 0.;
2057 retval.
resize (nc, b_nc);
2078 goto triangular_error;
2084 i <
cidx (kidx+1)-1;
i++)
2087 work[iidx] = work[iidx] - tmp *
data (
i);
2115 i <
cidx (iidx+1)-1;
i++)
2118 work[idx2] = work[idx2] - tmp *
data (
i);
2128 if (atmp > ainvnorm)
2131 rcond = 1. / ainvnorm / anorm;
2137 retval.
resize (nc, b_nc);
2154 goto triangular_error;
2162 work[iidx] = work[iidx] - tmp *
data (
i);
2168 retval.
xelem (
i, j) = work[
i];
2191 work[iidx] = work[iidx] - tmp *
data (
i);
2201 if (atmp > ainvnorm)
2204 rcond = 1. / ainvnorm / anorm;
2213 sing_handler (rcond);
2220 volatile double rcond_plus_one = rcond + 1.0;
2228 sing_handler (rcond);
2242 solve_singularity_handler sing_handler,
2243 bool calc_cond)
const
2252 if (nr != b.
rows ())
2254 (
"matrix dimension mismatch solution of linear equations");
2256 if (nr == 0 || nc == 0 || b.
cols () == 0)
2261 int typ = mattype.
type ();
2265 (*current_liboctave_error_handler) (
"incorrect matrix type");
2268 double ainvnorm = 0.;
2287 retval.
xcidx (0) = 0;
2317 goto triangular_error;
2323 i <
cidx (kidx+1)-1;
i++)
2326 work[iidx] = work[iidx] - tmp *
data (
i);
2338 if (ii + new_nnz > x_nz)
2347 if (work[rperm[
i]] != 0.)
2350 retval.
xdata (ii++) = work[rperm[
i]];
2352 retval.
xcidx (j+1) = ii;
2376 i <
cidx (iidx+1)-1;
i++)
2379 work[idx2] = work[idx2] - tmp *
data (
i);
2389 if (atmp > ainvnorm)
2392 rcond = 1. / ainvnorm / anorm;
2414 goto triangular_error;
2422 work[iidx] = work[iidx] - tmp *
data (
i);
2434 if (ii + new_nnz > x_nz)
2446 retval.
xdata (ii++) = work[
i];
2448 retval.
xcidx (j+1) = ii;
2473 work[iidx] = work[iidx] - tmp *
data (
i);
2483 if (atmp > ainvnorm)
2486 rcond = 1. / ainvnorm / anorm;
2495 sing_handler (rcond);
2502 volatile double rcond_plus_one = rcond + 1.0;
2510 sing_handler (rcond);
2524 solve_singularity_handler sing_handler,
2525 bool calc_cond)
const
2534 if (nr != b.
rows ())
2536 (
"matrix dimension mismatch solution of linear equations");
2538 if (nr == 0 || nc == 0 || b.
cols () == 0)
2543 int typ = mattype.
type ();
2547 (*current_liboctave_error_handler) (
"incorrect matrix type");
2550 double ainvnorm = 0.;
2569 retval.
resize (nc, b_nc);
2578 work[perm[
i]] =
b(
i,j);
2588 if (perm[
ridx (
i)] < minr)
2590 minr = perm[
ridx (
i)];
2594 if (minr !=
k ||
data (mini) == 0.)
2597 goto triangular_error;
2608 work[iidx] = work[iidx] - tmp *
data (
i);
2636 if (perm[
ridx (
i)] < minr)
2638 minr = perm[
ridx (
i)];
2651 work[iidx] = work[iidx] - tmp *
data (
i);
2662 if (atmp > ainvnorm)
2665 rcond = 1. / ainvnorm / anorm;
2671 retval.
resize (nc, b_nc, 0.);
2686 goto triangular_error;
2694 work[iidx] = work[iidx] - tmp *
data (
i);
2699 retval.
xelem (
i, j) = work[
i];
2723 work[iidx] = work[iidx] - tmp *
data (
i);
2733 if (atmp > ainvnorm)
2736 rcond = 1. / ainvnorm / anorm;
2744 sing_handler (rcond);
2751 volatile double rcond_plus_one = rcond + 1.0;
2759 sing_handler (rcond);
2773 solve_singularity_handler sing_handler,
2774 bool calc_cond)
const
2784 if (nr != b.
rows ())
2786 (
"matrix dimension mismatch solution of linear equations");
2788 if (nr == 0 || nc == 0 || b.
cols () == 0)
2793 int typ = mattype.
type ();
2797 (*current_liboctave_error_handler) (
"incorrect matrix type");
2800 double ainvnorm = 0.;
2819 retval.
xcidx (0) = 0;
2843 if (perm[
ridx (
i)] < minr)
2845 minr = perm[
ridx (
i)];
2849 if (minr !=
k ||
data (mini) == 0.)
2852 goto triangular_error;
2863 work[iidx] = work[iidx] - tmp *
data (
i);
2875 if (ii + new_nnz > x_nz)
2887 retval.
xdata (ii++) = work[
i];
2889 retval.
xcidx (j+1) = ii;
2913 if (perm[
ridx (
i)] < minr)
2915 minr = perm[
ridx (
i)];
2928 work[iidx] = work[iidx] - tmp *
data (
i);
2939 if (atmp > ainvnorm)
2942 rcond = 1. / ainvnorm / anorm;
2963 goto triangular_error;
2971 work[iidx] = work[iidx] - tmp *
data (
i);
2983 if (ii + new_nnz > x_nz)
2995 retval.
xdata (ii++) = work[
i];
2997 retval.
xcidx (j+1) = ii;
3023 work[iidx] = work[iidx] - tmp *
data (
i);
3033 if (atmp > ainvnorm)
3036 rcond = 1. / ainvnorm / anorm;
3045 sing_handler (rcond);
3052 volatile double rcond_plus_one = rcond + 1.0;
3060 sing_handler (rcond);
3074 solve_singularity_handler sing_handler,
3075 bool calc_cond)
const
3084 if (nr != b.
rows ())
3086 (
"matrix dimension mismatch solution of linear equations");
3088 if (nr == 0 || nc == 0 || b.
cols () == 0)
3093 int typ = mattype.
type ();
3097 (*current_liboctave_error_handler) (
"incorrect matrix type");
3100 double ainvnorm = 0.;
3119 retval.
resize (nc, b_nc);
3128 work[perm[
i]] =
b(
i,j);
3138 if (perm[
ridx (
i)] < minr)
3140 minr = perm[
ridx (
i)];
3144 if (minr !=
k ||
data (mini) == 0.)
3147 goto triangular_error;
3158 work[iidx] = work[iidx] - tmp *
data (
i);
3186 if (perm[
ridx (
i)] < minr)
3188 minr = perm[
ridx (
i)];
3201 work[iidx] = work[iidx] - tmp *
data (
i);
3212 if (atmp > ainvnorm)
3215 rcond = 1. / ainvnorm / anorm;
3221 retval.
resize (nc, b_nc, 0.);
3237 goto triangular_error;
3245 work[iidx] = work[iidx] - tmp *
data (
i);
3251 retval.
xelem (
i, j) = work[
i];
3275 work[iidx] = work[iidx] - tmp *
data (
i);
3285 if (atmp > ainvnorm)
3288 rcond = 1. / ainvnorm / anorm;
3297 sing_handler (rcond);
3304 volatile double rcond_plus_one = rcond + 1.0;
3312 sing_handler (rcond);
3326 solve_singularity_handler sing_handler,
3327 bool calc_cond)
const
3336 if (nr != b.
rows ())
3338 (
"matrix dimension mismatch solution of linear equations");
3340 if (nr == 0 || nc == 0 || b.
cols () == 0)
3345 int typ = mattype.
type ();
3349 (*current_liboctave_error_handler) (
"incorrect matrix type");
3352 double ainvnorm = 0.;
3371 retval.
xcidx (0) = 0;
3395 if (perm[
ridx (
i)] < minr)
3397 minr = perm[
ridx (
i)];
3401 if (minr !=
k ||
data (mini) == 0.)
3404 goto triangular_error;
3415 work[iidx] = work[iidx] - tmp *
data (
i);
3427 if (ii + new_nnz > x_nz)
3439 retval.
xdata (ii++) = work[
i];
3441 retval.
xcidx (j+1) = ii;
3465 if (perm[
ridx (
i)] < minr)
3467 minr = perm[
ridx (
i)];
3480 work[iidx] = work[iidx] - tmp *
data (
i);
3491 if (atmp > ainvnorm)
3494 rcond = 1. / ainvnorm / anorm;
3515 goto triangular_error;
3523 work[iidx] = work[iidx] - tmp *
data (
i);
3535 if (ii + new_nnz > x_nz)
3547 retval.
xdata (ii++) = work[
i];
3549 retval.
xcidx (j+1) = ii;
3575 work[iidx] = work[iidx] - tmp *
data (
i);
3585 if (atmp > ainvnorm)
3588 rcond = 1. / ainvnorm / anorm;
3597 sing_handler (rcond);
3604 volatile double rcond_plus_one = rcond + 1.0;
3612 sing_handler (rcond);
3626 solve_singularity_handler sing_handler,
3627 bool calc_cond)
const
3635 if (nr != nc || nr != b.
rows ())
3637 (
"matrix dimension mismatch solution of linear equations");
3639 if (nr == 0 || b.
cols () == 0)
3642 (*current_liboctave_error_handler)
3643 (
"calculation of condition number not implemented");
3647 volatile int typ = mattype.
type ();
3681 else if (
ridx (
i) == j + 1)
3717 DL[j] =
data (ii++);
3718 DU[j] =
data (ii++);
3720 D[nc-1] =
data (ii);
3737 else if (
ridx (
i) == j + 1)
3739 else if (
ridx (
i) == j - 1)
3759 sing_handler (rcond);
3770 (*current_liboctave_error_handler) (
"incorrect matrix type");
3779 solve_singularity_handler sing_handler,
3780 bool calc_cond)
const
3788 if (nr != nc || nr != b.
rows ())
3790 (
"matrix dimension mismatch solution of linear equations");
3792 if (nr == 0 || b.
cols () == 0)
3795 (*current_liboctave_error_handler)
3796 (
"calculation of condition number not implemented");
3800 int typ = mattype.
type ();
3821 DL[j] =
data (ii++);
3822 DU[j] =
data (ii++);
3824 D[nc-1] =
data (ii);
3841 else if (
ridx (
i) == j + 1)
3843 else if (
ridx (
i) == j - 1)
3858 sing_handler (rcond);
3870 retval.
xcidx (0) = 0;
3884 (F77_CONST_CHAR_ARG2 (&job, 1),
3888 F77_CHAR_ARG_LEN (1)));
3897 if (ii + new_nnz > x_nz)
3909 retval.
xdata (ii++) = work[
i];
3911 retval.
xcidx (j+1) = ii;
3918 (*current_liboctave_error_handler) (
"incorrect matrix type");
3927 solve_singularity_handler sing_handler,
3928 bool calc_cond)
const
3936 if (nr != nc || nr != b.
rows ())
3938 (
"matrix dimension mismatch solution of linear equations");
3940 if (nr == 0 || b.
cols () == 0)
3943 (*current_liboctave_error_handler)
3944 (
"calculation of condition number not implemented");
3948 volatile int typ = mattype.
type ();
3982 else if (
ridx (
i) == j + 1)
4019 DL[j] =
data (ii++);
4020 DU[j] =
data (ii++);
4022 D[nc-1] =
data (ii);
4039 else if (
ridx (
i) == j + 1)
4041 else if (
ridx (
i) == j - 1)
4064 sing_handler (rcond);
4072 (*current_liboctave_error_handler) (
"incorrect matrix type");
4082 solve_singularity_handler sing_handler,
4083 bool calc_cond)
const
4091 if (nr != nc || nr != b.
rows ())
4093 (
"matrix dimension mismatch solution of linear equations");
4095 if (nr == 0 || b.
cols () == 0)
4098 (*current_liboctave_error_handler)
4099 (
"calculation of condition number not implemented");
4103 int typ = mattype.
type ();
4124 DL[j] =
data (ii++);
4125 DU[j] =
data (ii++);
4127 D[nc-1] =
data (ii);
4144 else if (
ridx (
i) == j + 1)
4146 else if (
ridx (
i) == j - 1)
4161 sing_handler (rcond);
4181 retval.
xcidx (0) = 0;
4189 (F77_CONST_CHAR_ARG2 (&job, 1),
4193 F77_CHAR_ARG_LEN (1)));
4199 (*current_liboctave_error_handler)
4200 (
"SparseComplexMatrix::solve solve failed");
4213 if (ii + new_nnz > x_nz)
4225 retval.
xdata (ii++) = Bx[
i];
4228 retval.
xcidx (j+1) = ii;
4235 (*current_liboctave_error_handler) (
"incorrect matrix type");
4244 solve_singularity_handler sing_handler,
4245 bool calc_cond)
const
4253 if (nr != nc || nr != b.
rows ())
4255 (
"matrix dimension mismatch solution of linear equations");
4257 if (nr == 0 || b.
cols () == 0)
4262 volatile int typ = mattype.
type ();
4278 tmp_data[ii++] = 0.;
4286 m_band(ri - j, j) =
data (
i);
4295 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4297 F77_CHAR_ARG_LEN (1)));
4318 (F77_CONST_CHAR_ARG2 (&job, 1),
4321 F77_CHAR_ARG_LEN (1)));
4326 volatile double rcond_plus_one = rcond + 1.0;
4334 sing_handler (rcond);
4352 (F77_CONST_CHAR_ARG2 (&job, 1),
4355 F77_CHAR_ARG_LEN (1)));
4360 (*current_liboctave_error_handler)
4361 (
"SparseMatrix::solve solve failed");
4384 tmp_data[ii++] = 0.;
4389 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
4408 F77_XFCN (zgbtrf, ZGBTRF, (nr, nc, n_lower, n_upper,
4421 sing_handler (rcond);
4438 (F77_CONST_CHAR_ARG2 (&job, 1),
4441 F77_CHAR_ARG_LEN (1)));
4446 volatile double rcond_plus_one = rcond + 1.0;
4454 sing_handler (rcond);
4473 (F77_CONST_CHAR_ARG2 (&job, 1),
4476 F77_CHAR_ARG_LEN (1)));
4481 (*current_liboctave_error_handler) (
"incorrect matrix type");
4490 solve_singularity_handler sing_handler,
4491 bool calc_cond)
const
4499 if (nr != nc || nr != b.
rows ())
4501 (
"matrix dimension mismatch solution of linear equations");
4503 if (nr == 0 || b.
cols () == 0)
4508 volatile int typ = mattype.
type ();
4525 tmp_data[ii++] = 0.;
4533 m_band(ri - j, j) =
data (
i);
4542 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4544 F77_CHAR_ARG_LEN (1)));
4563 (F77_CONST_CHAR_ARG2 (&job, 1),
4566 F77_CHAR_ARG_LEN (1)));
4571 volatile double rcond_plus_one = rcond + 1.0;
4579 sing_handler (rcond);
4601 retval.
xcidx (0) = 0;
4608 (F77_CONST_CHAR_ARG2 (&job, 1),
4611 F77_CHAR_ARG_LEN (1)));
4616 (*current_liboctave_error_handler)
4617 (
"SparseComplexMatrix::solve solve failed");
4632 sz = (sz > 10 ? sz : 10) + x_nz;
4640 retval.
xcidx (j+1) = ii;
4664 tmp_data[ii++] = 0.;
4669 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
4688 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper,
4699 sing_handler (rcond);
4716 (F77_CONST_CHAR_ARG2 (&job, 1),
4719 F77_CHAR_ARG_LEN (1)));
4724 volatile double rcond_plus_one = rcond + 1.0;
4732 sing_handler (rcond);
4748 retval.
xcidx (0) = 0;
4762 (F77_CONST_CHAR_ARG2 (&job, 1),
4765 F77_CHAR_ARG_LEN (1)));
4774 if (ii + new_nnz > x_nz)
4786 retval.
xdata (ii++) = work[
i];
4788 retval.
xcidx (j+1) = ii;
4796 (*current_liboctave_error_handler) (
"incorrect matrix type");
4805 solve_singularity_handler sing_handler,
4806 bool calc_cond)
const
4814 if (nr != nc || nr != b.
rows ())
4816 (
"matrix dimension mismatch solution of linear equations");
4818 if (nr == 0 || b.
cols () == 0)
4823 volatile int typ = mattype.
type ();
4840 tmp_data[ii++] = 0.;
4848 m_band(ri - j, j) =
data (
i);
4857 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4859 F77_CHAR_ARG_LEN (1)));
4880 (F77_CONST_CHAR_ARG2 (&job, 1),
4883 F77_CHAR_ARG_LEN (1)));
4888 volatile double rcond_plus_one = rcond + 1.0;
4896 sing_handler (rcond);
4914 (F77_CONST_CHAR_ARG2 (&job, 1),
4917 F77_CHAR_ARG_LEN (1)));
4922 (*current_liboctave_error_handler)
4923 (
"SparseComplexMatrix::solve solve failed");
4946 tmp_data[ii++] = 0.;
4951 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
4970 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper,
4981 sing_handler (rcond);
4998 (F77_CONST_CHAR_ARG2 (&job, 1),
5001 F77_CHAR_ARG_LEN (1)));
5006 volatile double rcond_plus_one = rcond + 1.0;
5014 sing_handler (rcond);
5032 (F77_CONST_CHAR_ARG2 (&job, 1),
5035 F77_CHAR_ARG_LEN (1)));
5040 (*current_liboctave_error_handler) (
"incorrect matrix type");
5049 solve_singularity_handler sing_handler,
5050 bool calc_cond)
const
5058 if (nr != nc || nr != b.
rows ())
5060 (
"matrix dimension mismatch solution of linear equations");
5062 if (nr == 0 || b.
cols () == 0)
5067 volatile int typ = mattype.
type ();
5084 tmp_data[ii++] = 0.;
5092 m_band(ri - j, j) =
data (
i);
5101 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5103 F77_CHAR_ARG_LEN (1)));
5125 (F77_CONST_CHAR_ARG2 (&job, 1),
5128 F77_CHAR_ARG_LEN (1)));
5133 volatile double rcond_plus_one = rcond + 1.0;
5141 sing_handler (rcond);
5163 retval.
xcidx (0) = 0;
5171 (F77_CONST_CHAR_ARG2 (&job, 1),
5174 F77_CHAR_ARG_LEN (1)));
5179 (*current_liboctave_error_handler)
5180 (
"SparseMatrix::solve solve failed");
5192 if (ii + new_nnz > x_nz)
5204 retval.
xdata (ii++) = Bx[
i];
5207 retval.
xcidx (j+1) = ii;
5231 tmp_data[ii++] = 0.;
5236 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
5255 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper,
5266 sing_handler (rcond);
5283 (F77_CONST_CHAR_ARG2 (&job, 1),
5286 F77_CHAR_ARG_LEN (1)));
5291 volatile double rcond_plus_one = rcond + 1.0;
5299 sing_handler (rcond);
5315 retval.
xcidx (0) = 0;
5330 (F77_CONST_CHAR_ARG2 (&job, 1),
5333 F77_CHAR_ARG_LEN (1)));
5342 if (ii + new_nnz > x_nz)
5354 retval.
xdata (ii++) = Bx[
i];
5356 retval.
xcidx (j+1) = ii;
5364 (*current_liboctave_error_handler) (
"incorrect matrix type");
5373 solve_singularity_handler sing_handler,
5374 bool calc_cond)
const
5380 #if defined (HAVE_UMFPACK)
5383 Control =
Matrix (UMFPACK_CONTROL, 1);
5389 Control (UMFPACK_PRL) =
tmp;
5393 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
5394 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
5400 Control (UMFPACK_FIXQ) =
tmp;
5411 reinterpret_cast<const double *
> (Ax),
5415 Info =
Matrix (1, UMFPACK_INFO);
5418 reinterpret_cast<const double *
> (Ax),
5419 0, 0, &Symbolic, control, info);
5429 (*current_liboctave_error_handler)
5430 (
"SparseComplexMatrix::solve symbolic factorization failed");
5438 reinterpret_cast<const double *
> (Ax),
5439 0, Symbolic, &Numeric, control, info);
5443 rcond = Info (UMFPACK_RCOND);
5446 volatile double rcond_plus_one = rcond + 1.0;
5448 if (status == UMFPACK_WARNING_singular_matrix
5456 sing_handler (rcond);
5460 else if (status < 0)
5466 (*current_liboctave_error_handler)
5467 (
"SparseComplexMatrix::solve numeric factorization failed");
5482 octave_unused_parameter (rcond);
5483 octave_unused_parameter (Control);
5484 octave_unused_parameter (Info);
5485 octave_unused_parameter (sing_handler);
5486 octave_unused_parameter (calc_cond);
5488 (*current_liboctave_error_handler)
5489 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
5499 solve_singularity_handler sing_handler,
5500 bool calc_cond)
const
5508 if (nr != nc || nr != b.
rows ())
5510 (
"matrix dimension mismatch solution of linear equations");
5512 if (nr == 0 || b.
cols () == 0)
5517 volatile int typ = mattype.
type ();
5522 #if defined (HAVE_CHOLMOD)
5523 cholmod_common Common;
5524 cholmod_common *cm = &Common;
5527 CHOLMOD_NAME(
start) (cm);
5528 cm->prefer_zomplex =
false;
5534 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
5538 cm->print =
static_cast<int> (spu) + 2;
5539 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5543 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5544 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5546 cm->final_ll =
true;
5548 cholmod_sparse Astore;
5549 cholmod_sparse *
A = &Astore;
5560 #if defined (OCTAVE_ENABLE_64)
5561 A->itype = CHOLMOD_LONG;
5563 A->itype = CHOLMOD_INT;
5565 A->dtype = CHOLMOD_DOUBLE;
5567 A->xtype = CHOLMOD_COMPLEX;
5574 cholmod_dense Bstore;
5575 cholmod_dense *
B = &Bstore;
5576 B->nrow = b.
rows ();
5577 B->ncol = b.
cols ();
5579 B->nzmax = B->nrow * B->ncol;
5580 B->dtype = CHOLMOD_DOUBLE;
5581 B->xtype = CHOLMOD_REAL;
5582 if (nc < 1 || b.
cols () < 1)
5590 L = CHOLMOD_NAME(analyze) (
A, cm);
5593 rcond = CHOLMOD_NAME(rcond)(L, cm);
5606 volatile double rcond_plus_one = rcond + 1.0;
5614 sing_handler (rcond);
5625 X = CHOLMOD_NAME(
solve) (CHOLMOD_A, L,
B, cm);
5637 CHOLMOD_NAME(free_dense) (&X, cm);
5638 CHOLMOD_NAME(free_factor) (&L, cm);
5639 CHOLMOD_NAME(finish) (cm);
5640 static char tmp[] =
" ";
5641 CHOLMOD_NAME(print_common) (
tmp, cm);
5645 (*current_liboctave_warning_with_id_handler)
5646 (
"Octave:missing-dependency",
5647 "support for CHOLMOD was unavailable or disabled "
5648 "when liboctave was built");
5657 #if defined (HAVE_UMFPACK)
5659 void *Numeric =
factorize (err, rcond, Control, Info,
5660 sing_handler, calc_cond);
5672 #if defined (UMFPACK_SEPARATE_SPLIT)
5680 retval.
resize (b_nr, b_nc);
5685 #if defined (UMFPACK_SEPARATE_SPLIT)
5688 reinterpret_cast<const double *
> (Ax),
5690 reinterpret_cast<double *> (&Xx[iidx]),
5692 &Bx[iidx], Bz, Numeric,
5700 reinterpret_cast<const double *
> (Ax),
5702 reinterpret_cast<double *> (&Xx[iidx]),
5704 reinterpret_cast<const double *
> (Bz),
5714 (*current_liboctave_error_handler)
5715 (
"SparseComplexMatrix::solve solve failed");
5730 octave_unused_parameter (rcond);
5731 octave_unused_parameter (sing_handler);
5732 octave_unused_parameter (calc_cond);
5734 (*current_liboctave_error_handler)
5735 (
"support for UMFPACK was unavailable or disabled "
5736 "when liboctave was built");
5740 (*current_liboctave_error_handler) (
"incorrect matrix type");
5749 solve_singularity_handler sing_handler,
5750 bool calc_cond)
const
5758 if (nr != nc || nr != b.
rows ())
5760 (
"matrix dimension mismatch solution of linear equations");
5762 if (nr == 0 || b.
cols () == 0)
5767 volatile int typ = mattype.
type ();
5772 #if defined (HAVE_CHOLMOD)
5773 cholmod_common Common;
5774 cholmod_common *cm = &Common;
5777 CHOLMOD_NAME(
start) (cm);
5778 cm->prefer_zomplex =
false;
5784 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
5788 cm->print =
static_cast<int> (spu) + 2;
5789 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5793 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5794 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5796 cm->final_ll =
true;
5798 cholmod_sparse Astore;
5799 cholmod_sparse *
A = &Astore;
5810 #if defined (OCTAVE_ENABLE_64)
5811 A->itype = CHOLMOD_LONG;
5813 A->itype = CHOLMOD_INT;
5815 A->dtype = CHOLMOD_DOUBLE;
5817 A->xtype = CHOLMOD_COMPLEX;
5824 cholmod_sparse Bstore;
5825 cholmod_sparse *
B = &Bstore;
5826 B->nrow = b.
rows ();
5827 B->ncol = b.
cols ();
5830 B->nzmax = b.
nnz ();
5834 #if defined (OCTAVE_ENABLE_64)
5835 B->itype = CHOLMOD_LONG;
5837 B->itype = CHOLMOD_INT;
5839 B->dtype = CHOLMOD_DOUBLE;
5841 B->xtype = CHOLMOD_REAL;
5843 if (b.
rows () < 1 || b.
cols () < 1)
5850 L = CHOLMOD_NAME(analyze) (
A, cm);
5853 rcond = CHOLMOD_NAME(rcond)(L, cm);
5866 volatile double rcond_plus_one = rcond + 1.0;
5874 sing_handler (rcond);
5885 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L,
B, cm);
5889 (static_cast<octave_idx_type>(X->nrow),
5890 static_cast<octave_idx_type>(X->ncol),
5891 static_cast<octave_idx_type>(X->nzmax));
5893 j <= static_cast<octave_idx_type>(X->ncol); j++)
5896 j < static_cast<octave_idx_type>(X->nzmax); j++)
5903 CHOLMOD_NAME(free_sparse) (&X, cm);
5904 CHOLMOD_NAME(free_factor) (&L, cm);
5905 CHOLMOD_NAME(finish) (cm);
5906 static char tmp[] =
" ";
5907 CHOLMOD_NAME(print_common) (
tmp, cm);
5911 (*current_liboctave_warning_with_id_handler)
5912 (
"Octave:missing-dependency",
5913 "support for CHOLMOD was unavailable or disabled "
5914 "when liboctave was built");
5923 #if defined (HAVE_UMFPACK)
5925 void *Numeric =
factorize (err, rcond, Control, Info,
5926 sing_handler, calc_cond);
5939 #if defined (UMFPACK_SEPARATE_SPLIT)
5956 retval.
xcidx (0) = 0;
5960 #if defined (UMFPACK_SEPARATE_SPLIT)
5966 reinterpret_cast<const double *
> (Ax),
5968 reinterpret_cast<double *> (Xx),
5970 Bx, Bz, Numeric, control,
5977 reinterpret_cast<const double *
> (Ax),
5979 reinterpret_cast<double *> (Xx),
5981 reinterpret_cast<double *
> (Bz),
5991 (*current_liboctave_error_handler)
5992 (
"SparseComplexMatrix::solve solve failed");
6007 sz = (sz > 10 ? sz : 10) + x_nz;
6015 retval.
xcidx (j+1) = ii;
6028 octave_unused_parameter (rcond);
6029 octave_unused_parameter (sing_handler);
6030 octave_unused_parameter (calc_cond);
6032 (*current_liboctave_error_handler)
6033 (
"support for UMFPACK was unavailable or disabled "
6034 "when liboctave was built");
6038 (*current_liboctave_error_handler) (
"incorrect matrix type");
6047 solve_singularity_handler sing_handler,
6048 bool calc_cond)
const
6056 if (nr != nc || nr != b.
rows ())
6058 (
"matrix dimension mismatch solution of linear equations");
6060 if (nr == 0 || b.
cols () == 0)
6065 volatile int typ = mattype.
type ();
6070 #if defined (HAVE_CHOLMOD)
6071 cholmod_common Common;
6072 cholmod_common *cm = &Common;
6075 CHOLMOD_NAME(
start) (cm);
6076 cm->prefer_zomplex =
false;
6082 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
6086 cm->print =
static_cast<int> (spu) + 2;
6087 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6091 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6092 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6094 cm->final_ll =
true;
6096 cholmod_sparse Astore;
6097 cholmod_sparse *
A = &Astore;
6108 #if defined (OCTAVE_ENABLE_64)
6109 A->itype = CHOLMOD_LONG;
6111 A->itype = CHOLMOD_INT;
6113 A->dtype = CHOLMOD_DOUBLE;
6115 A->xtype = CHOLMOD_COMPLEX;
6122 cholmod_dense Bstore;
6123 cholmod_dense *
B = &Bstore;
6124 B->nrow = b.
rows ();
6125 B->ncol = b.
cols ();
6127 B->nzmax = B->nrow * B->ncol;
6128 B->dtype = CHOLMOD_DOUBLE;
6129 B->xtype = CHOLMOD_COMPLEX;
6130 if (nc < 1 || b.
cols () < 1)
6138 L = CHOLMOD_NAME(analyze) (
A, cm);
6141 rcond = CHOLMOD_NAME(rcond)(L, cm);
6154 volatile double rcond_plus_one = rcond + 1.0;
6162 sing_handler (rcond);
6173 X = CHOLMOD_NAME(
solve) (CHOLMOD_A, L,
B, cm);
6185 CHOLMOD_NAME(free_dense) (&X, cm);
6186 CHOLMOD_NAME(free_factor) (&L, cm);
6187 CHOLMOD_NAME(finish) (cm);
6188 static char tmp[] =
" ";
6189 CHOLMOD_NAME(print_common) (
tmp, cm);
6193 (*current_liboctave_warning_with_id_handler)
6194 (
"Octave:missing-dependency",
6195 "support for CHOLMOD was unavailable or disabled "
6196 "when liboctave was built");
6205 #if defined (HAVE_UMFPACK)
6207 void *Numeric =
factorize (err, rcond, Control, Info,
6208 sing_handler, calc_cond);
6222 retval.
resize (b_nr, b_nc);
6229 reinterpret_cast<const double *
> (Ax),
6231 reinterpret_cast<double *> (&Xx[iidx]),
6233 reinterpret_cast<const double *
> (&Bx[iidx]),
6234 0, Numeric, control, info);
6241 (*current_liboctave_error_handler)
6242 (
"SparseComplexMatrix::solve solve failed");
6257 octave_unused_parameter (rcond);
6258 octave_unused_parameter (sing_handler);
6259 octave_unused_parameter (calc_cond);
6261 (*current_liboctave_error_handler)
6262 (
"support for UMFPACK was unavailable or disabled "
6263 "when liboctave was built");
6267 (*current_liboctave_error_handler) (
"incorrect matrix type");
6276 solve_singularity_handler sing_handler,
6277 bool calc_cond)
const
6285 if (nr != nc || nr != b.
rows ())
6287 (
"matrix dimension mismatch solution of linear equations");
6289 if (nr == 0 || b.
cols () == 0)
6294 volatile int typ = mattype.
type ();
6299 #if defined (HAVE_CHOLMOD)
6300 cholmod_common Common;
6301 cholmod_common *cm = &Common;
6304 CHOLMOD_NAME(
start) (cm);
6305 cm->prefer_zomplex =
false;
6311 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
6315 cm->print =
static_cast<int> (spu) + 2;
6316 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6320 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6321 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6323 cm->final_ll =
true;
6325 cholmod_sparse Astore;
6326 cholmod_sparse *
A = &Astore;
6337 #if defined (OCTAVE_ENABLE_64)
6338 A->itype = CHOLMOD_LONG;
6340 A->itype = CHOLMOD_INT;
6342 A->dtype = CHOLMOD_DOUBLE;
6344 A->xtype = CHOLMOD_COMPLEX;
6351 cholmod_sparse Bstore;
6352 cholmod_sparse *
B = &Bstore;
6353 B->nrow = b.
rows ();
6354 B->ncol = b.
cols ();
6357 B->nzmax = b.
nnz ();
6361 #if defined (OCTAVE_ENABLE_64)
6362 B->itype = CHOLMOD_LONG;
6364 B->itype = CHOLMOD_INT;
6366 B->dtype = CHOLMOD_DOUBLE;
6368 B->xtype = CHOLMOD_COMPLEX;
6370 if (b.
rows () < 1 || b.
cols () < 1)
6377 L = CHOLMOD_NAME(analyze) (
A, cm);
6380 rcond = CHOLMOD_NAME(rcond)(L, cm);
6393 volatile double rcond_plus_one = rcond + 1.0;
6401 sing_handler (rcond);
6412 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L,
B, cm);
6416 (static_cast<octave_idx_type>(X->nrow),
6417 static_cast<octave_idx_type>(X->ncol),
6418 static_cast<octave_idx_type>(X->nzmax));
6420 j <= static_cast<octave_idx_type>(X->ncol); j++)
6423 j < static_cast<octave_idx_type>(X->nzmax); j++)
6430 CHOLMOD_NAME(free_sparse) (&X, cm);
6431 CHOLMOD_NAME(free_factor) (&L, cm);
6432 CHOLMOD_NAME(finish) (cm);
6433 static char tmp[] =
" ";
6434 CHOLMOD_NAME(print_common) (
tmp, cm);
6438 (*current_liboctave_warning_with_id_handler)
6439 (
"Octave:missing-dependency",
6440 "support for CHOLMOD was unavailable or disabled "
6441 "when liboctave was built");
6450 #if defined (HAVE_UMFPACK)
6452 void *Numeric =
factorize (err, rcond, Control, Info,
6453 sing_handler, calc_cond);
6476 retval.
xcidx (0) = 0;
6484 reinterpret_cast<const double *
> (Ax),
6486 reinterpret_cast<double *> (Xx),
6488 reinterpret_cast<double *
> (Bx),
6489 0, Numeric, control, info);
6496 (*current_liboctave_error_handler)
6497 (
"SparseComplexMatrix::solve solve failed");
6512 sz = (sz > 10 ? sz : 10) + x_nz;
6520 retval.
xcidx (j+1) = ii;
6525 rcond = Info (UMFPACK_RCOND);
6526 volatile double rcond_plus_one = rcond + 1.0;
6528 if (status == UMFPACK_WARNING_singular_matrix
6534 sing_handler (rcond);
6547 octave_unused_parameter (rcond);
6548 octave_unused_parameter (sing_handler);
6549 octave_unused_parameter (calc_cond);
6551 (*current_liboctave_error_handler)
6552 (
"support for UMFPACK was unavailable or disabled "
6553 "when liboctave was built");
6557 (*current_liboctave_error_handler) (
"incorrect matrix type");
6568 return solve (mattype, b, info, rcond, 0);
6576 return solve (mattype, b, info, rcond, 0);
6583 return solve (mattype, b, info, rcond, 0);
6589 solve_singularity_handler sing_handler,
6590 bool singular_fallback)
const
6593 int typ = mattype.
type (
false);
6596 typ = mattype.
type (*
this);
6599 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6601 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6603 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6605 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6608 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6610 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6612 (*current_liboctave_error_handler) (
"unknown matrix type");
6617 #if defined (USE_QRSOLVE)
6618 retval =
qrsolve (*
this, b, err);
6633 return solve (mattype, b, info, rcond, 0);
6641 return solve (mattype, b, info, rcond, 0);
6648 return solve (mattype, b, info, rcond, 0);
6654 solve_singularity_handler sing_handler,
6655 bool singular_fallback)
const
6658 int typ = mattype.
type (
false);
6661 typ = mattype.
type (*
this);
6664 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6666 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6668 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6670 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6673 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6675 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6677 (*current_liboctave_error_handler) (
"unknown matrix type");
6682 #if defined (USE_QRSOLVE)
6683 retval =
qrsolve (*
this, b, err);
6698 return solve (mattype, b, info, rcond, 0);
6706 return solve (mattype, b, info, rcond, 0);
6713 return solve (mattype, b, info, rcond, 0);
6719 solve_singularity_handler sing_handler,
6720 bool singular_fallback)
const
6723 int typ = mattype.
type (
false);
6726 typ = mattype.
type (*
this);
6729 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6731 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6733 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6735 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6738 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6740 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6742 (*current_liboctave_error_handler) (
"unknown matrix type");
6747 #if defined (USE_QRSOLVE)
6748 retval =
qrsolve (*
this, b, err);
6764 return solve (mattype, b, info, rcond, 0);
6772 return solve (mattype, b, info, rcond, 0);
6779 return solve (mattype, b, info, rcond, 0);
6785 solve_singularity_handler sing_handler,
6786 bool singular_fallback)
const
6789 int typ = mattype.
type (
false);
6792 typ = mattype.
type (*
this);
6795 retval =
dsolve (mattype, b, err, rcond, sing_handler,
false);
6797 retval =
utsolve (mattype, b, err, rcond, sing_handler,
false);
6799 retval =
ltsolve (mattype, b, err, rcond, sing_handler,
false);
6801 retval =
bsolve (mattype, b, err, rcond, sing_handler,
false);
6804 retval =
trisolve (mattype, b, err, rcond, sing_handler,
false);
6806 retval =
fsolve (mattype, b, err, rcond, sing_handler,
true);
6808 (*current_liboctave_error_handler) (
"unknown matrix type");
6813 #if defined (USE_QRSOLVE)
6814 retval =
qrsolve (*
this, b, err);
6817 SparseComplexMatrix> (*
this,
b,
err);
6828 return solve (mattype, b, info, rcond);
6836 return solve (mattype, b, info, rcond);
6843 return solve (mattype, b, info, rcond, 0);
6849 solve_singularity_handler sing_handler)
const
6852 return solve (mattype, tmp, info, rcond,
6853 sing_handler).
column (static_cast<octave_idx_type> (0));
6862 return solve (mattype, b, info, rcond, 0);
6870 return solve (mattype, b, info, rcond, 0);
6877 return solve (mattype, b, info, rcond, 0);
6883 solve_singularity_handler sing_handler)
const
6886 return solve (mattype, tmp, info, rcond,
6887 sing_handler).
column (static_cast<octave_idx_type> (0));
6895 return solve (b, info, rcond, 0);
6902 return solve (b, info, rcond, 0);
6907 double& rcond)
const
6909 return solve (b, info, rcond, 0);
6915 solve_singularity_handler sing_handler)
const
6918 return solve (mattype, b, err, rcond, sing_handler);
6926 return solve (b, info, rcond, 0);
6934 return solve (b, info, rcond, 0);
6941 return solve (b, info, rcond, 0);
6947 solve_singularity_handler sing_handler)
const
6950 return solve (mattype, b, err, rcond, sing_handler);
6958 return solve (b, info, rcond, 0);
6965 return solve (b, info, rcond, 0);
6971 solve_singularity_handler sing_handler)
const
6974 return solve (mattype, b, err, rcond, sing_handler);
6982 return solve (b, info, rcond, 0);
6990 return solve (b, info, rcond, 0);
6997 return solve (b, info, rcond, 0);
7003 solve_singularity_handler sing_handler)
const
7006 return solve (mattype, b, err, rcond, sing_handler);
7013 return solve (b, info, rcond);
7020 return solve (b, info, rcond);
7025 double& rcond)
const
7027 return solve (b, info, rcond, 0);
7033 solve_singularity_handler sing_handler)
const
7036 return solve (tmp, info, rcond,
7037 sing_handler).
column (static_cast<octave_idx_type> (0));
7045 return solve (b, info, rcond, 0);
7053 return solve (b, info, rcond, 0);
7058 double& rcond)
const
7060 return solve (b, info, rcond, 0);
7066 solve_singularity_handler sing_handler)
const
7069 return solve (tmp, info, rcond,
7070 sing_handler).
column (static_cast<octave_idx_type> (0));
7191 double r_val = val.real ();
7192 double i_val = val.imag ();
7194 if (r_val > max_val)
7197 if (i_val > max_val)
7200 if (r_val < min_val)
7203 if (i_val < min_val)
7249 if ((
rows () == 1 && dim == -1) || dim == 1)
7254 (
cidx (j+1) -
cidx (j) < nr ? 0.0 : 1.0), 1.0);
7268 Complex d = data (i); \
7269 tmp[ridx (i)] += d * conj (d)
7272 Complex d = data (i); \
7273 tmp[j] += d * conj (d)
7319 os << a.
ridx (
i) + 1 <<
" " << j + 1 <<
" ";
7333 return read_sparse_matrix<elt_type> (
is,
a, octave_read_value<Complex>);
7418 return do_mul_dm_sm<SparseComplexMatrix> (
d,
a);
7423 return do_mul_sm_dm<SparseComplexMatrix> (
a,
d);
7429 return do_mul_dm_sm<SparseComplexMatrix> (
d,
a);
7434 return do_mul_sm_dm<SparseComplexMatrix> (
a,
d);
7440 return do_mul_dm_sm<SparseComplexMatrix> (
d,
a);
7445 return do_mul_sm_dm<SparseComplexMatrix> (
a,
d);
7451 return do_add_dm_sm<SparseComplexMatrix> (
d,
a);
7456 return do_add_dm_sm<SparseComplexMatrix> (
d,
a);
7461 return do_add_dm_sm<SparseComplexMatrix> (
d,
a);
7466 return do_add_sm_dm<SparseComplexMatrix> (
a,
d);
7471 return do_add_sm_dm<SparseComplexMatrix> (
a,
d);
7476 return do_add_sm_dm<SparseComplexMatrix> (
a,
d);
7482 return do_sub_dm_sm<SparseComplexMatrix> (
d,
a);
7487 return do_sub_dm_sm<SparseComplexMatrix> (
d,
a);
7492 return do_sub_dm_sm<SparseComplexMatrix> (
d,
a);
7497 return do_sub_sm_dm<SparseComplexMatrix> (
a,
d);
7502 return do_sub_sm_dm<SparseComplexMatrix> (
a,
d);
7507 return do_sub_sm_dm<SparseComplexMatrix> (
a,
d);
7526 #define EMPTY_RETURN_CHECK(T) \
7527 if (nr == 0 || nc == 0) \
7570 if (a_nr == b_nr && a_nc == b_nc)
7580 bool ja_lt_max = ja < ja_max;
7584 bool jb_lt_max = jb < jb_max;
7586 while (ja_lt_max || jb_lt_max)
7589 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
7599 ja_lt_max= ja < ja_max;
7601 else if ((! ja_lt_max)
7602 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
7612 jb_lt_max= jb < jb_max;
7624 ja_lt_max= ja < ja_max;
7626 jb_lt_max= jb < jb_max;
7636 if (a_nr == 0 || a_nc == 0)
7638 else if (b_nr == 0 || b_nc == 0)
7687 if (a_nr == b_nr && a_nc == b_nc)
7697 bool ja_lt_max = ja < ja_max;
7701 bool jb_lt_max = jb < jb_max;
7703 while (ja_lt_max || jb_lt_max)
7706 if ((! jb_lt_max) || (ja_lt_max && (a.
ridx (ja) < b.
ridx (jb))))
7716 ja_lt_max= ja < ja_max;
7718 else if ((! ja_lt_max)
7719 || (jb_lt_max && (b.
ridx (jb) < a.
ridx (ja))))
7729 jb_lt_max= jb < jb_max;
7741 ja_lt_max= ja < ja_max;
7743 jb_lt_max= jb < jb_max;
7753 if (a_nr == 0 || a_nc == 0)
7755 else if (b_nr == 0 || b_nc == 0)
friend SparseComplexMatrix conj(const SparseComplexMatrix &a)
SparseComplexMatrix cumsum(int dim=-1) const
octave_idx_type * xridx(void)
ComplexMatrix bsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
octave_idx_type cols(void) const
bool operator==(const SparseComplexMatrix &a) const
MSparse< T > ipermute(const Array< octave_idx_type > &vec) const
std::istream & operator>>(std::istream &is, SparseComplexMatrix &a)
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
SparseComplexMatrix reshape(const dim_vector &new_dims) const
bool any_element_is_inf_or_nan(void) const
MSparse< T > & insert(const Sparse< T > &a, octave_idx_type r, octave_idx_type c)
SparseComplexMatrix cumprod(int dim=-1) const
octave_idx_type rows(void) const
#define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN)
SparseMatrix Q(void) const
ComplexMatrix solve(MatrixType &typ, const Matrix &b) const
ComplexMatrix matrix_value(void) const
const octave_base_value const Array< octave_idx_type > & ra_idx
SparseComplexMatrix max(const Complex &c, const SparseComplexMatrix &m)
Complex & elem(octave_idx_type n)
SparseComplexMatrix permute(const Array< octave_idx_type > &vec, bool inv=false) const
SparseMatrix transpose(void) const
template SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseComplexMatrix, SparseMatrix >(const SparseComplexMatrix &, const SparseMatrix &, octave_idx_type &)
dim_vector dims(void) const
#define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, INIT_VAL, MT_RESULT)
SparseMatrix Pc(void) const
std::ostream & operator<<(std::ostream &os, const SparseComplexMatrix &a)
ComplexMatrix dsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
identity matrix If supplied two scalar respectively For allows like xample val
SparseMatrix Pr(void) const
SparseComplexMatrix inverse(void) const
void SparseCholError(int status, char *file, int line, char *message)
MSparse< T > squeeze(void) const
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
#define F77_DBLE_CMPLX_ARG(x)
SM octinternal_do_mul_pm_sm(const PermMatrix &p, const SM &a)
SparseBoolMatrix operator!(void) const
ComplexColumnVector column(octave_idx_type i) const
#define SPARSE_SMSM_BOOL_OPS(M1, M2, ZERO)
void octave_write_complex(std::ostream &os, const Complex &c)
octave_idx_type * xcidx(void)
SparseComplexMatrix prod(int dim=-1) const
RowVector row(octave_idx_type i) const
SparseComplexMatrix concat(const SparseComplexMatrix &rb, const Array< octave_idx_type > &ra_idx)
SM octinternal_do_mul_sm_pm(const SM &a, const PermMatrix &p)
#define SPARSE_FULL_TRANS_MUL(RET_TYPE, EL_TYPE, ZERO, CONJ_OP)
SparseComplexMatrix operator*(const SparseComplexMatrix &m, const SparseMatrix &a)
octave_idx_type * cidx(void)
T & elem(octave_idx_type n)
bool is_hermitian(void) const
#define SPARSE_REDUCTION_OP(RET_TYPE, EL_TYPE, OP, INIT_VAL, MT_RESULT)
octave_idx_type columns(void) const
bool is_dense(void) const
void * factorize(octave_idx_type &err, double &rcond, Matrix &Control, Matrix &Info, solve_singularity_handler sing_handler, bool calc_cond) const
SparseBoolMatrix any(int dim=-1) const
#define FULL_SPARSE_MUL(RET_TYPE, EL_TYPE, ZERO)
#define UMFPACK_ZNAME(name)
#define F77_XFCN(f, F, args)
void err_nan_to_logical_conversion(void)
static double get_key(const std::string &key)
SparseComplexMatrix max(int dim=-1) const
SparseComplexMatrix min(int dim=-1) const
int first_non_singleton(int def=0) const
octave_idx_type rows(void) const
#define SPARSE_SMS_CMP_OPS(M, MZ, CM, S, SZ, CS)
#define SPARSE_FULL_MUL(RET_TYPE, EL_TYPE, ZERO)
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
ComplexMatrix mul_trans(const ComplexMatrix &m, const SparseComplexMatrix &a)
octave_idx_type * triangular_perm(void) const
octave_idx_type nnz(void) const
Actual number of nonzero terms.
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
OCTAVE_EXPORT octave_value_list search each directory of the loadpath for element of the cell array and return the first that matches If the second optional argument return a cell array containing the list of all files that have the same name in the path If no files are found
SparseComplexMatrix squeeze(void) const
SparseComplexMatrix transpose(void) const
#define SPARSE_ALL_OP(DIM)
SparseComplexMatrix(void)
bool mx_inline_all_real(size_t n, const std::complex< T > *x)
SparseComplexMatrix operator-(const ComplexDiagMatrix &d, const SparseMatrix &a)
#define SPARSE_ANY_OP(DIM)
static const Complex Complex_NaN_result(octave::numeric_limits< double >::NaN(), octave::numeric_limits< double >::NaN())
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
bool is_hermitian(void) const
template ComplexMatrix dmsolve< ComplexMatrix, SparseComplexMatrix, ComplexMatrix >(const SparseComplexMatrix &, const ComplexMatrix &, octave_idx_type &)
int type(bool quiet=true)
OCTAVE_EXPORT octave_value_list isdir nd deftypefn *std::string nm
ComplexRowVector row(octave_idx_type i) const
#define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN)
#define FULL_SPARSE_MUL_TRANS(RET_TYPE, EL_TYPE, ZERO, CONJ_OP)
ComplexMatrix herm_mul(const SparseComplexMatrix &m, const ComplexMatrix &a)
void change_capacity(octave_idx_type nz)
nd deftypefn *octave_map m
SparseComplexMatrix dinverse(MatrixType &mattyp, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
MatrixType transpose(void) const
#define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE)
Sparse< T > maybe_compress(bool remove_zeros=false)
void resize(octave_idx_type r, octave_idx_type c)
MSparse< T > permute(const Array< octave_idx_type > &vec, bool inv=false) const
SparseMatrix abs(void) const
ComplexMatrix trisolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
ComplexDET determinant(void) const
void resize(const dim_vector &dv, const T &rfv)
ComplexMatrix mul_herm(const ComplexMatrix &m, const SparseComplexMatrix &a)
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
MSparse< T > reshape(const dim_vector &new_dims) const
#define EMPTY_RETURN_CHECK(T)
bool too_large_for_float(void) const
void mark_as_unsymmetric(void)
SparseComplexMatrix sumsq(int dim=-1) const
SparseComplexMatrix min(const Complex &c, const SparseComplexMatrix &m)
SparseComplexMatrix hermitian(void) const
SparseComplexMatrix tinverse(MatrixType &mattyp, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
SparseComplexMatrix & insert(const SparseComplexMatrix &a, octave_idx_type r, octave_idx_type c)
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
void mark_as_rectangular(void)
octave_idx_type nnz(void) const
Count nonzero elements.
SparseBoolMatrix all(int dim=-1) const
With real return the complex result
#define SPARSE_SMSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2)
T & xelem(octave_idx_type n)
friend class ComplexMatrix
bool all_elements_are_real(void) const
octave_idx_type cols(void) const
SparseComplexMatrix operator+(const ComplexDiagMatrix &d, const SparseMatrix &a)
int SparseCholPrint(const char *fmt,...)
octave_idx_type * ridx(void)
MSparse< T > diag(octave_idx_type k=0) const
ComplexMatrix utsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
#define END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
bool any_element_is_nan(void) const
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
ComplexMatrix fsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
ComplexMatrix trans_mul(const SparseComplexMatrix &m, const ComplexMatrix &a)
#define SPARSE_SSM_CMP_OPS(S, SZ, SC, M, MZ, MC)
=val(i)}if ode{val(i)}occurs in table i
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
bool test_any(F fcn) const
SparseComplexMatrix sum(int dim=-1) const
octave_idx_type ndims(void) const
Number of dimensions.
SparseComplexMatrix conj(const SparseComplexMatrix &a)
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol NaN(Not a Number).NaN is the result of operations which do not produce a well defined 0 result.Common operations which produce a NaN are arithmetic with infinity ex($\infty-\infty $)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
SparseComplexMatrix diag(octave_idx_type k=0) const
bool xtoo_large_for_float(double x)
base_det< Complex > ComplexDET
OCTAVE_EXPORT octave_value_list error nd deftypefn *const octave_scalar_map err
#define SPARSE_SMS_BOOL_OPS(M, S, ZERO)
ComplexMatrix ltsolve(MatrixType &typ, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
#define BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
Matrix sum(int dim=-1) const
std::complex< double > Complex
const T * fortran_vec(void) const
#define SPARSE_SSM_BOOL_OPS(S, M, ZERO)
ColumnVector real(const ComplexColumnVector &a)
octave_idx_type cols(void) const
write the output to stdout if nargout is
SparseComplexMatrix ipermute(const Array< octave_idx_type > &vec) const
Vector representing the dimensions (size) of an Array.
bool operator!=(const SparseComplexMatrix &a) const
void warn_singular_matrix(double rcond)
bool all_integers(double &max_val, double &min_val) const
ComplexColumnVector column(octave_idx_type i) const
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)
octave_idx_type length(void) const
F77_RET_T const F77_INT F77_CMPLX * A
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
Array< T > array_value(void) const