25 #if defined (HAVE_CONFIG_H) 41 #include "mx-fcm-fs.h" 43 #include "mx-fs-fcm.h" 68 #if ! defined (USE_QRSOLVE) 120 if (nr != nr_a || nc != nc_a || nz != nz_a)
137 return !(*
this ==
a);
146 if (nr == nc && nr > 0)
188 return max (dummy_idx, dim);
212 if (nr == 0 || nc == 0 || dim >=
dv.
ndims ())
223 if (
ridx (
i) != idx_j)
276 if (nr == 0 || nc == 0 || dim >=
dv.
ndims ())
304 idx_arg.
elem (ir) = j;
310 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
320 if (idx_arg(j) == -1)
345 return min (dummy_idx, dim);
369 if (nr == 0 || nc == 0 || dim >=
dv.
ndims ())
380 if (
ridx (
i) != idx_j)
433 if (nr == 0 || nc == 0 || dim >=
dv.
ndims ())
461 idx_arg.
elem (ir) = j;
467 if (idx_arg.
elem (j) == -1 ||
elem (j, idx_arg.
elem (j)) != 0.)
477 if (idx_arg(j) == -1)
583 if (rb.
rows () > 0 && rb.
cols () > 0)
593 if (rb.
rows () > 0 && rb.
cols () > 0)
664 return inverse (mattype, info, rcond, 0, 0);
672 return inverse (mattype, info, rcond, 0, 0);
679 return inverse (mattype, info, rcond, 0, 0);
684 double& rcond,
const bool,
685 const bool calccond)
const 693 if (nr == 0 || nc == 0 || nr != nc)
694 (*current_liboctave_error_handler) (
"inverse requires square matrix");
697 int typ = mattype.
type ();
701 (*current_liboctave_error_handler) (
"incorrect matrix type");
734 double& rcond,
const bool,
735 const bool calccond)
const 743 if (nr == 0 || nc == 0 || nr != nc)
744 (*current_liboctave_error_handler) (
"inverse requires square matrix");
747 int typ = mattype.
type ();
752 (*current_liboctave_error_handler) (
"incorrect matrix type");
755 double ainvnorm = 0.;
786 retval.change_capacity (nz2);
804 (*current_liboctave_error_handler) (
"division by zero");
809 rpX =
retval.xridx (colXp);
823 while (rpX < j && rpU < j && colXp < cx && colUp < nz);
827 colUp =
cidx (j+1) - 1;
831 if (pivot == 0. ||
ridx (colUp) != j)
832 (*current_liboctave_error_handler) (
"division by zero");
839 retval.change_capacity (nz2);
843 retval.xdata (cx) = v / pivot;
851 colUp =
cidx (
i+1) - 1;
855 if (pivot == 0. ||
ridx (colUp) !=
i)
856 (*current_liboctave_error_handler) (
"division by zero");
860 retval.xdata (j) /= pivot;
920 (*current_liboctave_error_handler) (
"division by zero");
928 colUp =
cidx (perm[iidx]+1) - 1;
930 colUp =
cidx (perm[iidx]);
934 (*current_liboctave_error_handler) (
"division by zero");
947 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
948 retval.change_capacity (nz2);
956 retval.xdata (cx++) = work[j];
977 rcond = 1. / ainvnorm / anorm;
985 double& rcond,
bool,
bool calc_cond)
const 987 int typ = mattype.
type (
false);
991 typ = mattype.
type (*
this);
994 ret =
dinverse (mattype, info, rcond,
true, calc_cond);
1008 rcond = fact.
rcond ();
1016 ret =
Q * InvL.
hermitian () * InvL *
Q.transpose ();
1036 rcond = fact.
rcond ();
1072 #if defined (HAVE_UMFPACK) 1077 if (nr == 0 || nc == 0 || nr != nc)
1086 Matrix Control (UMFPACK_CONTROL, 1);
1092 Control (UMFPACK_PRL) =
tmp;
1097 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
1098 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
1104 Control (UMFPACK_FIXQ) =
tmp;
1107 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1118 reinterpret_cast<const double *
> (Ax),
1119 nullptr, 1, control);
1122 Matrix Info (1, UMFPACK_INFO);
1127 reinterpret_cast<const double *
> (Ax),
1128 nullptr,
nullptr, &Symbolic, control, info);
1137 (*current_liboctave_error_handler)
1138 (
"SparseComplexMatrix::determinant symbolic factorization failed");
1148 reinterpret_cast<const double *
> (Ax),
1149 nullptr, Symbolic, &Numeric, control, info);
1152 rcond = Info (UMFPACK_RCOND);
1161 (*current_liboctave_error_handler)
1162 (
"SparseComplexMatrix::determinant numeric factorization failed");
1170 status =
UMFPACK_ZNAME (get_determinant) (c10,
nullptr, &e10,
1178 (*current_liboctave_error_handler)
1179 (
"SparseComplexMatrix::determinant error calculating determinant");
1191 octave_unused_parameter (
err);
1192 octave_unused_parameter (rcond);
1194 (*current_liboctave_error_handler)
1195 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
1205 solve_singularity_handler,
bool calc_cond)
const 1214 if (nr !=
b.rows ())
1216 (
"matrix dimension mismatch solution of linear equations");
1218 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1223 int typ = mattype.
type ();
1227 (*current_liboctave_error_handler) (
"incorrect matrix type");
1252 rcond = dmin / dmax;
1264 solve_singularity_handler,
1265 bool calc_cond)
const 1274 if (nr !=
b.rows ())
1276 (
"matrix dimension mismatch solution of linear equations");
1278 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1283 int typ = mattype.
type ();
1287 (*current_liboctave_error_handler) (
"incorrect matrix type");
1300 if (
b.ridx (
i) >=
nm)
1315 for (
k =
b.cidx (j);
k <
b.cidx (j+1);
k++)
1342 rcond = dmin / dmax;
1354 solve_singularity_handler,
1355 bool calc_cond)
const 1364 if (nr !=
b.rows ())
1366 (
"matrix dimension mismatch solution of linear equations");
1368 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1373 int typ = mattype.
type ();
1377 (*current_liboctave_error_handler) (
"incorrect matrix type");
1402 rcond = dmin / dmax;
1414 solve_singularity_handler,
1415 bool calc_cond)
const 1424 if (nr !=
b.rows ())
1426 (
"matrix dimension mismatch solution of linear equations");
1428 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1433 int typ = mattype.
type ();
1437 (*current_liboctave_error_handler) (
"incorrect matrix type");
1450 if (
b.ridx (
i) >=
nm)
1465 for (
k =
b.cidx (j);
k <
b.cidx (j+1);
k++)
1492 rcond = dmin / dmax;
1504 solve_singularity_handler sing_handler,
1505 bool calc_cond)
const 1514 if (nr !=
b.rows ())
1516 (
"matrix dimension mismatch solution of linear equations");
1518 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1523 int typ = mattype.
type ();
1527 (*current_liboctave_error_handler) (
"incorrect matrix type");
1530 double ainvnorm = 0.;
1570 goto triangular_error;
1576 i <
cidx (kidx+1)-1;
i++)
1579 work[iidx] = work[iidx] -
tmp *
data (
i);
1607 i <
cidx (iidx+1)-1;
i++)
1610 work[idx2] = work[idx2] -
tmp *
data (
i);
1620 if (atmp > ainvnorm)
1623 rcond = 1. / ainvnorm / anorm;
1646 goto triangular_error;
1654 work[iidx] = work[iidx] -
tmp *
data (
i);
1683 work[iidx] = work[iidx] -
tmp *
data (
i);
1693 if (atmp > ainvnorm)
1696 rcond = 1. / ainvnorm / anorm;
1705 sing_handler (rcond);
1712 volatile double rcond_plus_one = rcond + 1.0;
1720 sing_handler (rcond);
1734 solve_singularity_handler sing_handler,
1735 bool calc_cond)
const 1744 if (nr !=
b.rows ())
1746 (
"matrix dimension mismatch solution of linear equations");
1748 if (nr == 0 || nc == 0 ||
b.cols () == 0)
1753 int typ = mattype.
type ();
1757 (*current_liboctave_error_handler) (
"incorrect matrix type");
1760 double ainvnorm = 0.;
1797 work[
b.ridx (
i)] =
b.data (
i);
1809 goto triangular_error;
1815 i <
cidx (kidx+1)-1;
i++)
1818 work[iidx] = work[iidx] -
tmp *
data (
i);
1830 if (ii + new_nnz > x_nz)
1839 if (work[rperm[
i]] != 0.)
1842 retval.xdata (ii++) = work[rperm[
i]];
1847 retval.maybe_compress ();
1868 i <
cidx (iidx+1)-1;
i++)
1871 work[idx2] = work[idx2] -
tmp *
data (
i);
1881 if (atmp > ainvnorm)
1884 rcond = 1. / ainvnorm / anorm;
1896 work[
b.ridx (
i)] =
b.data (
i);
1906 goto triangular_error;
1914 work[iidx] = work[iidx] -
tmp *
data (
i);
1926 if (ii + new_nnz > x_nz)
1938 retval.xdata (ii++) = work[
i];
1943 retval.maybe_compress ();
1965 work[iidx] = work[iidx] -
tmp *
data (
i);
1975 if (atmp > ainvnorm)
1978 rcond = 1. / ainvnorm / anorm;
1987 sing_handler (rcond);
1994 volatile double rcond_plus_one = rcond + 1.0;
2002 sing_handler (rcond);
2015 solve_singularity_handler sing_handler,
2016 bool calc_cond)
const 2025 if (nr !=
b.rows ())
2027 (
"matrix dimension mismatch solution of linear equations");
2029 if (nr == 0 || nc == 0 ||
b.cols () == 0)
2034 int typ = mattype.
type ();
2038 (*current_liboctave_error_handler) (
"incorrect matrix type");
2041 double ainvnorm = 0.;
2081 goto triangular_error;
2087 i <
cidx (kidx+1)-1;
i++)
2090 work[iidx] = work[iidx] -
tmp *
data (
i);
2118 i <
cidx (iidx+1)-1;
i++)
2121 work[idx2] = work[idx2] -
tmp *
data (
i);
2131 if (atmp > ainvnorm)
2134 rcond = 1. / ainvnorm / anorm;
2157 goto triangular_error;
2165 work[iidx] = work[iidx] -
tmp *
data (
i);
2194 work[iidx] = work[iidx] -
tmp *
data (
i);
2204 if (atmp > ainvnorm)
2207 rcond = 1. / ainvnorm / anorm;
2216 sing_handler (rcond);
2223 volatile double rcond_plus_one = rcond + 1.0;
2231 sing_handler (rcond);
2245 solve_singularity_handler sing_handler,
2246 bool calc_cond)
const 2255 if (nr !=
b.rows ())
2257 (
"matrix dimension mismatch solution of linear equations");
2259 if (nr == 0 || nc == 0 ||
b.cols () == 0)
2264 int typ = mattype.
type ();
2268 (*current_liboctave_error_handler) (
"incorrect matrix type");
2271 double ainvnorm = 0.;
2308 work[
b.ridx (
i)] =
b.data (
i);
2320 goto triangular_error;
2326 i <
cidx (kidx+1)-1;
i++)
2329 work[iidx] = work[iidx] -
tmp *
data (
i);
2341 if (ii + new_nnz > x_nz)
2350 if (work[rperm[
i]] != 0.)
2353 retval.xdata (ii++) = work[rperm[
i]];
2358 retval.maybe_compress ();
2379 i <
cidx (iidx+1)-1;
i++)
2382 work[idx2] = work[idx2] -
tmp *
data (
i);
2392 if (atmp > ainvnorm)
2395 rcond = 1. / ainvnorm / anorm;
2407 work[
b.ridx (
i)] =
b.data (
i);
2417 goto triangular_error;
2425 work[iidx] = work[iidx] -
tmp *
data (
i);
2437 if (ii + new_nnz > x_nz)
2449 retval.xdata (ii++) = work[
i];
2454 retval.maybe_compress ();
2476 work[iidx] = work[iidx] -
tmp *
data (
i);
2486 if (atmp > ainvnorm)
2489 rcond = 1. / ainvnorm / anorm;
2498 sing_handler (rcond);
2505 volatile double rcond_plus_one = rcond + 1.0;
2513 sing_handler (rcond);
2527 solve_singularity_handler sing_handler,
2528 bool calc_cond)
const 2537 if (nr !=
b.rows ())
2539 (
"matrix dimension mismatch solution of linear equations");
2541 if (nr == 0 || nc == 0 ||
b.cols () == 0)
2546 int typ = mattype.
type ();
2550 (*current_liboctave_error_handler) (
"incorrect matrix type");
2553 double ainvnorm = 0.;
2581 work[perm[
i]] =
b(
i,j);
2591 if (perm[
ridx (
i)] < minr)
2593 minr = perm[
ridx (
i)];
2597 if (minr !=
k ||
data (mini) == 0.)
2600 goto triangular_error;
2611 work[iidx] = work[iidx] -
tmp *
data (
i);
2639 if (perm[
ridx (
i)] < minr)
2641 minr = perm[
ridx (
i)];
2654 work[iidx] = work[iidx] -
tmp *
data (
i);
2665 if (atmp > ainvnorm)
2668 rcond = 1. / ainvnorm / anorm;
2689 goto triangular_error;
2697 work[iidx] = work[iidx] -
tmp *
data (
i);
2726 work[iidx] = work[iidx] -
tmp *
data (
i);
2736 if (atmp > ainvnorm)
2739 rcond = 1. / ainvnorm / anorm;
2747 sing_handler (rcond);
2754 volatile double rcond_plus_one = rcond + 1.0;
2762 sing_handler (rcond);
2776 solve_singularity_handler sing_handler,
2777 bool calc_cond)
const 2787 if (nr !=
b.rows ())
2789 (
"matrix dimension mismatch solution of linear equations");
2791 if (nr == 0 || nc == 0 ||
b.cols () == 0)
2796 int typ = mattype.
type ();
2800 (*current_liboctave_error_handler) (
"incorrect matrix type");
2803 double ainvnorm = 0.;
2836 work[perm[
b.ridx (
i)]] =
b.data (
i);
2846 if (perm[
ridx (
i)] < minr)
2848 minr = perm[
ridx (
i)];
2852 if (minr !=
k ||
data (mini) == 0.)
2855 goto triangular_error;
2866 work[iidx] = work[iidx] -
tmp *
data (
i);
2878 if (ii + new_nnz > x_nz)
2890 retval.xdata (ii++) = work[
i];
2895 retval.maybe_compress ();
2916 if (perm[
ridx (
i)] < minr)
2918 minr = perm[
ridx (
i)];
2931 work[iidx] = work[iidx] -
tmp *
data (
i);
2942 if (atmp > ainvnorm)
2945 rcond = 1. / ainvnorm / anorm;
2957 work[
b.ridx (
i)] =
b.data (
i);
2966 goto triangular_error;
2974 work[iidx] = work[iidx] -
tmp *
data (
i);
2986 if (ii + new_nnz > x_nz)
2998 retval.xdata (ii++) = work[
i];
3003 retval.maybe_compress ();
3026 work[iidx] = work[iidx] -
tmp *
data (
i);
3036 if (atmp > ainvnorm)
3039 rcond = 1. / ainvnorm / anorm;
3048 sing_handler (rcond);
3055 volatile double rcond_plus_one = rcond + 1.0;
3063 sing_handler (rcond);
3077 solve_singularity_handler sing_handler,
3078 bool calc_cond)
const 3087 if (nr !=
b.rows ())
3089 (
"matrix dimension mismatch solution of linear equations");
3091 if (nr == 0 || nc == 0 ||
b.cols () == 0)
3096 int typ = mattype.
type ();
3100 (*current_liboctave_error_handler) (
"incorrect matrix type");
3103 double ainvnorm = 0.;
3131 work[perm[
i]] =
b(
i,j);
3141 if (perm[
ridx (
i)] < minr)
3143 minr = perm[
ridx (
i)];
3147 if (minr !=
k ||
data (mini) == 0.)
3150 goto triangular_error;
3161 work[iidx] = work[iidx] -
tmp *
data (
i);
3189 if (perm[
ridx (
i)] < minr)
3191 minr = perm[
ridx (
i)];
3204 work[iidx] = work[iidx] -
tmp *
data (
i);
3215 if (atmp > ainvnorm)
3218 rcond = 1. / ainvnorm / anorm;
3240 goto triangular_error;
3248 work[iidx] = work[iidx] -
tmp *
data (
i);
3278 work[iidx] = work[iidx] -
tmp *
data (
i);
3288 if (atmp > ainvnorm)
3291 rcond = 1. / ainvnorm / anorm;
3300 sing_handler (rcond);
3307 volatile double rcond_plus_one = rcond + 1.0;
3315 sing_handler (rcond);
3329 solve_singularity_handler sing_handler,
3330 bool calc_cond)
const 3339 if (nr !=
b.rows ())
3341 (
"matrix dimension mismatch solution of linear equations");
3343 if (nr == 0 || nc == 0 ||
b.cols () == 0)
3348 int typ = mattype.
type ();
3352 (*current_liboctave_error_handler) (
"incorrect matrix type");
3355 double ainvnorm = 0.;
3388 work[perm[
b.ridx (
i)]] =
b.data (
i);
3398 if (perm[
ridx (
i)] < minr)
3400 minr = perm[
ridx (
i)];
3404 if (minr !=
k ||
data (mini) == 0.)
3407 goto triangular_error;
3418 work[iidx] = work[iidx] -
tmp *
data (
i);
3430 if (ii + new_nnz > x_nz)
3442 retval.xdata (ii++) = work[
i];
3447 retval.maybe_compress ();
3468 if (perm[
ridx (
i)] < minr)
3470 minr = perm[
ridx (
i)];
3483 work[iidx] = work[iidx] -
tmp *
data (
i);
3494 if (atmp > ainvnorm)
3497 rcond = 1. / ainvnorm / anorm;
3509 work[
b.ridx (
i)] =
b.data (
i);
3518 goto triangular_error;
3526 work[iidx] = work[iidx] -
tmp *
data (
i);
3538 if (ii + new_nnz > x_nz)
3550 retval.xdata (ii++) = work[
i];
3555 retval.maybe_compress ();
3578 work[iidx] = work[iidx] -
tmp *
data (
i);
3588 if (atmp > ainvnorm)
3591 rcond = 1. / ainvnorm / anorm;
3600 sing_handler (rcond);
3607 volatile double rcond_plus_one = rcond + 1.0;
3615 sing_handler (rcond);
3629 solve_singularity_handler sing_handler,
3630 bool calc_cond)
const 3638 if (nr != nc || nr !=
b.rows ())
3640 (
"matrix dimension mismatch solution of linear equations");
3642 if (nr == 0 ||
b.cols () == 0)
3645 (*current_liboctave_error_handler)
3646 (
"calculation of condition number not implemented");
3650 volatile int typ = mattype.
type ();
3684 else if (
ridx (
i) == j + 1)
3695 F77_INT tmp_nr = octave::to_f77_int (nr);
3728 DL[j] =
data (ii++);
3729 DU[j] =
data (ii++);
3731 D[nc-1] =
data (ii);
3748 else if (
ridx (
i) == j + 1)
3750 else if (
ridx (
i) == j - 1)
3761 F77_INT tmp_nr = octave::to_f77_int (nr);
3778 sing_handler (rcond);
3789 (*current_liboctave_error_handler) (
"incorrect matrix type");
3798 solve_singularity_handler sing_handler,
3799 bool calc_cond)
const 3807 if (nr != nc || nr !=
b.rows ())
3809 (
"matrix dimension mismatch solution of linear equations");
3811 if (nr == 0 ||
b.cols () == 0)
3814 (*current_liboctave_error_handler)
3815 (
"calculation of condition number not implemented");
3819 int typ = mattype.
type ();
3840 DL[j] =
data (ii++);
3841 DU[j] =
data (ii++);
3843 D[nc-1] =
data (ii);
3860 else if (
ridx (
i) == j + 1)
3862 else if (
ridx (
i) == j - 1)
3867 F77_INT tmp_nr = octave::to_f77_int (nr);
3886 sing_handler (rcond);
3910 work[
b.ridx (
i)] =
b.data (
i);
3913 (F77_CONST_CHAR_ARG2 (&job, 1),
3919 F77_CHAR_ARG_LEN (1)));
3930 if (ii + new_nnz > x_nz)
3942 retval.xdata (ii++) = work[
i];
3947 retval.maybe_compress ();
3951 (*current_liboctave_error_handler) (
"incorrect matrix type");
3960 solve_singularity_handler sing_handler,
3961 bool calc_cond)
const 3969 if (nr != nc || nr !=
b.rows ())
3971 (
"matrix dimension mismatch solution of linear equations");
3973 if (nr == 0 ||
b.cols () == 0)
3976 (*current_liboctave_error_handler)
3977 (
"calculation of condition number not implemented");
3981 volatile int typ = mattype.
type ();
4015 else if (
ridx (
i) == j + 1)
4028 F77_INT tmp_nr = octave::to_f77_int (nr);
4059 DL[j] =
data (ii++);
4060 DU[j] =
data (ii++);
4062 D[nc-1] =
data (ii);
4079 else if (
ridx (
i) == j + 1)
4081 else if (
ridx (
i) == j - 1)
4094 F77_INT tmp_nr = octave::to_f77_int (nr);
4111 sing_handler (rcond);
4119 (*current_liboctave_error_handler) (
"incorrect matrix type");
4129 solve_singularity_handler sing_handler,
4130 bool calc_cond)
const 4138 if (nr != nc || nr !=
b.rows ())
4140 (
"matrix dimension mismatch solution of linear equations");
4142 if (nr == 0 ||
b.cols () == 0)
4145 (*current_liboctave_error_handler)
4146 (
"calculation of condition number not implemented");
4150 int typ = mattype.
type ();
4171 DL[j] =
data (ii++);
4172 DU[j] =
data (ii++);
4174 D[nc-1] =
data (ii);
4191 else if (
ridx (
i) == j + 1)
4193 else if (
ridx (
i) == j - 1)
4198 F77_INT tmp_nr = octave::to_f77_int (nr);
4217 sing_handler (rcond);
4245 (F77_CONST_CHAR_ARG2 (&job, 1),
4251 F77_CHAR_ARG_LEN (1)));
4259 (*current_liboctave_error_handler)
4260 (
"SparseComplexMatrix::solve solve failed");
4273 if (ii + new_nnz > x_nz)
4291 retval.maybe_compress ();
4295 (*current_liboctave_error_handler) (
"incorrect matrix type");
4304 solve_singularity_handler sing_handler,
4305 bool calc_cond)
const 4313 if (nr != nc || nr !=
b.rows ())
4315 (
"matrix dimension mismatch solution of linear equations");
4317 if (nr == 0 ||
b.cols () == 0)
4322 volatile int typ = mattype.
type ();
4336 for (
F77_INT j = 0; j < ldm; j++)
4338 tmp_data[ii++] = 0.;
4346 m_band(ri - j, j) =
data (
i);
4354 F77_INT tmp_nr = octave::to_f77_int (nr);
4359 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4361 F77_CHAR_ARG_LEN (1)));
4384 (F77_CONST_CHAR_ARG2 (&job, 1),
4387 F77_CHAR_ARG_LEN (1)));
4394 volatile double rcond_plus_one = rcond + 1.0;
4402 sing_handler (rcond);
4421 (F77_CONST_CHAR_ARG2 (&job, 1),
4424 F77_CHAR_ARG_LEN (1)));
4431 (*current_liboctave_error_handler)
4432 (
"SparseMatrix::solve solve failed");
4444 F77_INT ldm = n_upper + 2 * n_lower + 1;
4453 for (
F77_INT j = 0; j < ldm; j++)
4455 tmp_data[ii++] = 0.;
4460 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
4479 F77_INT tmp_nr = octave::to_f77_int (nr);
4483 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4485 ldm, pipvt, tmp_err));
4498 sing_handler (rcond);
4514 F77_INT tmp_nc = octave::to_f77_int (nc);
4517 (F77_CONST_CHAR_ARG2 (&job, 1),
4520 F77_CHAR_ARG_LEN (1)));
4527 volatile double rcond_plus_one = rcond + 1.0;
4535 sing_handler (rcond);
4555 (F77_CONST_CHAR_ARG2 (&job, 1),
4558 F77_CHAR_ARG_LEN (1)));
4565 (*current_liboctave_error_handler) (
"incorrect matrix type");
4574 solve_singularity_handler sing_handler,
4575 bool calc_cond)
const 4583 if (nr != nc || nr !=
b.rows ())
4585 (
"matrix dimension mismatch solution of linear equations");
4587 if (nr == 0 ||
b.cols () == 0)
4592 volatile int typ = mattype.
type ();
4607 for (
F77_INT j = 0; j < ldm; j++)
4609 tmp_data[ii++] = 0.;
4617 m_band(ri - j, j) =
data (
i);
4625 F77_INT tmp_nr = octave::to_f77_int (nr);
4630 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4632 F77_CHAR_ARG_LEN (1)));
4653 (F77_CONST_CHAR_ARG2 (&job, 1),
4656 F77_CHAR_ARG_LEN (1)));
4663 volatile double rcond_plus_one = rcond + 1.0;
4671 sing_handler (rcond);
4697 Bx[
i] =
b.elem (
i, j);
4700 (F77_CONST_CHAR_ARG2 (&job, 1),
4703 F77_CHAR_ARG_LEN (1)));
4710 (*current_liboctave_error_handler)
4711 (
"SparseComplexMatrix::solve solve failed");
4727 sz = x_nz + (
sz > 100 ?
sz : 100);
4738 retval.maybe_compress ();
4748 F77_INT ldm = n_upper + 2 * n_lower + 1;
4757 for (
F77_INT j = 0; j < ldm; j++)
4759 tmp_data[ii++] = 0.;
4764 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
4783 F77_INT tmp_nr = octave::to_f77_int (nr);
4787 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4789 ldm, pipvt, tmp_err));
4800 sing_handler (rcond);
4816 F77_INT tmp_nc = octave::to_f77_int (nc);
4819 (F77_CONST_CHAR_ARG2 (&job, 1),
4822 F77_CHAR_ARG_LEN (1)));
4829 volatile double rcond_plus_one = rcond + 1.0;
4837 sing_handler (rcond);
4864 i <
b.cidx (j+1);
i++)
4865 work[
b.ridx (
i)] =
b.data (
i);
4868 (F77_CONST_CHAR_ARG2 (&job, 1),
4871 F77_CHAR_ARG_LEN (1)));
4882 if (ii + new_nnz > x_nz)
4894 retval.xdata (ii++) = work[
i];
4899 retval.maybe_compress ();
4904 (*current_liboctave_error_handler) (
"incorrect matrix type");
4913 solve_singularity_handler sing_handler,
4914 bool calc_cond)
const 4922 if (nr != nc || nr !=
b.rows ())
4924 (
"matrix dimension mismatch solution of linear equations");
4926 if (nr == 0 ||
b.cols () == 0)
4931 volatile int typ = mattype.
type ();
4946 for (
F77_INT j = 0; j < ldm; j++)
4948 tmp_data[ii++] = 0.;
4956 m_band(ri - j, j) =
data (
i);
4964 F77_INT tmp_nr = octave::to_f77_int (nr);
4969 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4971 F77_CHAR_ARG_LEN (1)));
4994 (F77_CONST_CHAR_ARG2 (&job, 1),
4997 F77_CHAR_ARG_LEN (1)));
5004 volatile double rcond_plus_one = rcond + 1.0;
5012 sing_handler (rcond);
5030 (F77_CONST_CHAR_ARG2 (&job, 1),
5033 F77_CHAR_ARG_LEN (1)));
5040 (*current_liboctave_error_handler)
5041 (
"SparseComplexMatrix::solve solve failed");
5053 F77_INT ldm = n_upper + 2 * n_lower + 1;
5062 for (
F77_INT j = 0; j < ldm; j++)
5064 tmp_data[ii++] = 0.;
5069 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
5088 F77_INT tmp_nr = octave::to_f77_int (nr);
5092 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5094 ldm, pipvt, tmp_err));
5105 sing_handler (rcond);
5121 F77_INT tmp_nc = octave::to_f77_int (nc);
5124 (F77_CONST_CHAR_ARG2 (&job, 1),
5127 F77_CHAR_ARG_LEN (1)));
5134 volatile double rcond_plus_one = rcond + 1.0;
5142 sing_handler (rcond);
5161 (F77_CONST_CHAR_ARG2 (&job, 1),
5164 F77_CHAR_ARG_LEN (1)));
5171 (*current_liboctave_error_handler) (
"incorrect matrix type");
5180 solve_singularity_handler sing_handler,
5181 bool calc_cond)
const 5189 if (nr != nc || nr !=
b.rows ())
5191 (
"matrix dimension mismatch solution of linear equations");
5193 if (nr == 0 ||
b.cols () == 0)
5198 volatile int typ = mattype.
type ();
5213 for (
F77_INT j = 0; j < ldm; j++)
5215 tmp_data[ii++] = 0.;
5223 m_band(ri - j, j) =
data (
i);
5231 F77_INT tmp_nr = octave::to_f77_int (nr);
5236 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5238 F77_CHAR_ARG_LEN (1)));
5262 (F77_CONST_CHAR_ARG2 (&job, 1),
5265 F77_CHAR_ARG_LEN (1)));
5272 volatile double rcond_plus_one = rcond + 1.0;
5280 sing_handler (rcond);
5310 (F77_CONST_CHAR_ARG2 (&job, 1),
5313 F77_CHAR_ARG_LEN (1)));
5320 (*current_liboctave_error_handler)
5321 (
"SparseMatrix::solve solve failed");
5333 if (ii + new_nnz > x_nz)
5351 retval.maybe_compress ();
5361 F77_INT ldm = n_upper + 2 * n_lower + 1;
5370 for (
F77_INT j = 0; j < ldm; j++)
5372 tmp_data[ii++] = 0.;
5377 m_band(
ridx (
i) - j + n_lower + n_upper, j) =
data (
i);
5396 F77_INT tmp_nr = octave::to_f77_int (nr);
5400 F77_XFCN (zgbtrf, ZGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5402 ldm, pipvt, tmp_err));
5413 sing_handler (rcond);
5429 F77_INT tmp_nc = octave::to_f77_int (nc);
5432 (F77_CONST_CHAR_ARG2 (&job, 1),
5435 F77_CHAR_ARG_LEN (1)));
5442 volatile double rcond_plus_one = rcond + 1.0;
5450 sing_handler (rcond);
5478 i <
b.cidx (j+1);
i++)
5479 Bx[
b.ridx (
i)] =
b.data (
i);
5482 (F77_CONST_CHAR_ARG2 (&job, 1),
5485 F77_CHAR_ARG_LEN (1)));
5496 if (ii + new_nnz > x_nz)
5513 retval.maybe_compress ();
5518 (*current_liboctave_error_handler) (
"incorrect matrix type");
5527 solve_singularity_handler sing_handler,
5528 bool calc_cond)
const 5531 void *Numeric =
nullptr;
5534 #if defined (HAVE_UMFPACK) 5537 Control =
Matrix (UMFPACK_CONTROL, 1);
5543 Control (UMFPACK_PRL) =
tmp;
5547 Control (UMFPACK_SYM_PIVOT_TOLERANCE) =
tmp;
5548 Control (UMFPACK_PIVOT_TOLERANCE) =
tmp;
5554 Control (UMFPACK_FIXQ) =
tmp;
5567 reinterpret_cast<const double *
> (Ax),
5568 nullptr, 1, control);
5571 Info =
Matrix (1, UMFPACK_INFO);
5576 reinterpret_cast<const double *
> (Ax),
5577 nullptr,
nullptr, &Symbolic, control, info);
5587 (*current_liboctave_error_handler)
5588 (
"SparseComplexMatrix::solve symbolic factorization failed");
5597 reinterpret_cast<const double *
> (Ax),
5598 nullptr, Symbolic, &Numeric, control, info);
5602 rcond = Info (UMFPACK_RCOND);
5605 volatile double rcond_plus_one = rcond + 1.0;
5607 if (status == UMFPACK_WARNING_singular_matrix
5615 sing_handler (rcond);
5619 else if (status < 0)
5625 (*current_liboctave_error_handler)
5626 (
"SparseComplexMatrix::solve numeric factorization failed");
5641 octave_unused_parameter (rcond);
5642 octave_unused_parameter (Control);
5643 octave_unused_parameter (Info);
5644 octave_unused_parameter (sing_handler);
5645 octave_unused_parameter (calc_cond);
5647 (*current_liboctave_error_handler)
5648 (
"support for UMFPACK was unavailable or disabled when liboctave was built");
5658 solve_singularity_handler sing_handler,
5659 bool calc_cond)
const 5667 if (nr != nc || nr !=
b.rows ())
5669 (
"matrix dimension mismatch solution of linear equations");
5671 if (nr == 0 ||
b.cols () == 0)
5676 volatile int typ = mattype.
type ();
5681 #if defined (HAVE_CHOLMOD) 5682 cholmod_common Common;
5683 cholmod_common *cm = &Common;
5687 cm->prefer_zomplex =
false;
5693 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
5697 cm->print =
static_cast<int> (spu) + 2;
5698 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5702 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5703 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5705 cm->final_ll =
true;
5707 cholmod_sparse Astore;
5708 cholmod_sparse *
A = &Astore;
5719 #if defined (OCTAVE_ENABLE_64) 5720 A->itype = CHOLMOD_LONG;
5722 A->itype = CHOLMOD_INT;
5724 A->dtype = CHOLMOD_DOUBLE;
5726 A->xtype = CHOLMOD_COMPLEX;
5729 if (
A->x ==
nullptr)
5732 cholmod_dense Bstore;
5733 cholmod_dense *
B = &Bstore;
5734 B->nrow =
b.rows ();
5735 B->ncol =
b.cols ();
5737 B->nzmax =
B->nrow *
B->ncol;
5738 B->dtype = CHOLMOD_DOUBLE;
5739 B->xtype = CHOLMOD_REAL;
5741 B->x =
const_cast<double *
>(
b.fortran_vec ());
5742 if (
B->x ==
nullptr)
5763 volatile double rcond_plus_one = rcond + 1.0;
5771 sing_handler (rcond);
5797 static char blank_name[] =
" ";
5802 (*current_liboctave_warning_with_id_handler)
5803 (
"Octave:missing-dependency",
5804 "support for CHOLMOD was unavailable or disabled " 5805 "when liboctave was built");
5814 #if defined (HAVE_UMFPACK) 5817 sing_handler, calc_cond);
5822 Control (UMFPACK_IRSTEP) = 1;
5831 #if defined (UMFPACK_SEPARATE_SPLIT) 5832 const double *Bx =
b.fortran_vec ();
5844 #if defined (UMFPACK_SEPARATE_SPLIT) 5848 reinterpret_cast<const double *
> (Ax),
5850 reinterpret_cast<double *> (&Xx[iidx]),
5852 &Bx[iidx], Bz, Numeric,
5856 Bz[
i] =
b.elem (
i, j);
5861 reinterpret_cast<const double *
> (Ax),
5863 reinterpret_cast<double *> (&Xx[iidx]),
5865 reinterpret_cast<const double *
> (Bz),
5875 (*current_liboctave_error_handler)
5876 (
"SparseComplexMatrix::solve solve failed");
5891 octave_unused_parameter (rcond);
5892 octave_unused_parameter (sing_handler);
5893 octave_unused_parameter (calc_cond);
5895 (*current_liboctave_error_handler)
5896 (
"support for UMFPACK was unavailable or disabled " 5897 "when liboctave was built");
5901 (*current_liboctave_error_handler) (
"incorrect matrix type");
5910 solve_singularity_handler sing_handler,
5911 bool calc_cond)
const 5919 if (nr != nc || nr !=
b.rows ())
5921 (
"matrix dimension mismatch solution of linear equations");
5923 if (nr == 0 ||
b.cols () == 0)
5928 volatile int typ = mattype.
type ();
5933 #if defined (HAVE_CHOLMOD) 5934 cholmod_common Common;
5935 cholmod_common *cm = &Common;
5939 cm->prefer_zomplex =
false;
5945 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
5949 cm->print =
static_cast<int> (spu) + 2;
5950 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
5954 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5955 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5957 cm->final_ll =
true;
5959 cholmod_sparse Astore;
5960 cholmod_sparse *
A = &Astore;
5971 #if defined (OCTAVE_ENABLE_64) 5972 A->itype = CHOLMOD_LONG;
5974 A->itype = CHOLMOD_INT;
5976 A->dtype = CHOLMOD_DOUBLE;
5978 A->xtype = CHOLMOD_COMPLEX;
5981 if (
A->x ==
nullptr)
5984 cholmod_sparse Bstore;
5985 cholmod_sparse *
B = &Bstore;
5986 B->nrow =
b.rows ();
5987 B->ncol =
b.cols ();
5990 B->nzmax =
b.nnz ();
5994 #if defined (OCTAVE_ENABLE_64) 5995 B->itype = CHOLMOD_LONG;
5997 B->itype = CHOLMOD_INT;
5999 B->dtype = CHOLMOD_DOUBLE;
6001 B->xtype = CHOLMOD_REAL;
6004 if (
B->x ==
nullptr)
6025 volatile double rcond_plus_one = rcond + 1.0;
6033 sing_handler (rcond);
6048 (static_cast<octave_idx_type>(X->nrow),
6049 static_cast<octave_idx_type>(X->ncol),
6050 static_cast<octave_idx_type>(X->nzmax));
6052 j <= static_cast<octave_idx_type>(X->ncol); j++)
6055 j < static_cast<octave_idx_type>(X->nzmax); j++)
6065 static char blank_name[] =
" ";
6070 (*current_liboctave_warning_with_id_handler)
6071 (
"Octave:missing-dependency",
6072 "support for CHOLMOD was unavailable or disabled " 6073 "when liboctave was built");
6082 #if defined (HAVE_UMFPACK) 6085 sing_handler, calc_cond);
6090 Control (UMFPACK_IRSTEP) = 1;
6100 #if defined (UMFPACK_SEPARATE_SPLIT) 6121 #if defined (UMFPACK_SEPARATE_SPLIT) 6123 Bx[
i] =
b.elem (
i, j);
6128 reinterpret_cast<const double *
> (Ax),
6130 reinterpret_cast<double *> (Xx),
6132 Bx, Bz, Numeric, control,
6136 Bz[
i] =
b.elem (
i, j);
6141 reinterpret_cast<const double *
> (Ax),
6143 reinterpret_cast<double *> (Xx),
6145 reinterpret_cast<double *
> (Bz),
6155 (*current_liboctave_error_handler)
6156 (
"SparseComplexMatrix::solve solve failed");
6173 sz = x_nz + (
sz > 100 ?
sz : 100);
6184 retval.maybe_compress ();
6194 octave_unused_parameter (rcond);
6195 octave_unused_parameter (sing_handler);
6196 octave_unused_parameter (calc_cond);
6198 (*current_liboctave_error_handler)
6199 (
"support for UMFPACK was unavailable or disabled " 6200 "when liboctave was built");
6204 (*current_liboctave_error_handler) (
"incorrect matrix type");
6213 solve_singularity_handler sing_handler,
6214 bool calc_cond)
const 6222 if (nr != nc || nr !=
b.rows ())
6224 (
"matrix dimension mismatch solution of linear equations");
6226 if (nr == 0 ||
b.cols () == 0)
6231 volatile int typ = mattype.
type ();
6236 #if defined (HAVE_CHOLMOD) 6237 cholmod_common Common;
6238 cholmod_common *cm = &Common;
6242 cm->prefer_zomplex =
false;
6248 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6252 cm->print =
static_cast<int> (spu) + 2;
6253 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6257 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6258 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6260 cm->final_ll =
true;
6262 cholmod_sparse Astore;
6263 cholmod_sparse *
A = &Astore;
6274 #if defined (OCTAVE_ENABLE_64) 6275 A->itype = CHOLMOD_LONG;
6277 A->itype = CHOLMOD_INT;
6279 A->dtype = CHOLMOD_DOUBLE;
6281 A->xtype = CHOLMOD_COMPLEX;
6284 if (
A->x ==
nullptr)
6287 cholmod_dense Bstore;
6288 cholmod_dense *
B = &Bstore;
6289 B->nrow =
b.rows ();
6290 B->ncol =
b.cols ();
6292 B->nzmax =
B->nrow *
B->ncol;
6293 B->dtype = CHOLMOD_DOUBLE;
6294 B->xtype = CHOLMOD_COMPLEX;
6296 B->x =
const_cast<Complex *
>(
b.fortran_vec ());
6297 if (
B->x ==
nullptr)
6318 volatile double rcond_plus_one = rcond + 1.0;
6326 sing_handler (rcond);
6352 static char blank_name[] =
" ";
6357 (*current_liboctave_warning_with_id_handler)
6358 (
"Octave:missing-dependency",
6359 "support for CHOLMOD was unavailable or disabled " 6360 "when liboctave was built");
6369 #if defined (HAVE_UMFPACK) 6372 sing_handler, calc_cond);
6377 Control (UMFPACK_IRSTEP) = 1;
6386 const Complex *Bx =
b.fortran_vec ();
6397 reinterpret_cast<const double *
> (Ax),
6399 reinterpret_cast<double *> (&Xx[iidx]),
6401 reinterpret_cast<const double *
> (&Bx[iidx]),
6402 nullptr, Numeric, control, info);
6409 (*current_liboctave_error_handler)
6410 (
"SparseComplexMatrix::solve solve failed");
6425 octave_unused_parameter (rcond);
6426 octave_unused_parameter (sing_handler);
6427 octave_unused_parameter (calc_cond);
6429 (*current_liboctave_error_handler)
6430 (
"support for UMFPACK was unavailable or disabled " 6431 "when liboctave was built");
6435 (*current_liboctave_error_handler) (
"incorrect matrix type");
6444 solve_singularity_handler sing_handler,
6445 bool calc_cond)
const 6453 if (nr != nc || nr !=
b.rows ())
6455 (
"matrix dimension mismatch solution of linear equations");
6457 if (nr == 0 ||
b.cols () == 0)
6462 volatile int typ = mattype.
type ();
6467 #if defined (HAVE_CHOLMOD) 6468 cholmod_common Common;
6469 cholmod_common *cm = &Common;
6473 cm->prefer_zomplex =
false;
6479 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
nullptr);
6483 cm->print =
static_cast<int> (spu) + 2;
6484 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &
SparseCholPrint);
6488 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6489 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6491 cm->final_ll =
true;
6493 cholmod_sparse Astore;
6494 cholmod_sparse *
A = &Astore;
6505 #if defined (OCTAVE_ENABLE_64) 6506 A->itype = CHOLMOD_LONG;
6508 A->itype = CHOLMOD_INT;
6510 A->dtype = CHOLMOD_DOUBLE;
6512 A->xtype = CHOLMOD_COMPLEX;
6515 if (
A->x ==
nullptr)
6518 cholmod_sparse Bstore;
6519 cholmod_sparse *
B = &Bstore;
6520 B->nrow =
b.rows ();
6521 B->ncol =
b.cols ();
6524 B->nzmax =
b.nnz ();
6528 #if defined (OCTAVE_ENABLE_64) 6529 B->itype = CHOLMOD_LONG;
6531 B->itype = CHOLMOD_INT;
6533 B->dtype = CHOLMOD_DOUBLE;
6535 B->xtype = CHOLMOD_COMPLEX;
6538 if (
B->x ==
nullptr)
6559 volatile double rcond_plus_one = rcond + 1.0;
6567 sing_handler (rcond);
6582 (static_cast<octave_idx_type>(X->nrow),
6583 static_cast<octave_idx_type>(X->ncol),
6584 static_cast<octave_idx_type>(X->nzmax));
6586 j <= static_cast<octave_idx_type>(X->ncol); j++)
6589 j < static_cast<octave_idx_type>(X->nzmax); j++)
6599 static char blank_name[] =
" ";
6604 (*current_liboctave_warning_with_id_handler)
6605 (
"Octave:missing-dependency",
6606 "support for CHOLMOD was unavailable or disabled " 6607 "when liboctave was built");
6616 #if defined (HAVE_UMFPACK) 6619 sing_handler, calc_cond);
6624 Control (UMFPACK_IRSTEP) = 1;
6653 reinterpret_cast<const double *
> (Ax),
6655 reinterpret_cast<double *> (Xx),
6657 reinterpret_cast<double *
> (Bx),
6658 nullptr, Numeric, control, info);
6665 (*current_liboctave_error_handler)
6666 (
"SparseComplexMatrix::solve solve failed");
6683 sz = x_nz + (
sz > 100 ?
sz : 100);
6694 retval.maybe_compress ();
6696 rcond = Info (UMFPACK_RCOND);
6697 volatile double rcond_plus_one = rcond + 1.0;
6699 if (status == UMFPACK_WARNING_singular_matrix
6705 sing_handler (rcond);
6718 octave_unused_parameter (rcond);
6719 octave_unused_parameter (sing_handler);
6720 octave_unused_parameter (calc_cond);
6722 (*current_liboctave_error_handler)
6723 (
"support for UMFPACK was unavailable or disabled " 6724 "when liboctave was built");
6728 (*current_liboctave_error_handler) (
"incorrect matrix type");
6739 return solve (mattype,
b, info, rcond,
nullptr);
6747 return solve (mattype,
b, info, rcond,
nullptr);
6754 return solve (mattype,
b, info, rcond,
nullptr);
6760 solve_singularity_handler sing_handler,
6761 bool singular_fallback)
const 6764 int typ = mattype.
type (
false);
6767 typ = mattype.
type (*
this);
6783 (*current_liboctave_error_handler) (
"unknown matrix type");
6788 #if defined (USE_QRSOLVE) 6804 return solve (mattype,
b, info, rcond,
nullptr);
6812 return solve (mattype,
b, info, rcond,
nullptr);
6819 return solve (mattype,
b, info, rcond,
nullptr);
6825 solve_singularity_handler sing_handler,
6826 bool singular_fallback)
const 6829 int typ = mattype.
type (
false);
6832 typ = mattype.
type (*
this);
6848 (*current_liboctave_error_handler) (
"unknown matrix type");
6853 #if defined (USE_QRSOLVE) 6869 return solve (mattype,
b, info, rcond,
nullptr);
6877 return solve (mattype,
b, info, rcond,
nullptr);
6884 return solve (mattype,
b, info, rcond,
nullptr);
6890 solve_singularity_handler sing_handler,
6891 bool singular_fallback)
const 6894 int typ = mattype.
type (
false);
6897 typ = mattype.
type (*
this);
6913 (*current_liboctave_error_handler) (
"unknown matrix type");
6918 #if defined (USE_QRSOLVE) 6935 return solve (mattype,
b, info, rcond,
nullptr);
6943 return solve (mattype,
b, info, rcond,
nullptr);
6950 return solve (mattype,
b, info, rcond,
nullptr);
6956 solve_singularity_handler sing_handler,
6957 bool singular_fallback)
const 6960 int typ = mattype.
type (
false);
6963 typ = mattype.
type (*
this);
6979 (*current_liboctave_error_handler) (
"unknown matrix type");
6984 #if defined (USE_QRSOLVE) 6999 return solve (mattype,
b, info, rcond);
7007 return solve (mattype,
b, info, rcond);
7014 return solve (mattype,
b, info, rcond,
nullptr);
7020 solve_singularity_handler sing_handler)
const 7023 return solve (mattype,
tmp, info, rcond,
7024 sing_handler).
column (static_cast<octave_idx_type> (0));
7033 return solve (mattype,
b, info, rcond,
nullptr);
7041 return solve (mattype,
b, info, rcond,
nullptr);
7048 return solve (mattype,
b, info, rcond,
nullptr);
7054 solve_singularity_handler sing_handler)
const 7057 return solve (mattype,
tmp, info, rcond,
7058 sing_handler).
column (static_cast<octave_idx_type> (0));
7066 return solve (
b, info, rcond,
nullptr);
7073 return solve (
b, info, rcond,
nullptr);
7078 double& rcond)
const 7080 return solve (
b, info, rcond,
nullptr);
7086 solve_singularity_handler sing_handler)
const 7089 return solve (mattype,
b,
err, rcond, sing_handler);
7097 return solve (
b, info, rcond,
nullptr);
7105 return solve (
b, info, rcond,
nullptr);
7112 return solve (
b, info, rcond,
nullptr);
7118 solve_singularity_handler sing_handler)
const 7121 return solve (mattype,
b,
err, rcond, sing_handler);
7129 return solve (
b, info, rcond,
nullptr);
7136 return solve (
b, info, rcond,
nullptr);
7142 solve_singularity_handler sing_handler)
const 7145 return solve (mattype,
b,
err, rcond, sing_handler);
7153 return solve (
b, info, rcond,
nullptr);
7161 return solve (
b, info, rcond,
nullptr);
7168 return solve (
b, info, rcond,
nullptr);
7174 solve_singularity_handler sing_handler)
const 7177 return solve (mattype,
b,
err, rcond, sing_handler);
7184 return solve (
b, info, rcond);
7191 return solve (
b, info, rcond);
7196 double& rcond)
const 7198 return solve (
b, info, rcond,
nullptr);
7204 solve_singularity_handler sing_handler)
const 7208 sing_handler).
column (static_cast<octave_idx_type> (0));
7216 return solve (
b, info, rcond,
nullptr);
7224 return solve (
b, info, rcond,
nullptr);
7229 double& rcond)
const 7231 return solve (
b, info, rcond,
nullptr);
7237 solve_singularity_handler sing_handler)
const 7241 sing_handler).
column (static_cast<octave_idx_type> (0));
7362 double r_val =
val.real ();
7363 double i_val =
val.imag ();
7365 if (r_val > max_val)
7368 if (i_val > max_val)
7371 if (r_val < min_val)
7374 if (i_val < min_val)
7420 if ((
rows () == 1 && dim == -1) || dim == 1)
7425 (
cidx (j+1) -
cidx (j) < nr ? 0.0 : 1.0), 1.0);
7439 Complex d = data (i); \ 7440 tmp[ridx (i)] += d * conj (d) 7443 Complex d = data (i); \ 7444 tmp[j] += d * conj (d) 7490 os <<
a.ridx (
i) + 1 <<
' ' << j + 1 <<
' ';
7504 return read_sparse_matrix<elt_type> (
is,
a, octave_read_value<Complex>);
7589 return do_mul_dm_sm<SparseComplexMatrix> (
d,
a);
7594 return do_mul_sm_dm<SparseComplexMatrix> (
a,
d);
7600 return do_mul_dm_sm<SparseComplexMatrix> (
d,
a);
7605 return do_mul_sm_dm<SparseComplexMatrix> (
a,
d);
7611 return do_mul_dm_sm<SparseComplexMatrix> (
d,
a);
7616 return do_mul_sm_dm<SparseComplexMatrix> (
a,
d);
7622 return do_add_dm_sm<SparseComplexMatrix> (
d,
a);
7627 return do_add_dm_sm<SparseComplexMatrix> (
d,
a);
7632 return do_add_dm_sm<SparseComplexMatrix> (
d,
a);
7637 return do_add_sm_dm<SparseComplexMatrix> (
a,
d);
7642 return do_add_sm_dm<SparseComplexMatrix> (
a,
d);
7647 return do_add_sm_dm<SparseComplexMatrix> (
a,
d);
7653 return do_sub_dm_sm<SparseComplexMatrix> (
d,
a);
7658 return do_sub_dm_sm<SparseComplexMatrix> (
d,
a);
7663 return do_sub_dm_sm<SparseComplexMatrix> (
d,
a);
7668 return do_sub_sm_dm<SparseComplexMatrix> (
a,
d);
7673 return do_sub_sm_dm<SparseComplexMatrix> (
a,
d);
7678 return do_sub_sm_dm<SparseComplexMatrix> (
a,
d);
7697 #define EMPTY_RETURN_CHECK(T) \ 7698 if (nr == 0 || nc == 0) \ 7751 bool ja_lt_max = ja < ja_max;
7755 bool jb_lt_max = jb < jb_max;
7757 while (ja_lt_max || jb_lt_max)
7760 if ((! jb_lt_max) || (ja_lt_max && (
a.ridx (ja) <
b.ridx (jb))))
7765 r.
ridx (jx) =
a.ridx (ja);
7770 ja_lt_max= ja < ja_max;
7772 else if ((! ja_lt_max)
7773 || (jb_lt_max && (
b.ridx (jb) <
a.ridx (ja))))
7778 r.
ridx (jx) =
b.ridx (jb);
7783 jb_lt_max= jb < jb_max;
7791 r.
ridx (jx) =
a.ridx (ja);
7795 ja_lt_max= ja < ja_max;
7797 jb_lt_max= jb < jb_max;
7868 bool ja_lt_max = ja < ja_max;
7872 bool jb_lt_max = jb < jb_max;
7874 while (ja_lt_max || jb_lt_max)
7877 if ((! jb_lt_max) || (ja_lt_max && (
a.ridx (ja) <
b.ridx (jb))))
7882 r.
ridx (jx) =
a.ridx (ja);
7887 ja_lt_max= ja < ja_max;
7889 else if ((! ja_lt_max)
7890 || (jb_lt_max && (
b.ridx (jb) <
a.ridx (ja))))
7895 r.
ridx (jx) =
b.ridx (jb);
7900 jb_lt_max= jb < jb_max;
7908 r.
ridx (jx) =
a.ridx (ja);
7912 ja_lt_max= ja < ja_max;
7914 jb_lt_max= jb < jb_max;
friend SparseComplexMatrix conj(const SparseComplexMatrix &a)
octave_idx_type rows(void) const
void resize(octave_idx_type r, octave_idx_type c)
bool all_elements_are_real(void) const
SparseComplexMatrix max(int dim=-1) const
is already an absolute the name is checked against the file system instead of Octave s loadpath In this if otherwise an empty string is returned If the first argument is a cell array of 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
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
std::istream & operator>>(std::istream &is, SparseComplexMatrix &a)
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
MSparse< T > & insert(const Sparse< T > &a, octave_idx_type r, octave_idx_type c)
#define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN)
const octave_base_value const Array< octave_idx_type > & ra_idx
SparseComplexMatrix transpose(void) const
SparseComplexMatrix max(const Complex &c, const SparseComplexMatrix &m)
SparseBoolMatrix any(int dim=-1) const
RowVector row(octave_idx_type i) const
template SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseComplexMatrix, SparseMatrix >(const SparseComplexMatrix &, const SparseMatrix &, octave_idx_type &)
#define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, INIT_VAL, MT_RESULT)
SparseComplexMatrix squeeze(void) const
ComplexMatrix fsolve(MatrixType &mattype, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
std::ostream & operator<<(std::ostream &os, const SparseComplexMatrix &a)
void * factorize(octave_idx_type &err, double &rcond, Matrix &Control, Matrix &Info, solve_singularity_handler sing_handler, bool calc_cond) const
SparseComplexMatrix ipermute(const Array< octave_idx_type > &vec) const
identity matrix If supplied two scalar respectively For allows like xample val
void SparseCholError(int status, char *file, int line, char *message)
bool ishermitian(void) const
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
octave_idx_type columns(void) const
#define F77_DBLE_CMPLX_ARG(x)
bool ishermitian(void) const
SM octinternal_do_mul_pm_sm(const PermMatrix &p, const SM &a)
octave_idx_type * cidx(void)
const T * fortran_vec(void) const
octave_idx_type * triangular_perm(void) const
#define SPARSE_SMSM_BOOL_OPS(M1, M2, ZERO)
void octave_write_complex(std::ostream &os, const Complex &c)
SparseMatrix transpose(void) const
SparseMatrix abs(void) 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)
bool test_any(F fcn) const
SparseComplexMatrix operator*(const SparseComplexMatrix &m, const SparseMatrix &a)
Matrix sum(int dim=-1) const
T & elem(octave_idx_type n)
ComplexMatrix matrix_value(void) const
#define SPARSE_REDUCTION_OP(RET_TYPE, EL_TYPE, OP, INIT_VAL, MT_RESULT)
MSparse< T > diag(octave_idx_type k=0) const
#define FULL_SPARSE_MUL(RET_TYPE, EL_TYPE, ZERO)
SparseMatrix Pr(void) const
nd example oindent opens the file binary numeric values will be read assuming they are stored in IEEE format with the least significant bit and then converted to the native representation Opening a file that is already open simply opens it again and returns a separate file id It is not an error to open a file several though writing to the same file through several different file ids may produce unexpected results The possible values of text mode reading and writing automatically converts linefeeds to the appropriate line end character for the you may append a you must also open the file in binary mode The parameter conversions are currently only supported for and permissions will be set to and then everything is written in a single operation This is very efficient and improves performance c
bool any_element_is_nan(void) const
#define UMFPACK_ZNAME(name)
#define F77_XFCN(f, F, args)
SparseComplexMatrix tinverse(MatrixType &mattype, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
void err_nan_to_logical_conversion(void)
static double get_key(const std::string &key)
MSparse< T > reshape(const dim_vector &new_dims) const
#define SPARSE_SMS_CMP_OPS(M, MZ, CM, S, SZ, CS)
#define SPARSE_FULL_MUL(RET_TYPE, EL_TYPE, ZERO)
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 const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
ComplexMatrix mul_trans(const ComplexMatrix &m, const SparseComplexMatrix &a)
SparseBoolMatrix all(int dim=-1) const
octave_idx_type cols(void) const
octave_value resize(const dim_vector &dv, bool fill=false) const
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
SparseComplexMatrix cumprod(int dim=-1) const
#define SPARSE_ALL_OP(DIM)
SparseComplexMatrix(void)
bool mx_inline_all_real(size_t n, const std::complex< T > *x)
octave_idx_type nnz(void) const
Actual number of nonzero terms.
ComplexColumnVector column(octave_idx_type i) const
SparseComplexMatrix operator-(const ComplexDiagMatrix &d, const SparseMatrix &a)
octave_idx_type nnz(void) const
Count nonzero elements.
#define SPARSE_ANY_OP(DIM)
static const Complex Complex_NaN_result(octave::numeric_limits< double >::NaN(), octave::numeric_limits< double >::NaN())
MSparse< T > squeeze(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
SparseComplexMatrix sum(int dim=-1) const
SparseComplexMatrix permute(const Array< octave_idx_type > &vec, bool inv=false) const
F77_RET_T const F77_INT F77_CMPLX * A
#define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN)
Complex & elem(octave_idx_type n)
#define FULL_SPARSE_MUL_TRANS(RET_TYPE, EL_TYPE, ZERO, CONJ_OP)
#define CHOLMOD_NAME(name)
SparseMatrix Pc(void) const
ComplexMatrix herm_mul(const SparseComplexMatrix &m, const ComplexMatrix &a)
#define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE)
SparseComplexMatrix dinverse(MatrixType &mattype, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
int first_non_singleton(int def=0) const
SparseBoolMatrix operator!(void) const
SparseComplexMatrix hermitian(void) const
void resize(const dim_vector &dv, const T &rfv)
Resizing (with fill).
SparseComplexMatrix diag(octave_idx_type k=0) const
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)
#define EMPTY_RETURN_CHECK(T)
void mark_as_unsymmetric(void)
octave_idx_type * ridx(void)
SparseComplexMatrix min(const Complex &c, const SparseComplexMatrix &m)
SparseComplexMatrix cumsum(int dim=-1) const
bool operator!=(const SparseComplexMatrix &a) const
SparseComplexMatrix & insert(const SparseComplexMatrix &a, octave_idx_type r, octave_idx_type c)
Array< T > array_value(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
ComplexMatrix utsolve(MatrixType &mattype, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
void mark_as_rectangular(void)
With real return the complex result
#define SPARSE_SMSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2)
SparseComplexMatrix reshape(const dim_vector &new_dims) const
ComplexMatrix trisolve(MatrixType &mattype, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
friend class ComplexMatrix
MSparse< T > permute(const Array< octave_idx_type > &vec, bool inv=false) const
ComplexMatrix dsolve(MatrixType &mattype, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
SparseComplexMatrix prod(int dim=-1) const
octave_idx_type cols(void) const
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
SparseComplexMatrix operator+(const ComplexDiagMatrix &d, const SparseMatrix &a)
int SparseCholPrint(const char *fmt,...)
bool is_dense(void) const
#define END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
octave_f77_int_type F77_INT
ComplexMatrix trans_mul(const SparseComplexMatrix &m, const ComplexMatrix &a)
#define SPARSE_SSM_CMP_OPS(S, SZ, SC, M, MZ, MC)
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$)
bool all_integers(double &max_val, double &min_val) const
bool any_element_is_inf_or_nan(void) const
ComplexMatrix solve(MatrixType &mattype, const Matrix &b) const
ComplexDET determinant(void) const
SparseComplexMatrix conj(const SparseComplexMatrix &a)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Sparse< T > maybe_compress(bool remove_zeros=false)
SparseComplexMatrix inverse(void) const
bool xtoo_large_for_float(double x)
MSparse< T > ipermute(const Array< octave_idx_type > &vec) const
ComplexColumnVector column(octave_idx_type i) const
ComplexRowVector row(octave_idx_type i) const
dim_vector dims(void) const
MatrixType transpose(void) const
SparseComplexMatrix min(int dim=-1) const
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)
#define BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
octave_idx_type ndims(void) const
Number of dimensions.
SparseMatrix Q(void) const
std::complex< double > Complex
bool operator==(const SparseComplexMatrix &a) const
#define SPARSE_SSM_BOOL_OPS(S, M, ZERO)
ColumnVector real(const ComplexColumnVector &a)
write the output to stdout if nargout is
Vector representing the dimensions (size) of an Array.
void warn_singular_matrix(double rcond)
template ComplexMatrix dmsolve< ComplexMatrix, SparseComplexMatrix, Matrix >(const SparseComplexMatrix &, const Matrix &, octave_idx_type &)
ComplexMatrix ltsolve(MatrixType &mattype, const Matrix &b, octave_idx_type &info, double &rcond, solve_singularity_handler sing_handler, bool calc_cond=false) const
RT dmsolve(const ST &a, const T &b, octave_idx_type &info)
octave_idx_type rows(void) const
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
SparseComplexMatrix sumsq(int dim=-1) const
bool too_large_for_float(void) const