7SEXP
mkDet(
double modulus,
int logarithm,
int sign)
9 SEXP nms = PROTECT(allocVector(STRSXP, 2)),
10 cl = PROTECT(mkString(
"det")),
11 det = PROTECT(allocVector(VECSXP, 2)),
12 det0 = PROTECT(ScalarReal((logarithm) ? modulus : exp(modulus))),
13 det1 = PROTECT(ScalarInteger(sign)),
14 det0a = PROTECT(ScalarLogical(logarithm));
15 SET_STRING_ELT(nms, 0, mkChar(
"modulus"));
16 SET_STRING_ELT(nms, 1, mkChar(
"sign"));
17 setAttrib(det, R_NamesSymbol, nms);
18 setAttrib(det, R_ClassSymbol,
cl);
19 setAttrib(det0, install(
"logarithm"), det0a);
20 SET_VECTOR_ELT(det, 0, det0);
21 SET_VECTOR_ELT(det, 1, det1);
29#define DETERMINANT_START \
30 SEXP dim = GET_SLOT(obj, Matrix_DimSym); \
31 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1]; \
33 error(_("determinant of non-square matrix is undefined")); \
34 int givelog = asLogical(logarithm) != 0; \
40 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
44 R_xlen_t n1a = (R_xlen_t) n + 1;
45 if (TYPEOF(x) == CPLXSXP) {
46 Rcomplex *px = COMPLEX(x);
47 for (j = 0; j < n; ++j) {
48 modulus += log(hypot((*px).r, (*px).i));
53 int *ppivot = INTEGER(pivot);
55 for (j = 0; j < n; ++j) {
56 if (ISNAN(*px) || *px >= 0.0) {
61 modulus += log(-(*px));
72 return mkDet(modulus, givelog, sign);
80 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
84 int *ppivot = INTEGER(pivot);
87 char ul = *CHAR(STRING_ELT(uplo, 0));
90 XLENGTH(x) == (R_xlen_t) m * m;
91 R_xlen_t n1a = (R_xlen_t) n + 1;
92 if (TYPEOF(x) == CPLXSXP) {
93 Rcomplex *px = COMPLEX(x), a, b,
c;
96 modulus += log(hypot((*px).r, (*px).i));
97 px += (unpacked) ? n1a : ((ul ==
'U') ? j + 2 : n - j);
102 px += (unpacked) ? n1a : j + 2;
105 px += (unpacked) ? n1a : j + 3;
108 px += (unpacked) ? n1a : n - j;
110 px += (unpacked) ? n1a : n - j - 1;
112 modulus += log(hypot(a.r * b.r - a.i * b.i -
113 c.r *
c.r +
c.i *
c.i,
114 a.r * b.i + a.i * b.r -
120 double *px = REAL(x), a, b,
c, logab, logcc;
126 modulus += log(-(*px));
129 px += (unpacked) ? n1a : ((ul ==
'U') ? j + 2 : n - j);
134 px += (unpacked) ? n1a : j + 2;
137 px += (unpacked) ? n1a : j + 3;
140 px += (unpacked) ? n1a : n - j;
142 px += (unpacked) ? n1a : n - j - 1;
144 logab = log(fabs(a)) + log(fabs(b));
145 logcc = 2.0 * log(fabs(
c));
146 if ((a < 0.0) != (b < 0.0)) {
148 modulus += logspace_add(logab, logcc);
150 }
else if (logab < logcc) {
152 modulus += logspace_sub(logcc, logab);
156 modulus += logspace_sub(logab, logcc);
165 return mkDet(modulus, givelog, sign);
173 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
177 char ul = *CHAR(STRING_ELT(uplo, 0));
180 XLENGTH(x) == (R_xlen_t) m * m;
181 R_xlen_t n1a = (R_xlen_t) n + 1;
182 if (TYPEOF(x) == CPLXSXP) {
183 Rcomplex *px = COMPLEX(x);
184 for (j = 0; j < n; ++j) {
185 modulus += log(hypot((*px).r, (*px).i));
186 px += (unpacked) ? n1a : ((ul ==
'U') ? j + 2 : n - j);
189 double *px = REAL(x);
190 for (j = 0; j < n; ++j) {
191 if (ISNAN(*px) || *px >= 0.0)
194 modulus += log(-(*px));
197 px += (unpacked) ? n1a : ((ul ==
'U') ? j + 2 : n - j);
204 return mkDet(modulus, givelog, sign);
213 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
218 int *pp = INTEGER(p) + 1, *pi = INTEGER(i), j, k = 0, kend;
219 if (TYPEOF(x) == CPLXSXP) {
220 Rcomplex *px = COMPLEX(x);
221 for (j = 0; j < n; ++j) {
223 if (k < kend && pi[kend - 1] == j)
224 modulus += log(hypot(px[kend - 1].r, px[kend - 1].i));
227 return mkDet(R_NegInf, givelog, 1);
232 double *px = REAL(x);
233 for (j = 0; j < n; ++j) {
235 if (k < kend && pi[kend - 1] == j) {
236 if (ISNAN(px[kend - 1]) || px[kend - 1] >= 0.0)
237 modulus += log(px[kend - 1]);
239 modulus += log(-px[kend - 1]);
244 return mkDet(R_NegInf, givelog, 1);
250 int signPerm(
const int *,
int,
int);
253 if (
signPerm(INTEGER(p), LENGTH(p), 0) < 0)
256 if (
signPerm(INTEGER(p), LENGTH(p), 0) < 0)
263 return mkDet(modulus, givelog, sign);
272 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
275 if (INTEGER(dim)[0] > n)
276 error(
_(
"%s(<%s>) does not support structurally rank deficient case"),
277 "determinant",
"sparseQR");
282 int *pp = INTEGER(p) + 1, *pi = INTEGER(i), j, k = 0, kend;
283 if (TYPEOF(x) == CPLXSXP) {
284 Rcomplex *px = COMPLEX(x);
285 for (j = 0; j < n; ++j) {
287 if (k < kend && pi[kend - 1] == j)
288 modulus += log(hypot(px[kend - 1].r, px[kend - 1].i));
291 return mkDet(R_NegInf, givelog, 1);
296 double *px = REAL(x);
297 for (j = 0; j < n; ++j) {
299 if (k < kend && pi[kend - 1] == j) {
300 if (ISNAN(px[kend - 1]) || px[kend - 1] >= 0.0)
301 modulus += log(px[kend - 1]);
303 modulus += log(-px[kend - 1]);
308 return mkDet(R_NegInf, givelog, 1);
314 int signPerm(
const int *,
int,
int);
317 if (
signPerm(INTEGER(p), LENGTH(p), 0) < 0)
320 if (
signPerm(INTEGER(p), LENGTH(p), 0) < 0)
329 return mkDet(modulus, givelog, sign);
336 cholmod_factor *L =
M2CHF(obj, 1);
337 int sign = (L->xtype == CHOLMOD_COMPLEX) ? NA_INTEGER : 1;
340 int j, sqrt_ = asLogical(sqrt);
343 nsuper = (int) L->nsuper,
344 *psuper = (
int *) L->super,
345 *ppi = (
int *) L->pi,
346 *ppx = (
int *) L->px;
348 if (L->xtype == CHOLMOD_COMPLEX) {
349 Rcomplex *px = (Rcomplex *) L->x, *px_;
350 for (k = 0; k < nsuper; ++k) {
351 nc = psuper[k + 1] - psuper[k];
352 nr1a = (R_xlen_t) (ppi[k + 1] - ppi[k]) + 1;
354 for (j = 0; j < nc; ++j) {
355 modulus += log(hypot((*px_).r, (*px_).i));
360 double *px = (
double *) L->x, *px_;
361 for (k = 0; k < nsuper; ++k) {
362 nc = psuper[k + 1] - psuper[k];
363 nr1a = (R_xlen_t) (ppi[k + 1] - ppi[k]) + 1;
365 for (j = 0; j < nc; ++j) {
366 modulus += log(*px_);
373 int *pp = (
int *) L->p;
374 if (L->xtype == CHOLMOD_COMPLEX) {
375 Rcomplex *px = (Rcomplex *) L->x;
376 for (j = 0; j < n; ++j)
377 modulus += log(hypot(px[pp[j]].r, px[pp[j]].i));
381 double *px = (
double *) L->x;
383 for (j = 0; j < n; ++j)
384 modulus += log(px[pp[j]]);
387 for (j = 0; j < n; ++j) {
388 if (ISNAN(px[pp[j]]) || px[pp[j]] >= 0.0)
389 modulus += log(px[pp[j]]);
392 return mkDet(R_NaN, givelog, 1);
393 modulus += log(-px[pp[j]]);
404 return mkDet(modulus, givelog, sign);