8SEXP
det(
double modulus,
int logarithm,
int sign)
10 SEXP s_modulus = PROTECT(Rf_ScalarReal((logarithm) ? modulus : exp(modulus))),
11 s_logarithm = PROTECT(Rf_ScalarLogical(logarithm)),
12 s_sign = PROTECT(Rf_ScalarInteger(sign)),
13 ans = PROTECT(Rf_allocVector(VECSXP, 2)),
14 class = PROTECT(Rf_allocVector(STRSXP, 1)),
15 names = PROTECT(Rf_allocVector(STRSXP, 2));
16 SET_STRING_ELT(
class, 0, Rf_mkChar(
"det"));
17 SET_STRING_ELT(names, 0, Rf_mkChar(
"modulus"));
18 SET_STRING_ELT(names, 1, Rf_mkChar(
"sign"));
19 Rf_setAttrib(ans, R_ClassSymbol,
class);
20 Rf_setAttrib(ans, R_NamesSymbol, names);
22 SET_VECTOR_ELT(ans, 0, s_modulus);
23 SET_VECTOR_ELT(ans, 1, s_sign);
31#define DETERMINANT_START(_F_) \
32 int *pdim = DIM(_F_), n = pdim[1]; \
34 Rf_error(_("matrix is not square")); \
35 int givelog = Rf_asLogical(s_logarithm); \
41 int sign = (
TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
45 R_xlen_t n1a = (R_xlen_t) n + 1;
46 if (
TYPEOF(x) == CPLXSXP) {
47 Rcomplex *px = COMPLEX(x);
48 for (j = 0; j < n; ++j) {
49 modulus += log(hypot((*px).r, (*px).i));
54 int *ppivot = INTEGER(pivot);
56 for (j = 0; j < n; ++j) {
57 if (ISNAN(*px) || *px >= 0.0) {
62 modulus += log(-(*px));
73 return det(modulus, givelog, sign);
83 char ct = (
TYPEOF(x) == CPLXSXP) ?
TRANS(s_trf) :
'C';
88 char ul =
UPLO(s_trf);
91 int *ppivot = INTEGER(pivot);
93 int j = 0, packed = XLENGTH(x) != (int_fast64_t) n * n;
94 R_xlen_t n1a = (R_xlen_t) n + 1;
96 if (
TYPEOF(x) == CPLXSXP) {
97 Rcomplex *px = COMPLEX(x), a, b,
c;
102 modulus += log(hypot(a.r, a.i));
103 else if (ISNAN(a.r) || a.r >= 0.0)
106 modulus += log(-a.r);
109 px += (!packed) ? n1a : ((ul ==
'U') ? j + 2 : n - j);
113 px += (!packed) ? n1a : j + 2;
116 px += (!packed) ? n1a : j + 3;
119 px += (!packed) ? n1a : n - j;
121 px += (!packed) ? n1a : n - j - 1;
124 modulus += log(hypot(a.r * b.r - a.i * b.i -
125 c.r *
c.r +
c.i *
c.i,
126 a.r * b.i + a.i * b.r -
129 logab = log(fabs(a.r)) + log(fabs(b.r));
130 logcc = 2.0 * log(hypot(
c.r,
c.i));
131 if ((a.r < 0.0) != (b.r < 0.0)) {
133 modulus += Rf_logspace_add(logab, logcc);
135 }
else if (logab < logcc) {
137 modulus += Rf_logspace_sub(logcc, logab);
141 modulus += Rf_logspace_sub(logab, logcc);
148 double *px = REAL(x), a, b,
c;
158 px += (!packed) ? n1a : ((ul ==
'U') ? j + 2 : n - j);
162 px += (!packed) ? n1a : j + 2;
165 px += (!packed) ? n1a : j + 3;
168 px += (!packed) ? n1a : n - j;
170 px += (!packed) ? n1a : n - j - 1;
172 logab = log(fabs(a)) + log(fabs(b));
173 logcc = 2.0 * log(fabs(
c));
174 if ((a < 0.0) != (b < 0.0)) {
176 modulus += Rf_logspace_add(logab, logcc);
178 }
else if (logab < logcc) {
180 modulus += Rf_logspace_sub(logcc, logab);
184 modulus += Rf_logspace_sub(logab, logcc);
193 return det(modulus, givelog, sign);
204 char ul =
UPLO(s_trf);
206 int j, packed = XLENGTH(x) != (int_fast64_t) n * n;
207 R_xlen_t n1a = (R_xlen_t) n + 1;
208 if (
TYPEOF(x) == CPLXSXP) {
209 Rcomplex *px = COMPLEX(x);
210 for (j = 0; j < n; ++j) {
211 modulus += log((*px).r);
212 px += (!packed) ? n1a : ((ul ==
'U') ? j + 2 : n - j);
215 double *px = REAL(x);
216 for (j = 0; j < n; ++j) {
218 px += (!packed) ? n1a : ((ul ==
'U') ? j + 2 : n - j);
225 return det(modulus, givelog, sign);
237 Rf_error(
_(
"%s(<%s>) does not support structurally rank deficient case"),
238 "determinant",
"sparseQR");
243 int *pp = INTEGER(p) + 1, *pi = INTEGER(i), j, k, kend;
244 if (
TYPEOF(x) == CPLXSXP) {
245 Rcomplex *px = COMPLEX(x);
246 for (j = 0, k = 0; j < n; ++j) {
248 if (k < kend && pi[kend - 1] == j) {
249 if (ISNAN(px[kend - 1].r) || px[kend - 1].r >= 0.0)
250 modulus += log(px[kend - 1].r);
252 modulus += log(-px[kend - 1].r);
257 return det(R_NegInf, givelog, 1);
262 double *px = REAL(x);
263 for (j = 0, k = 0; j < n; ++j) {
265 if (k < kend && pi[kend - 1] == j) {
266 if (ISNAN(px[kend - 1]) || px[kend - 1] >= 0.0)
267 modulus += log(px[kend - 1]);
269 modulus += log(-px[kend - 1]);
274 return det(R_NegInf, givelog, 1);
282 int signPerm(
const int *,
int,
int);
285 if (
signPerm(INTEGER(p), LENGTH(p), 0) < 0)
288 if (
signPerm(INTEGER(p), LENGTH(p), 0) < 0)
295 return det(modulus, givelog, sign);
304 int sign = (
TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
309 int *pp = INTEGER(p) + 1, *pi = INTEGER(i), j, k, kend;
310 if (
TYPEOF(x) == CPLXSXP) {
311 Rcomplex *px = COMPLEX(x);
312 for (j = 0, k = 0; j < n; ++j) {
314 if (k < kend && pi[kend - 1] == j)
315 modulus += log(hypot(px[kend - 1].r, px[kend - 1].i));
318 return det(R_NegInf, givelog, 1);
323 double *px = REAL(x);
324 for (j = 0, k = 0; j < n; ++j) {
326 if (k < kend && pi[kend - 1] == j) {
327 if (ISNAN(px[kend - 1]) || px[kend - 1] >= 0.0)
328 modulus += log(px[kend - 1]);
330 modulus += log(-px[kend - 1]);
335 return det(R_NegInf, givelog, 1);
342 if (
TYPEOF(x) != CPLXSXP) {
344 int signPerm(
const int *,
int,
int);
347 if (
signPerm(INTEGER(p), LENGTH(p), 0) < 0)
350 if (
signPerm(INTEGER(p), LENGTH(p), 0) < 0)
356 return det(modulus, givelog, sign);
363 cholmod_factor *L =
M2CHF(s_trf, 1);
367 int j, root = Rf_asLogical(s_root);
370 nsuper = (int) L->nsuper,
371 *psuper = (
int *) L->super,
372 *ppi = (
int *) L->pi,
373 *ppx = (
int *) L->px;
375 if (L->xtype == CHOLMOD_COMPLEX) {
376 Rcomplex *px = (Rcomplex *) L->x, *py;
377 for (k = 0; k < nsuper; ++k) {
378 nc = psuper[k + 1] - psuper[k];
379 nr1a = (R_xlen_t) (ppi[k + 1] - ppi[k]) + 1;
381 for (j = 0; j < nc; ++j) {
382 modulus += log((*py).r);
387 double *px = (
double *) L->x, *py;
388 for (k = 0; k < nsuper; ++k) {
389 nc = psuper[k + 1] - psuper[k];
390 nr1a = (R_xlen_t) (ppi[k + 1] - ppi[k]) + 1;
392 for (j = 0; j < nc; ++j) {
400 int *pp = (
int *) L->p;
401 if (L->xtype == CHOLMOD_COMPLEX) {
402 Rcomplex *px = (Rcomplex *) L->x;
404 for (j = 0; j < n; ++j)
405 modulus += log(px[pp[j]].r);
408 for (j = 0; j < n; ++j) {
409 if (ISNAN(px[pp[j]].r) || px[pp[j]].r >= 0.0)
410 modulus += log(px[pp[j]].r);
413 return det(R_NaN, givelog, 1);
414 modulus += log(-px[pp[j]].r);
420 double *px = (
double *) L->x;
422 for (j = 0; j < n; ++j)
423 modulus += log(px[pp[j]]);
426 for (j = 0; j < n; ++j) {
427 if (ISNAN(px[pp[j]]) || px[pp[j]] >= 0.0)
428 modulus += log(px[pp[j]]);
431 return det(R_NaN, givelog, 1);
432 modulus += log(-px[pp[j]]);
443 return det(modulus, givelog, sign);