Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
dsCMatrix.c
Go to the documentation of this file.
1 #include "dsCMatrix.h"
2 
3 static int chk_nm(const char *nm, int perm, int LDL, int super)
4 {
5  if (strlen(nm) != 11) return 0;
6  if (strcmp(nm + 3, "Cholesky")) return 0;
7  if (super > 0 && nm[0] != 'S') return 0;
8  if (super == 0 && nm[0] != 's') return 0;
9  if (perm > 0 && nm[1] != 'P') return 0;
10  if (perm == 0 && nm[1] != 'p') return 0;
11  if (LDL > 0 && nm[2] != 'D') return 0;
12  if (LDL == 0 && nm[2] != 'd') return 0;
13  return 1;
14 }
15 
16 SEXP R_chkName_Cholesky(SEXP nm, SEXP perm, SEXP LDL, SEXP super)
17 {
18  return ScalarLogical(chk_nm(CHAR(asChar(nm)), asLogical(perm),
19  asLogical(LDL), asLogical(super)));
20 }
21 
22 // must be called with 'nm' a string of length 11
23 static void chm_factor_name(char* nm, int perm, int LDL, int super) {
24  if (strlen(nm) != 11) {
25  error(_("chm_factor_name(): did not get string of length 11"));
26  return;
27  }
28  nm[0] = (super > 0) ? 'S' : 's';
29  nm[1] = (perm == 0) ? 'p' : 'P';
30  nm[2] = (LDL == 0) ? 'd' : 'D';
31  return;
32 }
33 
34 // must be called with 'nm' a string of length 11
35 SEXP R_chm_factor_name(SEXP perm, SEXP LDL, SEXP super)
36 {
37  char nm[12] = "...Cholesky";// 11 + final \0
38  chm_factor_name(nm, asLogical(perm), asLogical(LDL), asLogical(super));
39  return mkString(nm);
40 }
41 
42 
43 
61 static CHM_FR
62 internal_chm_factor(SEXP Ap, int perm, int LDL, int super, double Imult)
63 {
64  SEXP facs = GET_SLOT(Ap, Matrix_factorSym);
65  SEXP nms = getAttrib(facs, R_NamesSymbol);
66  CHM_FR L;
67  CHM_SP A = AS_CHM_SP__(Ap);
68  R_CheckStack();
69 
70  CHM_store_common(); /* save settings from c */
71  if (LENGTH(facs)) {
72  for (int i = 0; i < LENGTH(nms); i++) { /* look for a match in cache */
73  if (chk_nm(CHAR(STRING_ELT(nms, i)), perm, LDL, super)) {
74  L = AS_CHM_FR(VECTOR_ELT(facs, i));
75  R_CheckStack();
76  /* copy the factor so later it can safely be cholmod_free'd */
77  L = cholmod_copy_factor(L, &c);
78  if (Imult) cholmod_factorize_p(A, &Imult, (int*)NULL, 0, L, &c);
79  return L;
80  }
81  }
82  }
83  /* Else: No cached factor - create one */
84 
85  c.final_ll = (LDL == 0) ? 1 : 0;
86  c.supernodal = (super > 0) ? CHOLMOD_SUPERNODAL :
87  ((super < 0) ? CHOLMOD_AUTO :
88  /* super == 0 */ CHOLMOD_SIMPLICIAL);
89 
90  if (perm) { /* obtain fill-reducing permutation */
91  L = cholmod_analyze(A, &c);
92  } else { /* require identity permutation */
93  c.nmethods = 1; c.method[0].ordering = CHOLMOD_NATURAL; c.postorder = FALSE;
94  // *_restore_*() below or in R_cholmod_error() will restore c.<foo>
95  L = cholmod_analyze(A, &c);
96  }
97  if (!cholmod_factorize_p(A, &Imult, (int*)NULL, 0 /*fsize*/, L, &c))
98  // have never seen this, rather R_cholmod_error(status, ..) is called :
99  error(_("Cholesky factorization failed; unusually, please report to Matrix-authors"));
100 
101  if (!Imult) { /* cache the factor */
102  if(!chm_factor_ok(L)) {
103  cholmod_free_factor(&L, &c);// <- do not leak!
105  error(_("internal_chm_factor: Cholesky factorization failed"));
106  }
107 
108  /* now that we allow (super, LDL) to be "< 0", be careful :*/
109  if(super < 0) super = L->is_super ? 1 : 0;
110  if(LDL < 0) LDL = L->is_ll ? 0 : 1;
111 
112  char fnm[12] = "...Cholesky";// 11 + final \0
113  chm_factor_name(fnm, perm, LDL, super);
114 
115  set_factors(Ap, chm_factor_to_SEXP(L, 0), fnm);
116  }
118  return L;
119 }
120 
121 SEXP dsCMatrix_chol(SEXP x, SEXP pivot)
122 {
123  int pivP = asLogical(pivot);
124  CHM_FR L = internal_chm_factor(x, pivP, /*LDL = */ 0, /* super = */ 0,
125  /* Imult = */ 0.);
126  CHM_SP R, Rt;
127  SEXP ans;
128 
129  Rt = cholmod_factor_to_sparse(L, &c);
130  R = cholmod_transpose(Rt, /*values*/ 1, &c);
131  cholmod_free_sparse(&Rt, &c);
132  ans = PROTECT(chm_sparse_to_SEXP(R, 1/*do_free*/, 1/*uploT*/, 0/*Rkind*/,
133  "N"/*diag*/, GET_SLOT(x, Matrix_DimNamesSym)));
134 
135  if (pivP) {
136  SEXP piv = PROTECT(allocVector(INTSXP, L->n));
137  int *dest = INTEGER(piv), *src = (int*)L->Perm;
138 
139  for (int i = 0; i < L->n; i++) dest[i] = src[i] + 1;
140  setAttrib(ans, install("pivot"), piv);
141  setAttrib(ans, install("rank"), ScalarInteger((size_t) L->minor));
142  UNPROTECT(1);
143  }
144  cholmod_free_factor(&L, &c);
145  UNPROTECT(1);
146  return ans;
147 }
148 
149 SEXP dsCMatrix_Cholesky(SEXP Ap, SEXP perm, SEXP LDL, SEXP super, SEXP Imult)
150 {
151  int iSuper = asLogical(super),
152  iPerm = asLogical(perm),
153  iLDL = asLogical(LDL);
154 
155  /* When parameter is set to NA in R, let CHOLMOD choose */
156  if(iSuper == NA_LOGICAL) iSuper = -1;
157  /* if(iPerm == NA_LOGICAL) iPerm = -1; */
158  if(iLDL == NA_LOGICAL) iLDL = -1;
159  SEXP r = chm_factor_to_SEXP(internal_chm_factor(Ap, iPerm, iLDL, iSuper,
160  asReal(Imult)),
161  1 /* dofree */);
162  return r;
163 }
164 
178 SEXP dsCMatrix_LDL_D(SEXP Ap, SEXP permP, SEXP resultKind)
179 {
180  CHM_FR L;
181  SEXP ans;
182  L = internal_chm_factor(Ap, asLogical(permP),
183  /*LDL*/ 1, /*super*/0, /*Imult*/0.);
184  // ./Csparse.c :
185  ans = PROTECT(diag_tC_ptr(L->n,
186  L->p,
187  L->x,
188  /* is_U = */ FALSE,
189  L->Perm,
190  resultKind));
191  cholmod_free_factor(&L, &c);
192  UNPROTECT(1);
193  return(ans);
194 }
195 
196 // using cholmod_spsolve() --> sparse result
197 SEXP dsCMatrix_Csparse_solve(SEXP a, SEXP b, SEXP LDL)
198 {
199  int iLDL = asLogical(LDL);
200  // When parameter is set to NA in R, let CHOLMOD choose
201  if(iLDL == NA_LOGICAL) iLDL = -1;
202 
203  CHM_FR L = internal_chm_factor(a, /*perm*/-1, iLDL, /*super*/-1, /*Imult*/0.);
204  if(!chm_factor_ok(L)) {
205  cholmod_free_factor(&L, &c);
206  return R_NilValue;// == "CHOLMOD factorization failed"
207  }
208  CHM_SP cb = AS_CHM_SP(b), cx;
209  R_CheckStack();
210 
211  cx = cholmod_spsolve(CHOLMOD_A, L, cb, &c);
212  cholmod_free_factor(&L, &c);
213  return chm_sparse_to_SEXP(cx, /*do_free*/ 1, /*uploT*/ 0,
214  /*Rkind*/ 0, /*diag*/ "N",
215  /*dimnames = */ R_NilValue);
216 }
217 
218 // using cholmod_solve() --> dense result
219 SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b, SEXP LDL)
220 {
221  int iLDL = asLogical(LDL);
222  // When parameter is set to NA in R, let CHOLMOD choose
223  if(iLDL == NA_LOGICAL) iLDL = -1;
224 
225  CHM_FR L = internal_chm_factor(a, /*perm*/-1, iLDL, /*super*/-1, /*Imult*/0.);
226  if(!chm_factor_ok(L)) {
227  cholmod_free_factor(&L, &c);
228  return R_NilValue;// == "CHOLMOD factorization failed"
229  }
230 
231  CHM_DN cx, cb = AS_CHM_DN(PROTECT(mMatrix_as_dgeMatrix(b)));
232  R_CheckStack();
233  cx = cholmod_solve(CHOLMOD_A, L, cb, &c);
234  cholmod_free_factor(&L, &c);
235  UNPROTECT(1);
236  return chm_dense_to_SEXP(cx, 1, 0, /*dimnames = */ R_NilValue, /* transp: */ FALSE);
237 }
238 
239 /* Needed for printing dsCMatrix objects */
240 /* FIXME: Create a more general version of this operation: also for lsC, (dsR?),..
241 * e.g. make compressed_to_dgTMatrix() in ./dgCMatrix.c work for dsC */
243 {
244  CHM_SP A = AS_CHM_SP__(x);
245  CHM_SP Afull = cholmod_copy(A, /*stype*/ 0, /*mode*/ 1, &c);
246  CHM_TR At = cholmod_sparse_to_triplet(Afull, &c);
247  R_CheckStack();
248 
249  if (!A->stype)
250  error(_("Non-symmetric matrix passed to dsCMatrix_to_dgTMatrix"));
251  cholmod_free_sparse(&Afull, &c);
252  return chm_triplet_to_SEXP(At, 1, /*uploT*/ 0, /*Rkind*/ 0, "",
253  GET_SLOT(x, Matrix_DimNamesSym));
254 }
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
static int chk_nm(const char *nm, int perm, int LDL, int super)
Definition: dsCMatrix.c:3
static CHM_FR internal_chm_factor(SEXP Ap, int perm, int LDL, int super, double Imult)
Return a CHOLMOD copy of the cached Cholesky decomposition with the required perm, LDL and super attributes.
Definition: dsCMatrix.c:62
void CHM_store_common()
Definition: chm_common.c:25
SEXP dsCMatrix_Csparse_solve(SEXP a, SEXP b, SEXP LDL)
Definition: dsCMatrix.c:197
cholmod_sparse * CHM_SP
Definition: chm_common.h:25
SEXP diag_tC_ptr(int n, int *x_p, double *x_x, Rboolean is_U, int *perm, SEXP resultKind)
SEXP dsCMatrix_Cholesky(SEXP Ap, SEXP perm, SEXP LDL, SEXP super, SEXP Imult)
Definition: dsCMatrix.c:149
SEXP Matrix_factorSym
Definition: Syms.h:2
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
cholmod_triplet * CHM_TR
Definition: chm_common.h:27
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
static R_INLINE SEXP mMatrix_as_dgeMatrix(SEXP A)
Definition: Mutils.h:334
SEXP dsCMatrix_chol(SEXP x, SEXP pivot)
Definition: dsCMatrix.c:121
SEXP dsCMatrix_LDL_D(SEXP Ap, SEXP permP, SEXP resultKind)
Fast version of getting at the diagonal matrix D of the (generalized) simplicial Cholesky LDL' decomp...
Definition: dsCMatrix.c:178
SEXP R_chm_factor_name(SEXP perm, SEXP LDL, SEXP super)
Definition: dsCMatrix.c:35
SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b, SEXP LDL)
Definition: dsCMatrix.c:219
SEXP set_factors(SEXP obj, SEXP val, char *nm)
Caches 'val' in the 'factors' slot of obj, i.e.
Definition: Mutils.c:130
static R_INLINE Rboolean chm_factor_ok(CHM_FR f)
Definition: chm_common.h:57
cholmod_factor * CHM_FR
Definition: chm_common.h:23
#define _(String)
Definition: Mutils.h:32
cholmod_dense * CHM_DN
Definition: chm_common.h:21
SEXP dsCMatrix_to_dgTMatrix(SEXP x)
Definition: dsCMatrix.c:242
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
SEXP R_chkName_Cholesky(SEXP nm, SEXP perm, SEXP LDL, SEXP super)
Definition: dsCMatrix.c:16
#define AS_CHM_DN(x)
Definition: chm_common.h:43
static void chm_factor_name(char *nm, int perm, int LDL, int super)
Definition: dsCMatrix.c:23
#define AS_CHM_FR(x)
Definition: chm_common.h:45
void CHM_restore_common()
Definition: chm_common.c:50
cholmod_common c
Definition: chm_common.c:15
SEXP chm_triplet_to_SEXP(CHM_TR a, int dofree, int uploT, int Rkind, const char *diag, SEXP dn)
Copy the contents of a to an appropriate TsparseMatrix object and, optionally, free a or free both a ...
Definition: chm_common.c:573
#define AS_CHM_SP(x)
Definition: chm_common.h:46