Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
CHMfactor.c
Go to the documentation of this file.
1  /* CHOLMOD factors */
2 #include "CHMfactor.h"
3 
4 SEXP CHMfactor_to_sparse(SEXP x)
5 {
6  CHM_FR L = AS_CHM_FR(x), Lcp;
7  CHM_SP Lm;
8  R_CheckStack();
9 
10  /* cholmod_factor_to_sparse changes its first argument. Make a copy */
11  Lcp = cholmod_copy_factor(L, &c);
12  if (!(Lcp->is_ll))
13  if (!cholmod_change_factor(Lcp->xtype, 1, 0, 1, 1, Lcp, &c))
14  error(_("cholmod_change_factor failed with status %d"), c.status);
15  Lm = cholmod_factor_to_sparse(Lcp, &c); cholmod_free_factor(&Lcp, &c);
16  return chm_sparse_to_SEXP(Lm, 1/*do_free*/, -1/*uploT*/, 0/*Rkind*/,
17  "N"/*non_unit*/, R_NilValue/*dimNames*/);
18 }
19 
20 SEXP CHMfactor_solve(SEXP a, SEXP b, SEXP system)
21 {
22  CHM_FR L = AS_CHM_FR(a);
23  SEXP bb = PROTECT(dup_mMatrix_as_dgeMatrix(b));
24  CHM_DN B = AS_CHM_DN(bb), X;
25  int sys = asInteger(system);
26  R_CheckStack();
27 
28  if (!(sys--)) /* align with CHOLMOD defs: R's {1:9} --> {0:8},
29  see ./CHOLMOD/Cholesky/cholmod_solve.c */
30  error(_("system argument is not valid"));
31 
32  X = cholmod_solve(sys, L, B, &c);
33  UNPROTECT(1);
34  return chm_dense_to_SEXP(X, 1/*do_free*/, 0/*Rkind*/,
35  GET_SLOT(bb, Matrix_DimNamesSym), FALSE);
36 }
37 
38 SEXP CHMfactor_updown(SEXP upd, SEXP C_, SEXP L_)
39 {
40  CHM_FR L = AS_CHM_FR(L_), Lcp;
41  CHM_SP C = AS_CHM_SP__(C_);
42  int update = asInteger(upd);
43  R_CheckStack();
44 
45  Lcp = cholmod_copy_factor(L, &c);
46  int r = cholmod_updown(update, C, Lcp, &c);
47  if(!r) error(_("cholmod_updown() returned %d"), r);
48  return chm_factor_to_SEXP(Lcp, 1);
49 }
50 
51 SEXP CHMfactor_spsolve(SEXP a, SEXP b, SEXP system)
52 {
53  CHM_FR L = AS_CHM_FR(a);
54  CHM_SP B = AS_CHM_SP__(b);
55  int sys = asInteger(system);
56  R_CheckStack();
57 
58  if (!(sys--)) /* align with CHOLMOD defs: R's {1:9} --> {0:8},
59  see ./CHOLMOD/Cholesky/cholmod_solve.c */
60  error(_("system argument is not valid"));
61 
62  // dimnames:
63  SEXP dn = PROTECT(allocVector(VECSXP, 2));
64  // none from a: our CHMfactor objects have no dimnames
65  SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));
66  UNPROTECT(1);
67 
68  return chm_sparse_to_SEXP(cholmod_spsolve(sys, L, B, &c),
69  1/*do_free*/, 0/*uploT*/, 0/*Rkind*/, "", dn);
70 }
71 
81 {
82  int i, j, p;
83  double ans = 0;
84 
85  if (f->is_super) {
86  int *lpi = (int*)(f->pi), *lsup = (int*)(f->super);
87  for (i = 0; i < f->nsuper; i++) { /* supernodal block i */
88  int nrp1 = 1 + lpi[i + 1] - lpi[i],
89  nc = lsup[i + 1] - lsup[i];
90  double *x = (double*)(f->x) + ((int*)(f->px))[i];
91 
92  for (j = 0; j < nc; j++) {
93  ans += 2 * log(fabs(x[j * nrp1]));
94  }
95  }
96  } else {
97  int *li = (int*)(f->i), *lp = (int*)(f->p);
98  double *lx = (double *)(f->x);
99 
100  for (j = 0; j < f->n; j++) {
101  for (p = lp[j]; li[p] != j && p < lp[j + 1]; p++) {};
102  if (li[p] != j) {
103  error(_("diagonal element %d of Cholesky factor is missing"), j);
104  break; /* -Wall */
105  }
106  ans += log(lx[p] * ((f->is_ll) ? lx[p] : 1.));
107  }
108  }
109  return ans;
110 }
111 
112 SEXP CHMfactor_ldetL2(SEXP x)
113 {
114  CHM_FR L = AS_CHM_FR(x); R_CheckStack();
115 
116  return ScalarReal(chm_factor_ldetL2(L));
117 }
118 
134 {
135  int ll = f->is_ll;
136  double mm[2] = {0, 0};
137  mm[0] = mult;
138  // NB: Result depends if A is "dsC" or "dgC"; the latter case assumes we mean AA' !!!
139  if (!cholmod_factorize_p(A, mm, (int*)NULL, 0 /*fsize*/, f, &c))
140  /* -> ./CHOLMOD/Cholesky/cholmod_factorize.c */
141  error(_("cholmod_factorize_p failed: status %d, minor %d of ncol %d"),
142  c.status, f->minor, f->n);
143  if (f->is_ll != ll)
144  if(!cholmod_change_factor(f->xtype, ll, f->is_super, 1 /*to_packed*/,
145  1 /*to_monotonic*/, f, &c))
146  error(_("cholmod_change_factor failed"));
147  return f;
148 }
149 
150 SEXP CHMfactor_update(SEXP object, SEXP parent, SEXP mult)
151 {
152  CHM_FR L = AS_CHM_FR(object), Lcp;
153  CHM_SP A = AS_CHM_SP__(parent);
154  R_CheckStack();
155 
156  Lcp = cholmod_copy_factor(L, &c);
157  return chm_factor_to_SEXP(chm_factor_update(Lcp, A, asReal(mult)), 1);
158 }
159 
160 SEXP destructive_CHM_update(SEXP object, SEXP parent, SEXP mult)
161 {
162  CHM_FR L = AS_CHM_FR(object);
163  CHM_SP A = AS_CHM_SP__(parent);
164  R_CheckStack();
165 
166  return chm_factor_to_SEXP(chm_factor_update(L, A, asReal(mult)), 0);
167 }
168 
169 SEXP CHMfactor_ldetL2up(SEXP x, SEXP parent, SEXP mult)
170 {
171  SEXP ans = PROTECT(duplicate(mult));
172  int i, nmult = LENGTH(mult);
173  double *aa = REAL(ans), *mm = REAL(mult);
174  CHM_FR L = AS_CHM_FR(x), Lcp;
175  CHM_SP A = AS_CHM_SP__(parent);
176  R_CheckStack();
177 
178  Lcp = cholmod_copy_factor(L, &c);
179  for (i = 0; i < nmult; i++)
180  aa[i] = chm_factor_ldetL2(chm_factor_update(Lcp, A, mm[i]));
181  cholmod_free_factor(&Lcp, &c);
182  UNPROTECT(1);
183  return ans;
184 }
SEXP CHMfactor_ldetL2up(SEXP x, SEXP parent, SEXP mult)
Definition: CHMfactor.c:169
SEXP chm_dense_to_SEXP(CHM_DN a, int dofree, int Rkind, SEXP dn, Rboolean transp)
Copy the contents of a to an appropriate denseMatrix object and, optionally, free a or free both a an...
Definition: chm_common.c:835
#define AS_CHM_SP__(x)
Definition: chm_common.h:49
SEXP CHMfactor_updown(SEXP upd, SEXP C_, SEXP L_)
Definition: CHMfactor.c:38
cholmod_sparse * CHM_SP
Definition: chm_common.h:25
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
SEXP CHMfactor_to_sparse(SEXP x)
Definition: CHMfactor.c:4
SEXP chm_sparse_to_SEXP(CHM_SP a, int dofree, int uploT, int Rkind, const char *diag, SEXP dn)
Copy the contents of a to an appropriate CsparseMatrix object and, optionally, free a or free both a ...
Definition: chm_common.c:335
SEXP dup_mMatrix_as_dgeMatrix(SEXP A)
Definition: Mutils.c:852
SEXP CHMfactor_spsolve(SEXP a, SEXP b, SEXP system)
Definition: CHMfactor.c:51
SEXP destructive_CHM_update(SEXP object, SEXP parent, SEXP mult)
Definition: CHMfactor.c:160
cholmod_factor * CHM_FR
Definition: chm_common.h:23
double chm_factor_ldetL2(CHM_FR f)
Evaluate the logarithm of the square of the determinant of L.
Definition: CHMfactor.c:80
#define _(String)
Definition: Mutils.h:32
cholmod_dense * CHM_DN
Definition: chm_common.h:21
SEXP CHMfactor_update(SEXP object, SEXP parent, SEXP mult)
Definition: CHMfactor.c:150
SEXP chm_factor_to_SEXP(CHM_FR f, int dofree)
Copy the contents of f to an appropriate CHMfactor object and, optionally, free f or free both f and ...
Definition: chm_common.c:1106
#define AS_CHM_DN(x)
Definition: chm_common.h:43
#define AS_CHM_FR(x)
Definition: chm_common.h:45
cholmod_common c
Definition: chm_common.c:15
SEXP CHMfactor_ldetL2(SEXP x)
Definition: CHMfactor.c:112
SEXP CHMfactor_solve(SEXP a, SEXP b, SEXP system)
Definition: CHMfactor.c:20
CHM_FR chm_factor_update(CHM_FR f, CHM_SP A, double mult)
Update the numerical values in the factor f as A + mult * I, if A is symmetric, otherwise AA' + mult ...
Definition: CHMfactor.c:133