Matrix r4655
Loading...
Searching...
No Matches
determinant.c
Go to the documentation of this file.
1#include <Rmath.h> /* math.h, logspace_add, logspace_sub */
2#include "Mdefines.h"
3#include "cholmod-etc.h"
4#include "determinant.h"
5
6static
7SEXP mkDet(double modulus, int logarithm, int sign)
8{
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);
22 UNPROTECT(6);
23 return det;
24}
25
26SEXP denseLU_determinant(SEXP obj, SEXP logarithm)
27{
28
29#define DETERMINANT_START \
30 SEXP dim = GET_SLOT(obj, Matrix_DimSym); \
31 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1]; \
32 if (m != n) \
33 error(_("determinant of non-square matrix is undefined")); \
34 int givelog = asLogical(logarithm) != 0; \
35 double modulus = 0.0; /* result for n == 0 */
36
37 DETERMINANT_START;
38
39 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym));
40 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
41
42 if (n > 0) {
43 int j;
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));
49 px += n1a;
50 }
51 } else {
52 SEXP pivot = GET_SLOT(obj, Matrix_permSym);
53 int *ppivot = INTEGER(pivot);
54 double *px = REAL(x);
55 for (j = 0; j < n; ++j) {
56 if (ISNAN(*px) || *px >= 0.0) {
57 modulus += log(*px);
58 if (*ppivot != j + 1)
59 sign = -sign;
60 } else {
61 modulus += log(-(*px));
62 if (*ppivot == j + 1)
63 sign = -sign;
64 }
65 px += n1a;
66 ppivot += 1;
67 }
68 }
69 }
70
71 UNPROTECT(1); /* x */
72 return mkDet(modulus, givelog, sign);
73}
74
75SEXP BunchKaufman_determinant(SEXP obj, SEXP logarithm)
76{
77 DETERMINANT_START;
78
79 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym));
80 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
81
82 if (n > 0) {
83 SEXP pivot = GET_SLOT(obj, Matrix_permSym);
84 int *ppivot = INTEGER(pivot);
85
86 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
87 char ul = *CHAR(STRING_ELT(uplo, 0));
88
89 int j = 0, unpacked = (Matrix_int_fast64_t) n * n <= R_XLEN_T_MAX &&
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;
94 while (j < n) {
95 if (ppivot[j] > 0) {
96 modulus += log(hypot((*px).r, (*px).i));
97 px += (unpacked) ? n1a : ((ul == 'U') ? j + 2 : n - j);
98 j += 1;
99 } else {
100 a = *px;
101 if (ul == 'U') {
102 px += (unpacked) ? n1a : j + 2;
103 b = *px;
104 c = *(px - 1);
105 px += (unpacked) ? n1a : j + 3;
106 } else {
107 c = *(px + 1);
108 px += (unpacked) ? n1a : n - j;
109 b = *px;
110 px += (unpacked) ? n1a : n - j - 1;
111 }
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 -
115 2.0 * c.r * c.i));
116 j += 2;
117 }
118 }
119 } else {
120 double *px = REAL(x), a, b, c, logab, logcc;
121 while (j < n) {
122 if (ppivot[j] > 0) {
123 if (*px >= 0.0)
124 modulus += log(*px);
125 else {
126 modulus += log(-(*px));
127 sign = -sign;
128 }
129 px += (unpacked) ? n1a : ((ul == 'U') ? j + 2 : n - j);
130 j += 1;
131 } else {
132 a = *px;
133 if (ul == 'U') {
134 px += (unpacked) ? n1a : j + 2;
135 b = *px;
136 c = *(px - 1);
137 px += (unpacked) ? n1a : j + 3;
138 } else {
139 c = *(px + 1);
140 px += (unpacked) ? n1a : n - j;
141 b = *px;
142 px += (unpacked) ? n1a : n - j - 1;
143 }
144 logab = log(fabs(a)) + log(fabs(b));
145 logcc = 2.0 * log(fabs(c));
146 if ((a < 0.0) != (b < 0.0)) {
147 /* det = ab - cc = -(abs(ab) + cc) < 0 */
148 modulus += logspace_add(logab, logcc);
149 sign = -sign;
150 } else if (logab < logcc) {
151 /* det = ab - cc = -(cc - ab) < 0 */
152 modulus += logspace_sub(logcc, logab);
153 sign = -sign;
154 } else {
155 /* det = ab - cc > 0 */
156 modulus += logspace_sub(logab, logcc);
157 }
158 j += 2;
159 }
160 }
161 }
162 }
163
164 UNPROTECT(1); /* x */
165 return mkDet(modulus, givelog, sign);
166}
167
168SEXP Cholesky_determinant(SEXP obj, SEXP logarithm)
169{
170 DETERMINANT_START;
171
172 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym));
173 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
174
175 if (n > 0) {
176 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
177 char ul = *CHAR(STRING_ELT(uplo, 0));
178
179 int j, unpacked = (Matrix_int_fast64_t) n * n <= R_XLEN_T_MAX &&
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);
187 }
188 } else {
189 double *px = REAL(x);
190 for (j = 0; j < n; ++j) {
191 if (ISNAN(*px) || *px >= 0.0)
192 modulus += log(*px);
193 else {
194 modulus += log(-(*px));
195 sign = -sign;
196 }
197 px += (unpacked) ? n1a : ((ul == 'U') ? j + 2 : n - j);
198 }
199 }
200 modulus *= 2.0;
201 }
202
203 UNPROTECT(1); /* x */
204 return mkDet(modulus, givelog, sign);
205}
206
207SEXP sparseLU_determinant(SEXP obj, SEXP logarithm)
208{
209 DETERMINANT_START;
210
211 SEXP U = PROTECT(GET_SLOT(obj, Matrix_USym)),
212 x = PROTECT(GET_SLOT(U, Matrix_xSym));
213 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
214
215 if (n > 0) {
216 SEXP p = PROTECT(GET_SLOT(U, Matrix_pSym)),
217 i = PROTECT(GET_SLOT(U, Matrix_iSym));
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) {
222 kend = pp[j];
223 if (k < kend && pi[kend - 1] == j)
224 modulus += log(hypot(px[kend - 1].r, px[kend - 1].i));
225 else {
226 UNPROTECT(4); /* i, p, x, U */
227 return mkDet(R_NegInf, givelog, 1);
228 }
229 k = kend;
230 }
231 } else {
232 double *px = REAL(x);
233 for (j = 0; j < n; ++j) {
234 kend = pp[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]);
238 else {
239 modulus += log(-px[kend - 1]);
240 sign = -sign;
241 }
242 } else {
243 UNPROTECT(4); /* i, p, x, U */
244 return mkDet(R_NegInf, givelog, 1);
245 }
246 k = kend;
247 }
248
249 /* defined in ./perm.c : */
250 int signPerm(const int *, int, int);
251
252 p = GET_SLOT(obj, Matrix_pSym);
253 if (signPerm(INTEGER(p), LENGTH(p), 0) < 0)
254 sign = -sign;
255 p = GET_SLOT(obj, Matrix_qSym);
256 if (signPerm(INTEGER(p), LENGTH(p), 0) < 0)
257 sign = -sign;
258 }
259 UNPROTECT(2); /* i, p */
260 }
261
262 UNPROTECT(2); /* x, U */
263 return mkDet(modulus, givelog, sign);
264}
265
266SEXP sparseQR_determinant(SEXP obj, SEXP logarithm)
267{
268 DETERMINANT_START;
269
270 SEXP R = PROTECT(GET_SLOT(obj, Matrix_RSym)),
271 x = PROTECT(GET_SLOT(R, Matrix_xSym));
272 int sign = (TYPEOF(x) == CPLXSXP) ? NA_INTEGER : 1;
273
274 dim = GET_SLOT(R, Matrix_DimSym);
275 if (INTEGER(dim)[0] > n)
276 error(_("%s(<%s>) does not support structurally rank deficient case"),
277 "determinant", "sparseQR");
278
279 if (n > 0) {
280 SEXP p = PROTECT(GET_SLOT(R, Matrix_pSym)),
281 i = PROTECT(GET_SLOT(R, Matrix_iSym));
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) {
286 kend = pp[j];
287 if (k < kend && pi[kend - 1] == j)
288 modulus += log(hypot(px[kend - 1].r, px[kend - 1].i));
289 else {
290 UNPROTECT(4); /* i, p, x, U */
291 return mkDet(R_NegInf, givelog, 1);
292 }
293 k = kend;
294 }
295 } else {
296 double *px = REAL(x);
297 for (j = 0; j < n; ++j) {
298 kend = pp[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]);
302 else {
303 modulus += log(-px[kend - 1]);
304 sign = -sign;
305 }
306 } else {
307 UNPROTECT(4); /* i, p, x, R */
308 return mkDet(R_NegInf, givelog, 1);
309 }
310 k = kend;
311 }
312
313 /* defined in ./perm.c : */
314 int signPerm(const int *, int, int);
315
316 p = GET_SLOT(obj, Matrix_pSym);
317 if (signPerm(INTEGER(p), LENGTH(p), 0) < 0)
318 sign = -sign;
319 p = GET_SLOT(obj, Matrix_qSym);
320 if (signPerm(INTEGER(p), LENGTH(p), 0) < 0)
321 sign = -sign;
322 if (n % 2)
323 sign = -sign;
324 }
325 UNPROTECT(2); /* i, p */
326 }
327
328 UNPROTECT(2); /* x, R */
329 return mkDet(modulus, givelog, sign);
330}
331
332SEXP CHMfactor_determinant(SEXP obj, SEXP logarithm, SEXP sqrt)
333{
334 DETERMINANT_START;
335
336 cholmod_factor *L = M2CHF(obj, 1);
337 int sign = (L->xtype == CHOLMOD_COMPLEX) ? NA_INTEGER : 1;
338
339 if (n > 0) {
340 int j, sqrt_ = asLogical(sqrt);
341 if (L->is_super) {
342 int k, nc,
343 nsuper = (int) L->nsuper,
344 *psuper = (int *) L->super,
345 *ppi = (int *) L->pi,
346 *ppx = (int *) L->px;
347 R_xlen_t nr1a;
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;
353 px_ = px + ppx[k];
354 for (j = 0; j < nc; ++j) {
355 modulus += log(hypot((*px_).r, (*px_).i));
356 px_ += nr1a;
357 }
358 }
359 } else {
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;
364 px_ = px + ppx[k];
365 for (j = 0; j < nc; ++j) {
366 modulus += log(*px_);
367 px_ += nr1a;
368 }
369 }
370 }
371 modulus *= 2.0;
372 } else {
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));
378 if (L->is_ll)
379 modulus *= 2.0;
380 } else {
381 double *px = (double *) L->x;
382 if (L->is_ll) {
383 for (j = 0; j < n; ++j)
384 modulus += log(px[pp[j]]);
385 modulus *= 2.0;
386 } else {
387 for (j = 0; j < n; ++j) {
388 if (ISNAN(px[pp[j]]) || px[pp[j]] >= 0.0)
389 modulus += log(px[pp[j]]);
390 else {
391 if (sqrt_)
392 return mkDet(R_NaN, givelog, 1);
393 modulus += log(-px[pp[j]]);
394 sign = -sign;
395 }
396 }
397 }
398 }
399 }
400 if (sqrt_)
401 modulus *= 0.5;
402 }
403
404 return mkDet(modulus, givelog, sign);
405}
Matrix_int_fast64_t
long long Matrix_int_fast64_t
Definition Mdefines.h:27
_
#define _(String)
Definition Mdefines.h:44
GET_SLOT
#define GET_SLOT(x, what)
Definition Mdefines.h:85
Matrix_permSym
SEXP Matrix_permSym
Definition Msymbols.h:18
Matrix_DimSym
SEXP Matrix_DimSym
Definition Msymbols.h:3
Matrix_RSym
SEXP Matrix_RSym
Definition Msymbols.h:6
Matrix_xSym
SEXP Matrix_xSym
Definition Msymbols.h:22
Matrix_iSym
SEXP Matrix_iSym
Definition Msymbols.h:13
Matrix_uploSym
SEXP Matrix_uploSym
Definition Msymbols.h:21
Matrix_qSym
SEXP Matrix_qSym
Definition Msymbols.h:19
Matrix_USym
SEXP Matrix_USym
Definition Msymbols.h:8
Matrix_pSym
SEXP Matrix_pSym
Definition Msymbols.h:17
M2CHF
cholmod_factor * M2CHF(SEXP obj, int values)
Definition cholmod-etc.c:8
c
cholmod_common c
Definition cholmod-etc.c:5
cl
cholmod_common cl
Definition cholmod-etc.c:6
denseLU_determinant
SEXP denseLU_determinant(SEXP obj, SEXP logarithm)
Definition determinant.c:26
sparseLU_determinant
SEXP sparseLU_determinant(SEXP obj, SEXP logarithm)
Definition determinant.c:207
sparseQR_determinant
SEXP sparseQR_determinant(SEXP obj, SEXP logarithm)
Definition determinant.c:266
Cholesky_determinant
SEXP Cholesky_determinant(SEXP obj, SEXP logarithm)
Definition determinant.c:168
DETERMINANT_START
#define DETERMINANT_START
mkDet
static SEXP mkDet(double modulus, int logarithm, int sign)
Definition determinant.c:7
CHMfactor_determinant
SEXP CHMfactor_determinant(SEXP obj, SEXP logarithm, SEXP sqrt)
Definition determinant.c:332
BunchKaufman_determinant
SEXP BunchKaufman_determinant(SEXP obj, SEXP logarithm)
Definition determinant.c:75
signPerm
int signPerm(const int *p, int n, int off)
Definition perm.c:23
Generated by doxygen 1.9.8