Matrix r5059
Loading...
Searching...
No Matches
determinant.c
Go to the documentation of this file.
1/* C implementation of methods for determinant */
2
3#include "cholmod-etc.h"
4#include "Mdefines.h"
5#include <Rmath.h> /* logspace_add, logspace_sub */
6
7static
8SEXP det(double modulus, int logarithm, int sign)
9{
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);
21 Rf_setAttrib(s_modulus, Matrix_logarithmSym, s_logarithm);
22 SET_VECTOR_ELT(ans, 0, s_modulus);
23 SET_VECTOR_ELT(ans, 1, s_sign);
24 UNPROTECT(6);
25 return ans;
26}
27
28SEXP denseLU_determinant(SEXP s_trf, SEXP s_logarithm)
29{
30
31#define DETERMINANT_START(_F_) \
32 int *pdim = DIM(_F_), n = pdim[1]; \
33 if (pdim[0] != n) \
34 Rf_error(_("matrix is not square")); \
35 int givelog = Rf_asLogical(s_logarithm); \
36 double modulus = 0.0; /* result for n == 0 */
37
38 DETERMINANT_START(s_trf);
39
40 SEXP x = PROTECT(GET_SLOT(s_trf, Matrix_xSym));
41 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
42
43 if (n > 0) {
44 int j;
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));
50 px += n1a;
51 }
52 } else {
53 SEXP pivot = GET_SLOT(s_trf, Matrix_permSym);
54 int *ppivot = INTEGER(pivot);
55 double *px = REAL(x);
56 for (j = 0; j < n; ++j) {
57 if (ISNAN(*px) || *px >= 0.0) {
58 modulus += log(*px);
59 if (*ppivot != j + 1)
60 sign = -sign;
61 } else {
62 modulus += log(-(*px));
63 if (*ppivot == j + 1)
64 sign = -sign;
65 }
66 px += n1a;
67 ppivot += 1;
68 }
69 }
70 }
71
72 UNPROTECT(1); /* x */
73 return det(modulus, givelog, sign);
74}
75
76SEXP denseBunchKaufman_determinant(SEXP s_trf, SEXP s_logarithm)
77{
78 DETERMINANT_START(s_trf);
79
80 SEXP x = PROTECT(GET_SLOT(s_trf, Matrix_xSym));
81 int sign = 1;
82
83 char ct = (TYPEOF(x) == CPLXSXP) ? TRANS(s_trf) : 'C';
84 if (ct != 'C')
85 sign = NA_INTEGER;
86
87 if (n > 0) {
88 char ul = UPLO(s_trf);
89
90 SEXP pivot = GET_SLOT(s_trf, Matrix_permSym);
91 int *ppivot = INTEGER(pivot);
92
93 int j = 0, packed = XLENGTH(x) != (int_fast64_t) n * n;
94 R_xlen_t n1a = (R_xlen_t) n + 1;
95 double logab, logcc;
96 if (TYPEOF(x) == CPLXSXP) {
97 Rcomplex *px = COMPLEX(x), a, b, c;
98 while (j < n) {
99 a = *px;
100 if (ppivot[j] > 0) {
101 if (ct != 'C')
102 modulus += log(hypot(a.r, a.i));
103 else if (ISNAN(a.r) || a.r >= 0.0)
104 modulus += log(a.r);
105 else {
106 modulus += log(-a.r);
107 sign = -sign;
108 }
109 px += (!packed) ? n1a : ((ul == 'U') ? j + 2 : n - j);
110 j += 1;
111 } else {
112 if (ul == 'U') {
113 px += (!packed) ? n1a : j + 2;
114 b = *px;
115 c = *(px - 1);
116 px += (!packed) ? n1a : j + 3;
117 } else {
118 c = *(px + 1);
119 px += (!packed) ? n1a : n - j;
120 b = *px;
121 px += (!packed) ? n1a : n - j - 1;
122 }
123 if (ct != 'C')
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 -
127 2.0 * c.r * c.i));
128 else {
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)) {
132 /* det = ab - cc = -(abs(ab) + cc) < 0 */
133 modulus += Rf_logspace_add(logab, logcc);
134 sign = -sign;
135 } else if (logab < logcc) {
136 /* det = ab - cc = -(cc - ab) < 0 */
137 modulus += Rf_logspace_sub(logcc, logab);
138 sign = -sign;
139 } else {
140 /* det = ab - cc > 0 */
141 modulus += Rf_logspace_sub(logab, logcc);
142 }
143 }
144 j += 2;
145 }
146 }
147 } else {
148 double *px = REAL(x), a, b, c;
149 while (j < n) {
150 a = *px;
151 if (ppivot[j] > 0) {
152 if (*px >= 0.0)
153 modulus += log(a);
154 else {
155 modulus += log(-a);
156 sign = -sign;
157 }
158 px += (!packed) ? n1a : ((ul == 'U') ? j + 2 : n - j);
159 j += 1;
160 } else {
161 if (ul == 'U') {
162 px += (!packed) ? n1a : j + 2;
163 b = *px;
164 c = *(px - 1);
165 px += (!packed) ? n1a : j + 3;
166 } else {
167 c = *(px + 1);
168 px += (!packed) ? n1a : n - j;
169 b = *px;
170 px += (!packed) ? n1a : n - j - 1;
171 }
172 logab = log(fabs(a)) + log(fabs(b));
173 logcc = 2.0 * log(fabs(c));
174 if ((a < 0.0) != (b < 0.0)) {
175 /* det = ab - cc = -(abs(ab) + cc) < 0 */
176 modulus += Rf_logspace_add(logab, logcc);
177 sign = -sign;
178 } else if (logab < logcc) {
179 /* det = ab - cc = -(cc - ab) < 0 */
180 modulus += Rf_logspace_sub(logcc, logab);
181 sign = -sign;
182 } else {
183 /* det = ab - cc > 0 */
184 modulus += Rf_logspace_sub(logab, logcc);
185 }
186 j += 2;
187 }
188 }
189 }
190 }
191
192 UNPROTECT(1); /* x */
193 return det(modulus, givelog, sign);
194}
195
196SEXP denseCholesky_determinant(SEXP s_trf, SEXP s_logarithm)
197{
198 DETERMINANT_START(s_trf);
199
200 SEXP x = PROTECT(GET_SLOT(s_trf, Matrix_xSym));
201 int sign = 1;
202
203 if (n > 0) {
204 char ul = UPLO(s_trf);
205
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);
213 }
214 } else {
215 double *px = REAL(x);
216 for (j = 0; j < n; ++j) {
217 modulus += log(*px);
218 px += (!packed) ? n1a : ((ul == 'U') ? j + 2 : n - j);
219 }
220 }
221 modulus *= 2.0;
222 }
223
224 UNPROTECT(1); /* x */
225 return det(modulus, givelog, sign);
226}
227
228SEXP sparseQR_determinant(SEXP orf, SEXP s_logarithm)
229{
230 DETERMINANT_START(orf);
231
232 SEXP R = PROTECT(GET_SLOT(orf, Matrix_RSym)),
233 x = PROTECT(GET_SLOT(R, Matrix_xSym));
234 int sign = 1;
235
236 if (DIM(R)[0] > n)
237 Rf_error(_("%s(<%s>) does not support structurally rank deficient case"),
238 "determinant", "sparseQR");
239
240 if (n > 0) {
241 SEXP p = PROTECT(GET_SLOT(R, Matrix_pSym)),
242 i = PROTECT(GET_SLOT(R, Matrix_iSym));
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) {
247 kend = pp[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);
251 else {
252 modulus += log(-px[kend - 1].r);
253 sign = -sign;
254 }
255 } else {
256 UNPROTECT(4); /* i, p, x, R */
257 return det(R_NegInf, givelog, 1);
258 }
259 k = kend;
260 }
261 } else {
262 double *px = REAL(x);
263 for (j = 0, k = 0; j < n; ++j) {
264 kend = pp[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]);
268 else {
269 modulus += log(-px[kend - 1]);
270 sign = -sign;
271 }
272 } else {
273 UNPROTECT(4); /* i, p, x, R */
274 return det(R_NegInf, givelog, 1);
275 }
276 k = kend;
277 }
278 }
279 UNPROTECT(2); /* i, p */
280
281 /* defined in ./perm.c : */
282 int signPerm(const int *, int, int);
283
284 p = GET_SLOT(orf, Matrix_pSym);
285 if (signPerm(INTEGER(p), LENGTH(p), 0) < 0)
286 sign = -sign;
287 p = GET_SLOT(orf, Matrix_qSym);
288 if (signPerm(INTEGER(p), LENGTH(p), 0) < 0)
289 sign = -sign;
290 if (n % 2)
291 sign = -sign;
292 }
293
294 UNPROTECT(2); /* x, R */
295 return det(modulus, givelog, sign);
296}
297
298SEXP sparseLU_determinant(SEXP s_trf, SEXP s_logarithm)
299{
300 DETERMINANT_START(s_trf);
301
302 SEXP U = PROTECT(GET_SLOT(s_trf, Matrix_USym)),
303 x = PROTECT(GET_SLOT(U, Matrix_xSym));
304 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
305
306 if (n > 0) {
307 SEXP p = PROTECT(GET_SLOT(U, Matrix_pSym)),
308 i = PROTECT(GET_SLOT(U, Matrix_iSym));
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) {
313 kend = pp[j];
314 if (k < kend && pi[kend - 1] == j)
315 modulus += log(hypot(px[kend - 1].r, px[kend - 1].i));
316 else {
317 UNPROTECT(4); /* i, p, x, U */
318 return det(R_NegInf, givelog, 1);
319 }
320 k = kend;
321 }
322 } else {
323 double *px = REAL(x);
324 for (j = 0, k = 0; j < n; ++j) {
325 kend = pp[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]);
329 else {
330 modulus += log(-px[kend - 1]);
331 sign = -sign;
332 }
333 } else {
334 UNPROTECT(4); /* i, p, x, U */
335 return det(R_NegInf, givelog, 1);
336 }
337 k = kend;
338 }
339 }
340 UNPROTECT(2); /* i, p */
341
342 if (TYPEOF(x) != CPLXSXP) {
343 /* defined in ./perm.c : */
344 int signPerm(const int *, int, int);
345
346 p = GET_SLOT(s_trf, Matrix_pSym);
347 if (signPerm(INTEGER(p), LENGTH(p), 0) < 0)
348 sign = -sign;
349 p = GET_SLOT(s_trf, Matrix_qSym);
350 if (signPerm(INTEGER(p), LENGTH(p), 0) < 0)
351 sign = -sign;
352 }
353 }
354
355 UNPROTECT(2); /* x, U */
356 return det(modulus, givelog, sign);
357}
358
359SEXP sparseCholesky_determinant(SEXP s_trf, SEXP s_logarithm, SEXP s_root)
360{
361 DETERMINANT_START(s_trf);
362
363 cholmod_factor *L = M2CHF(s_trf, 1);
364 int sign = 1;
365
366 if (n > 0) {
367 int j, root = Rf_asLogical(s_root);
368 if (L->is_super) {
369 int k, nc,
370 nsuper = (int) L->nsuper,
371 *psuper = (int *) L->super,
372 *ppi = (int *) L->pi,
373 *ppx = (int *) L->px;
374 R_xlen_t nr1a;
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;
380 py = px + ppx[k];
381 for (j = 0; j < nc; ++j) {
382 modulus += log((*py).r);
383 py += nr1a;
384 }
385 }
386 } else {
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;
391 py = px + ppx[k];
392 for (j = 0; j < nc; ++j) {
393 modulus += log(*py);
394 py += nr1a;
395 }
396 }
397 }
398 modulus *= 2.0;
399 } else {
400 int *pp = (int *) L->p;
401 if (L->xtype == CHOLMOD_COMPLEX) {
402 Rcomplex *px = (Rcomplex *) L->x;
403 if (L->is_ll) {
404 for (j = 0; j < n; ++j)
405 modulus += log(px[pp[j]].r);
406 modulus *= 2.0;
407 } else {
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);
411 else {
412 if (root)
413 return det(R_NaN, givelog, 1);
414 modulus += log(-px[pp[j]].r);
415 sign = -sign;
416 }
417 }
418 }
419 } else {
420 double *px = (double *) L->x;
421 if (L->is_ll) {
422 for (j = 0; j < n; ++j)
423 modulus += log(px[pp[j]]);
424 modulus *= 2.0;
425 } else {
426 for (j = 0; j < n; ++j) {
427 if (ISNAN(px[pp[j]]) || px[pp[j]] >= 0.0)
428 modulus += log(px[pp[j]]);
429 else {
430 if (root)
431 return det(R_NaN, givelog, 1);
432 modulus += log(-px[pp[j]]);
433 sign = -sign;
434 }
435 }
436 }
437 }
438 }
439 if (root)
440 modulus *= 0.5;
441 }
442
443 return det(modulus, givelog, sign);
444}
_
#define _(String)
Definition Mdefines.h:66
UPLO
#define UPLO(x)
Definition Mdefines.h:101
TRANS
#define TRANS(x)
Definition Mdefines.h:106
DIM
#define DIM(x)
Definition Mdefines.h:85
GET_SLOT
#define GET_SLOT(x, name)
Definition Mdefines.h:72
TYPEOF
#define TYPEOF(s)
Definition Mdefines.h:123
M2CHF
cholmod_factor * M2CHF(SEXP obj, int values)
Definition cholmod-etc.c:39
c
cholmod_common c
Definition cholmod-etc.c:5
sparseCholesky_determinant
SEXP sparseCholesky_determinant(SEXP s_trf, SEXP s_logarithm, SEXP s_root)
Definition determinant.c:359
denseCholesky_determinant
SEXP denseCholesky_determinant(SEXP s_trf, SEXP s_logarithm)
Definition determinant.c:196
sparseLU_determinant
SEXP sparseLU_determinant(SEXP s_trf, SEXP s_logarithm)
Definition determinant.c:298
DETERMINANT_START
#define DETERMINANT_START(_F_)
denseLU_determinant
SEXP denseLU_determinant(SEXP s_trf, SEXP s_logarithm)
Definition determinant.c:28
det
static SEXP det(double modulus, int logarithm, int sign)
Definition determinant.c:8
denseBunchKaufman_determinant
SEXP denseBunchKaufman_determinant(SEXP s_trf, SEXP s_logarithm)
Definition determinant.c:76
sparseQR_determinant
SEXP sparseQR_determinant(SEXP orf, SEXP s_logarithm)
Definition determinant.c:228
Matrix_permSym
SEXP Matrix_permSym
Definition init.c:623
Matrix_logarithmSym
SEXP Matrix_logarithmSym
Definition init.c:613
Matrix_RSym
SEXP Matrix_RSym
Definition init.c:600
Matrix_xSym
SEXP Matrix_xSym
Definition init.c:635
Matrix_iSym
SEXP Matrix_iSym
Definition init.c:607
Matrix_qSym
SEXP Matrix_qSym
Definition init.c:627
Matrix_USym
SEXP Matrix_USym
Definition init.c:601
Matrix_pSym
SEXP Matrix_pSym
Definition init.c:622
signPerm
int signPerm(const int *p, int n, int off)
Definition perm.c:22
Generated by doxygen 1.12.0