Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
dtCMatrix.c
Go to the documentation of this file.
1  /* Sparse triangular numeric matrices */
2 #include "dtCMatrix.h"
3 #include "cs_utils.h"
4 
5 #define RETURN(_CH_) UNPROTECT(1); return (_CH_);
6 
7 /* This is used for *BOTH* triangular and symmetric Csparse: */
8 SEXP tCMatrix_validate(SEXP x)
9 {
10  SEXP val = xCMatrix_validate(x);/* checks x slot */
11  if(isString(val))
12  return(val);
13  else {
14  SEXP
15  islot = GET_SLOT(x, Matrix_iSym),
16  pslot = GET_SLOT(x, Matrix_pSym);
17  int uploT = (*uplo_P(x) == 'U'),
18  k, nnz = length(islot),
19  *xi = INTEGER(islot),
20  *xj = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
21 
22  expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xj);
23 
24  /* Maybe FIXME: ">" should be ">=" for diag = 'U' (uplo = 'U') */
25  if(uploT) {
26  for (k = 0; k < nnz; k++)
27  if(xi[k] > xj[k]) {
28  RETURN(mkString(_("uplo='U' must not have sparse entries below the diagonal")));
29  }
30  }
31  else {
32  for (k = 0; k < nnz; k++)
33  if(xi[k] < xj[k]) {
34  RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
35  }
36  }
37 
38  RETURN(ScalarLogical(1));
39  }
40 }
41 
42 /* This is used for *BOTH* triangular and symmetric Rsparse: */
43 SEXP tRMatrix_validate(SEXP x)
44 {
45  SEXP val = xRMatrix_validate(x);/* checks x slot */
46  if(isString(val))
47  return(val);
48  else {
49  SEXP
50  jslot = GET_SLOT(x, Matrix_jSym),
51  pslot = GET_SLOT(x, Matrix_pSym);
52  int uploT = (*uplo_P(x) == 'U'),
53  k, nnz = length(jslot),
54  *xj = INTEGER(jslot),
55  *xi = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
56 
57  expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xi);
58 
59  /* Maybe FIXME: ">" should be ">=" for diag = 'U' (uplo = 'U') */
60  if(uploT) {
61  for (k = 0; k < nnz; k++)
62  if(xi[k] > xj[k]) {
63  RETURN(mkString(_("uplo='U' must not have sparse entries below the diagonal")));
64  }
65  }
66  else {
67  for (k = 0; k < nnz; k++)
68  if(xi[k] < xj[k]) {
69  RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
70  }
71  }
72 
73  RETURN(ScalarLogical(1));
74  }
75 }
76 
77 SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
78 {
79  int cl = asLogical(classed);
80  SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
81  CSP A = AS_CSP(a);
82  int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
83  *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
84  getAttrib(b, R_DimSymbol));
85  int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
86  double *bx;
87  R_CheckStack();
88 
89  if (adims[0] != n || n != adims[1])
90  error(_("Dimensions of system to be solved are inconsistent"));
91  Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
92  // dimnames:
93  SEXP dn = PROTECT(allocVector(VECSXP, 2)), dn2;
94  SET_VECTOR_ELT(dn, 0, duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1)));
95  if(!cl) {
96  dn2 = getAttrib(b, R_DimNamesSymbol);
97  if(dn2 != R_NilValue) // either NULL or list(<dn1>, <dn2>)
98  dn2 = VECTOR_ELT(dn2, 1);
99  }
100  SET_VECTOR_ELT(dn, 1, duplicate(cl // b can be "Matrix" or not:
101  ? VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)
102  : dn2));
103  SET_SLOT(ans, Matrix_DimNamesSym, dn);
104  UNPROTECT(1);
105  if(n >= 1 && nrhs >=1) {
106  bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
107  REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
108  for (j = 0; j < nrhs; j++)
109  lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
110  }
111  RETURN(ans);
112 }
113 
114 SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)
115 {
116  SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
117  CSP A = AS_CSP(a), B = AS_CSP(b);
118  R_CheckStack();
119  if (A->m != A->n || B->n < 1 || A->n < 1 || A->n != B->m)
120  error(_("Dimensions of system to be solved are inconsistent"));
121  // *before* Calloc()ing below [memory leak]! -- FIXME: 0-extent should work
122 
123  int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
124  xnz = 10 * B->p[B->n]; /* initial estimate of nnz in x */
125  int k, lo = uplo_P(a)[0] == 'L', pos = 0;
126  int *ti = Calloc(xnz, int), *xi = Calloc(2*A->n, int); /* for cs_reach */
127  double *tx = Calloc(xnz, double), *wrk = Calloc( A->n, double);
128 
129  slot_dup(ans, b, Matrix_DimSym);
130 
131  xp[0] = 0;
132  for (k = 0; k < B->n; k++) {
133  int top = cs_spsolve (A, B, k, xi, wrk, (int *)NULL, lo);
134  int nz = A->n - top;
135 
136  xp[k + 1] = nz + xp[k];
137  if (xp[k + 1] > xnz) {
138  while (xp[k + 1] > xnz) xnz *= 2;
139  ti = Realloc(ti, xnz, int);
140  tx = Realloc(tx, xnz, double);
141  }
142  if (lo) /* increasing row order */
143  for(int p = top; p < A->n; p++, pos++) {
144  ti[pos] = xi[p];
145  tx[pos] = wrk[xi[p]];
146  }
147  else /* decreasing order, reverse copy */
148  for(int p = A->n - 1; p >= top; p--, pos++) {
149  ti[pos] = xi[p];
150  tx[pos] = wrk[xi[p]];
151  }
152  }
153  xnz = xp[B->n];
154  Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, xnz)), ti, xnz);
155  Memcpy( REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, xnz)), tx, xnz);
156 
157  Free(ti); Free(tx);
158  Free(wrk); Free(xi);
159 
160  // dimnames:
161  SEXP dn = PROTECT(allocVector(VECSXP, 2));
162  SET_VECTOR_ELT(dn, 0, duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1)));
163  SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));
164  SET_SLOT(ans, Matrix_DimNamesSym, dn);
165  UNPROTECT(1);
166 
167  RETURN(ans);
168 }
169 #undef RETURN
170 
SEXP Matrix_DimSym
Definition: Syms.h:2
SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)
Definition: dtCMatrix.c:114
SEXP Matrix_xSym
Definition: Syms.h:2
#define slot_dup(dest, src, sym)
Definition: Mutils.h:149
csi n
Definition: cs.h:37
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
csi m
Definition: cs.h:36
cholmod_common cl
Definition: chm_common.c:16
#define AS_CSP(x)
Definition: cs_utils.h:12
SEXP Matrix_jSym
Definition: Syms.h:2
#define RETURN(_CH_)
Definition: dtCMatrix.c:5
SEXP xCMatrix_validate(SEXP x)
Definition: dgCMatrix.c:11
SEXP xRMatrix_validate(SEXP x)
Definition: dgCMatrix.c:23
static R_INLINE int * expand_cmprPt(int ncol, const int mp[], int mj[])
Expand compressed pointers in the array mp into a full set of indices in the array mj...
Definition: Mutils.h:259
#define uplo_P(_x_)
Definition: Mutils.h:174
#define _(String)
Definition: Mutils.h:32
csi cs_usolve(const cs *U, double *x)
Definition: cs.c:1890
SEXP Matrix_iSym
Definition: Syms.h:2
Definition: cs.h:33
SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
Definition: dtCMatrix.c:77
csi cs_spsolve(cs *G, const cs *B, csi k, csi *xi, double *x, const csi *pinv, csi lo)
Definition: cs.c:1656
SEXP Matrix_pSym
Definition: Syms.h:2
SEXP tCMatrix_validate(SEXP x)
Definition: dtCMatrix.c:8
csi cs_lsolve(const cs *L, double *x)
Definition: cs.c:985
static R_INLINE SEXP ALLOC_SLOT(SEXP obj, SEXP nm, SEXPTYPE type, int length)
Allocate an SEXP of given type and length, assign it as slot nm in the object, and return the SEXP...
Definition: Mutils.h:240
SEXP tRMatrix_validate(SEXP x)
Definition: dtCMatrix.c:43