Annotation of /pkg/src/Csparse.c

Parent Directory | Revision Log

---

Revision 2137 - (view) (download) (as text)

1 :	bates	1218	/* Sparse matrices in compressed column-oriented form */
2 :	bates	922	#include "Csparse.h"
3 :	maechler	2120	#include "Tsparse.h"
4 :	bates	922	#include "chm_common.h"
5 :
6 :			SEXP Csparse_validate(SEXP x)
7 :			{
8 :	maechler	1575	/* NB: we do *NOT* check a potential 'x' slot here, at all */
9 :	bates	922	SEXP pslot = GET_SLOT(x, Matrix_pSym),
10 :			islot = GET_SLOT(x, Matrix_iSym);
11 :	maechler	1893	Rboolean sorted, strictly;
12 :			int j, k,
13 :	bates	922	*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
14 :	maechler	1660	nrow = dims[0],
15 :			ncol = dims[1],
16 :	maechler	1654	*xp = INTEGER(pslot),
17 :	bates	922	*xi = INTEGER(islot);
18 :
19 :	maechler	1654	if (length(pslot) != dims[1] + 1)
20 :			return mkString(_("slot p must have length = ncol(.) + 1"));
21 :	bates	922	if (xp[0] != 0)
22 :			return mkString(_("first element of slot p must be zero"));
23 :	maechler	1893	if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/
24 :	bates	1555	return
25 :			mkString(_("last element of slot p must match length of slots i and x"));
26 :			for (j = 0; j < length(islot); j++) {
27 :			if (xi[j] < 0 || xi[j] >= nrow)
28 :			return mkString(_("all row indices must be between 0 and nrow-1"));
29 :			}
30 :	maechler	1893	sorted = TRUE; strictly = TRUE;
31 :	bates	922	for (j = 0; j < ncol; j++) {
32 :			if (xp[j] > xp[j+1])
33 :			return mkString(_("slot p must be non-decreasing"));
34 :	maechler	1893	if(sorted)
35 :			for (k = xp[j] + 1; k < xp[j + 1]; k++) {
36 :			if (xi[k] < xi[k - 1])
37 :			sorted = FALSE;
38 :			else if (xi[k] == xi[k - 1])
39 :			strictly = FALSE;
40 :			}
41 :	bates	922	}
42 :	maechler	1654	if (!sorted) {
43 :	maechler	1960	CHM_SP chx = AS_CHM_SP(x);
44 :			R_CheckStack();
45 :
46 :	maechler	1654	cholmod_sort(chx, &c);
47 :	maechler	1893	/* Now re-check that row indices are *strictly* increasing
48 :			* (and not just increasing) within each column : */
49 :			for (j = 0; j < ncol; j++) {
50 :			for (k = xp[j] + 1; k < xp[j + 1]; k++)
51 :			if (xi[k] == xi[k - 1])
52 :			return mkString(_("slot i is not *strictly* increasing inside a column (even after cholmod_sort)"));
53 :			}
54 :
55 :			} else if(!strictly) { /* sorted, but not strictly */
56 :			return mkString(_("slot i is not *strictly* increasing inside a column"));
57 :	maechler	1654	}
58 :	bates	922	return ScalarLogical(1);
59 :			}
60 :
61 :	maechler	1968	SEXP Rsparse_validate(SEXP x)
62 :			{
63 :			/* NB: we do *NOT* check a potential 'x' slot here, at all */
64 :			SEXP pslot = GET_SLOT(x, Matrix_pSym),
65 :			jslot = GET_SLOT(x, Matrix_jSym);
66 :			Rboolean sorted, strictly;
67 :			int i, k,
68 :			*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
69 :			nrow = dims[0],
70 :			ncol = dims[1],
71 :			*xp = INTEGER(pslot),
72 :			*xj = INTEGER(jslot);
73 :
74 :			if (length(pslot) != dims[0] + 1)
75 :			return mkString(_("slot p must have length = nrow(.) + 1"));
76 :			if (xp[0] != 0)
77 :			return mkString(_("first element of slot p must be zero"));
78 :			if (length(jslot) < xp[nrow]) /* allow larger slots from over-allocation!*/
79 :			return
80 :			mkString(_("last element of slot p must match length of slots j and x"));
81 :			for (i = 0; i < length(jslot); i++) {
82 :			if (xj[i] < 0 || xj[i] >= ncol)
83 :			return mkString(_("all column indices must be between 0 and ncol-1"));
84 :			}
85 :			sorted = TRUE; strictly = TRUE;
86 :			for (i = 0; i < nrow; i++) {
87 :			if (xp[i] > xp[i+1])
88 :			return mkString(_("slot p must be non-decreasing"));
89 :			if(sorted)
90 :			for (k = xp[i] + 1; k < xp[i + 1]; k++) {
91 :			if (xj[k] < xj[k - 1])
92 :			sorted = FALSE;
93 :			else if (xj[k] == xj[k - 1])
94 :			strictly = FALSE;
95 :			}
96 :			}
97 :			if (!sorted)
98 :			/* cannot easily use cholmod_sort(.) ... -> "error out" :*/
99 :			return mkString(_("slot j is not increasing inside a column"));
100 :			else if(!strictly) /* sorted, but not strictly */
101 :			return mkString(_("slot j is not *strictly* increasing inside a column"));
102 :
103 :			return ScalarLogical(1);
104 :			}
105 :
106 :
107 :	maechler	1751	/* Called from ../R/Csparse.R : */
108 :			/* Can only return [dln]geMatrix (no symm/triang);
109 :			* FIXME: replace by non-CHOLMOD code ! */
110 :	bates	1059	SEXP Csparse_to_dense(SEXP x)
111 :			{
112 :	maechler	1960	CHM_SP chxs = AS_CHM_SP(x);
113 :	maechler	1751	/* This loses the symmetry property, since cholmod_dense has none,
114 :			* BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices
115 :			* to numeric (CHOLMOD_REAL) ones : */
116 :	maechler	1960	CHM_DN chxd = cholmod_sparse_to_dense(chxs, &c);
117 :	maechler	1751	int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x);
118 :	maechler	1960	R_CheckStack();
119 :	bates	1059
120 :	maechler	1736	return chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym));
121 :	bates	1059	}
122 :
123 :	maechler	1548	SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri)
124 :	bates	1371	{
125 :	maechler	1960	CHM_SP chxs = AS_CHM_SP(x);
126 :			CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c);
127 :	bates	1867	int tr = asLogical(tri);
128 :	maechler	1960	R_CheckStack();
129 :	bates	1371
130 :	maechler	1960	return chm_sparse_to_SEXP(chxcp, 1/*do_free*/,
131 :	bates	1867	tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
132 :			0, tr ? diag_P(x) : "",
133 :	maechler	1548	GET_SLOT(x, Matrix_DimNamesSym));
134 :	bates	1371	}
135 :
136 :	bates	1366	SEXP Csparse_to_matrix(SEXP x)
137 :	bates	922	{
138 :	maechler	1960	return chm_dense_to_matrix(cholmod_sparse_to_dense(AS_CHM_SP(x), &c),
139 :			1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym));
140 :	bates	1366	}
141 :
142 :			SEXP Csparse_to_Tsparse(SEXP x, SEXP tri)
143 :			{
144 :	maechler	1960	CHM_SP chxs = AS_CHM_SP(x);
145 :			CHM_TR chxt = cholmod_sparse_to_triplet(chxs, &c);
146 :	bates	1867	int tr = asLogical(tri);
147 :	maechler	1736	int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
148 :	maechler	1960	R_CheckStack();
149 :	bates	922
150 :	bates	1867	return chm_triplet_to_SEXP(chxt, 1,
151 :			tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
152 :			Rkind, tr ? diag_P(x) : "",
153 :	bates	1371	GET_SLOT(x, Matrix_DimNamesSym));
154 :	bates	922	}
155 :
156 :	bates	1448	/* this used to be called sCMatrix_to_gCMatrix(..) [in ./dsCMatrix.c ]: */
157 :	bates	1371	SEXP Csparse_symmetric_to_general(SEXP x)
158 :			{
159 :	maechler	1960	CHM_SP chx = AS_CHM_SP(x), chgx;
160 :	maechler	1736	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
161 :	maechler	1960	R_CheckStack();
162 :	bates	1371
163 :			if (!(chx->stype))
164 :	maechler	1548	error(_("Nonsymmetric matrix in Csparse_symmetric_to_general"));
165 :	maechler	1375	chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
166 :			/* xtype: pattern, "real", complex or .. */
167 :	maechler	1548	return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
168 :	bates	1371	GET_SLOT(x, Matrix_DimNamesSym));
169 :			}
170 :
171 :	maechler	1618	SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo)
172 :	maechler	1598	{
173 :	maechler	1960	CHM_SP chx = AS_CHM_SP(x), chgx;
174 :	maechler	2113	int uploT = (*CHAR(STRING_ELT(uplo,0)) == 'U') ? 1 : -1;
175 :	maechler	1736	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
176 :	maechler	1960	R_CheckStack();
177 :	maechler	1598
178 :	maechler	1618	chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c);
179 :	maechler	1598	/* xtype: pattern, "real", complex or .. */
180 :			return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
181 :			GET_SLOT(x, Matrix_DimNamesSym));
182 :			}
183 :
184 :	bates	1369	SEXP Csparse_transpose(SEXP x, SEXP tri)
185 :	bates	922	{
186 :	maechler	1921	/* TODO: lgCMatrix & igC* currently go via double prec. cholmod -
187 :			* since cholmod (& cs) lacks sparse 'int' matrices */
188 :	maechler	1960	CHM_SP chx = AS_CHM_SP(x);
189 :	maechler	1736	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
190 :	maechler	1960	CHM_SP chxt = cholmod_transpose(chx, chx->xtype, &c);
191 :	bates	1366	SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp;
192 :	bates	1867	int tr = asLogical(tri);
193 :	maechler	1960	R_CheckStack();
194 :	bates	1369
195 :	bates	1366	tmp = VECTOR_ELT(dn, 0); /* swap the dimnames */
196 :			SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1));
197 :			SET_VECTOR_ELT(dn, 1, tmp);
198 :			UNPROTECT(1);
199 :	bates	1867	return chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */
200 :			tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0,
201 :			Rkind, tr ? diag_P(x) : "", dn);
202 :	bates	922	}
203 :
204 :			SEXP Csparse_Csparse_prod(SEXP a, SEXP b)
205 :			{
206 :	maechler	2120	CHM_SP
207 :			cha = AS_CHM_SP(Csparse_diagU2N(a)),
208 :			chb = AS_CHM_SP(Csparse_diagU2N(b)),
209 :	maechler	2125	chc = cholmod_ssmult(cha, chb, /*out_stype:*/ 0,
210 :			cha->xtype, /*out sorted:*/ 1, &c);
211 :			const char *cl_a = class_P(a), *cl_b = class_P(b);
212 :			char diag[] = {'\0', '\0'};
213 :			int uploT = 0;
214 :	bates	1366	SEXP dn = allocVector(VECSXP, 2);
215 :	maechler	1960	R_CheckStack();
216 :	bates	922
217 :	maechler	2125	/* Preserve triangularity and even unit-triangularity if appropriate.
218 :			* Note that in that case, the multiplication itself should happen
219 :			* faster. But there's no support for that in CHOLMOD */
220 :
221 :			/* UGLY hack -- rather should have (fast!) C-level version of
222 :			* is(a, "triangularMatrix") etc */
223 :			if (cl_a[1] == 't' && cl_b[1] == 't')
224 :			/* FIXME: fails for "Cholesky","BunchKaufmann"..*/
225 :			if(*uplo_P(a) == *uplo_P(b)) { /* both upper, or both lower tri. */
226 :			uploT = (*uplo_P(a) == 'U') ? 1 : -1;
227 :			if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */
228 :			/* "remove the diagonal entries": */
229 :			chm_diagN2U(chc, uploT, /* do_realloc */ FALSE);
230 :			diag[0]= 'U';
231 :			}
232 :			else diag[0]= 'N';
233 :			}
234 :	bates	1366	SET_VECTOR_ELT(dn, 0, /* establish dimnames */
235 :			duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
236 :			SET_VECTOR_ELT(dn, 1,
237 :			duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));
238 :	maechler	2125	return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
239 :	bates	922	}
240 :
241 :	maechler	1659	SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans)
242 :	bates	1657	{
243 :	maechler	1659	int tr = asLogical(trans);
244 :	maechler	2120	CHM_SP
245 :			cha = AS_CHM_SP(Csparse_diagU2N(a)),
246 :			chb = AS_CHM_SP(Csparse_diagU2N(b)),
247 :			chTr, chc;
248 :	maechler	2125	const char *cl_a = class_P(a), *cl_b = class_P(b);
249 :			char diag[] = {'\0', '\0'};
250 :			int uploT = 0;
251 :	bates	1657	SEXP dn = allocVector(VECSXP, 2);
252 :	maechler	1960	R_CheckStack();
253 :	bates	1657
254 :	maechler	1960	chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c);
255 :	maechler	1659	chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb,
256 :	maechler	2125	/*out_stype:*/ 0, cha->xtype, /*out sorted:*/ 1, &c);
257 :	maechler	1960	cholmod_free_sparse(&chTr, &c);
258 :	maechler	1659
259 :	maechler	2125	/* Preserve triangularity and unit-triangularity if appropriate;
260 :			* see Csparse_Csparse_prod() for comments */
261 :			if (cl_a[1] == 't' && cl_b[1] == 't')
262 :			if(*uplo_P(a) != *uplo_P(b)) { /* one 'U', the other 'L' */
263 :			uploT = (*uplo_P(b) == 'U') ? 1 : -1;
264 :			if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */
265 :			chm_diagN2U(chc, uploT, /* do_realloc */ FALSE);
266 :			diag[0]= 'U';
267 :			}
268 :			else diag[0]= 'N';
269 :			}
270 :
271 :	bates	1657	SET_VECTOR_ELT(dn, 0, /* establish dimnames */
272 :	maechler	1659	duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), (tr) ? 0 : 1)));
273 :	bates	1657	SET_VECTOR_ELT(dn, 1,
274 :	maechler	1659	duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), (tr) ? 0 : 1)));
275 :	maechler	2125	return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
276 :	bates	1657	}
277 :
278 :	bates	922	SEXP Csparse_dense_prod(SEXP a, SEXP b)
279 :			{
280 :	maechler	2120	CHM_SP cha = AS_CHM_SP(Csparse_diagU2N(a));
281 :	maechler	1660	SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));
282 :	maechler	1960	CHM_DN chb = AS_CHM_DN(b_M);
283 :			CHM_DN chc = cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow,
284 :			chb->xtype, &c);
285 :			SEXP dn = PROTECT(allocVector(VECSXP, 2));
286 :			double one[] = {1,0}, zero[] = {0,0};
287 :			R_CheckStack();
288 :	bates	922
289 :	maechler	1960	cholmod_sdmult(cha, 0, one, zero, chb, chc, &c);
290 :	maechler	1660	SET_VECTOR_ELT(dn, 0, /* establish dimnames */
291 :			duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
292 :			SET_VECTOR_ELT(dn, 1,
293 :			duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1)));
294 :	maechler	1960	UNPROTECT(2);
295 :	maechler	1660	return chm_dense_to_SEXP(chc, 1, 0, dn);
296 :	bates	922	}
297 :	maechler	925
298 :	bates	1067	SEXP Csparse_dense_crossprod(SEXP a, SEXP b)
299 :			{
300 :	maechler	2120	CHM_SP cha = AS_CHM_SP(Csparse_diagU2N(a));
301 :	maechler	1660	SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));
302 :	maechler	1960	CHM_DN chb = AS_CHM_DN(b_M);
303 :			CHM_DN chc = cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol,
304 :			chb->xtype, &c);
305 :			SEXP dn = PROTECT(allocVector(VECSXP, 2));
306 :			double one[] = {1,0}, zero[] = {0,0};
307 :			R_CheckStack();
308 :	bates	1067
309 :	maechler	1960	cholmod_sdmult(cha, 1, one, zero, chb, chc, &c);
310 :	maechler	1660	SET_VECTOR_ELT(dn, 0, /* establish dimnames */
311 :			duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1)));
312 :			SET_VECTOR_ELT(dn, 1,
313 :			duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1)));
314 :	maechler	1960	UNPROTECT(2);
315 :	maechler	1660	return chm_dense_to_SEXP(chc, 1, 0, dn);
316 :	bates	1067	}
317 :
318 :	maechler	2125	/* Computes x'x or x x' -- *also* for Tsparse (triplet = TRUE)
319 :			see Csparse_Csparse_crossprod above for x'y and x y' */
320 :	bates	928	SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet)
321 :	bates	922	{
322 :	maechler	957	int trip = asLogical(triplet),
323 :			tr = asLogical(trans); /* gets reversed because _aat is tcrossprod */
324 :	maechler	2120	CHM_TR cht = trip ? AS_CHM_TR(Tsparse_diagU2N(x)) : (CHM_TR) NULL;
325 :	maechler	1960	CHM_SP chcp, chxt,
326 :	maechler	2120	chx = (trip ?
327 :			cholmod_triplet_to_sparse(cht, cht->nnz, &c) :
328 :			AS_CHM_SP(Csparse_diagU2N(x)));
329 :	bates	1366	SEXP dn = PROTECT(allocVector(VECSXP, 2));
330 :	maechler	1960	R_CheckStack();
331 :	bates	922
332 :	maechler	1960	if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c);
333 :	bates	928	chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);
334 :	maechler	2120	if(!chcp) {
335 :			UNPROTECT(1);
336 :			error(_("Csparse_crossprod(): error return from cholmod_aat()"));
337 :			}
338 :	bates	1360	cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c);
339 :			chcp->stype = 1;
340 :	maechler	1960	if (trip) cholmod_free_sparse(&chx, &c);
341 :	bates	923	if (!tr) cholmod_free_sparse(&chxt, &c);
342 :	maechler	1960	SET_VECTOR_ELT(dn, 0, /* establish dimnames */
343 :	bates	1366	duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym),
344 :	maechler	1660	(tr) ? 0 : 1)));
345 :	bates	1366	SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0)));
346 :			UNPROTECT(1);
347 :	maechler	1548	return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);
348 :	bates	922	}
349 :	bates	923
350 :	maechler	1618	SEXP Csparse_drop(SEXP x, SEXP tol)
351 :			{
352 :	maechler	1960	CHM_SP chx = AS_CHM_SP(x);
353 :			CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c);
354 :	maechler	1618	double dtol = asReal(tol);
355 :	maechler	1736	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
356 :	maechler	1960	R_CheckStack();
357 :	maechler	1618
358 :			if(!cholmod_drop(dtol, ans, &c))
359 :			error(_("cholmod_drop() failed"));
360 :	maechler	1736	return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "",
361 :			GET_SLOT(x, Matrix_DimNamesSym));
362 :	maechler	1618	}
363 :
364 :	bates	1218	SEXP Csparse_horzcat(SEXP x, SEXP y)
365 :			{
366 :	maechler	1960	CHM_SP chx = AS_CHM_SP(x), chy = AS_CHM_SP(y);
367 :	maechler	1548	int Rkind = 0; /* only for "d" - FIXME */
368 :	maechler	1960	R_CheckStack();
369 :	maechler	1375
370 :	bates	1366	/* FIXME: currently drops dimnames */
371 :	maechler	1960	return chm_sparse_to_SEXP(cholmod_horzcat(chx, chy, 1, &c),
372 :			1, 0, Rkind, "", R_NilValue);
373 :	bates	1218	}
374 :
375 :			SEXP Csparse_vertcat(SEXP x, SEXP y)
376 :			{
377 :	maechler	1960	CHM_SP chx = AS_CHM_SP(x), chy = AS_CHM_SP(y);
378 :	maechler	1548	int Rkind = 0; /* only for "d" - FIXME */
379 :	maechler	1960	R_CheckStack();
380 :	maechler	1375
381 :	bates	1366	/* FIXME: currently drops dimnames */
382 :	maechler	1960	return chm_sparse_to_SEXP(cholmod_vertcat(chx, chy, 1, &c),
383 :			1, 0, Rkind, "", R_NilValue);
384 :	bates	1218	}
385 :	bates	1265
386 :			SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2)
387 :			{
388 :	maechler	1960	CHM_SP chx = AS_CHM_SP(x);
389 :	maechler	1736	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
390 :	maechler	1960	CHM_SP ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c);
391 :			R_CheckStack();
392 :	bates	1265
393 :	maechler	1736	return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "",
394 :			GET_SLOT(x, Matrix_DimNamesSym));
395 :	bates	1265	}
396 :	bates	1366
397 :			SEXP Csparse_diagU2N(SEXP x)
398 :			{
399 :	maechler	2120	const char *cl = class_P(x);
400 :			/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
401 :			if (cl[1] != 't' || *diag_P(x) != 'U') {
402 :	maechler	2125	/* "trivially fast" when not triangular (<==> no 'diag' slot),
403 :			or not *unit* triangular */
404 :	maechler	1708	return (x);
405 :			}
406 :	maechler	2125	else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */
407 :	maechler	1960	CHM_SP chx = AS_CHM_SP(x);
408 :			CHM_SP eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c);
409 :	maechler	1708	double one[] = {1, 0};
410 :	maechler	1960	CHM_SP ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c);
411 :	maechler	1710	int uploT = (*uplo_P(x) == 'U') ? 1 : -1;
412 :	maechler	1736	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
413 :	bates	1366
414 :	maechler	1960	R_CheckStack();
415 :			cholmod_free_sparse(&eye, &c);
416 :	maechler	1708	return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N",
417 :	maechler	1736	GET_SLOT(x, Matrix_DimNamesSym));
418 :	maechler	1708	}
419 :	bates	1366	}
420 :
421 :	maechler	2125	SEXP Csparse_diagN2U(SEXP x)
422 :			{
423 :			const char *cl = class_P(x);
424 :			/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
425 :			if (cl[1] != 't' || *diag_P(x) != 'N') {
426 :			/* "trivially fast" when not triangular (<==> no 'diag' slot),
427 :			or already *unit* triangular */
428 :			return (x);
429 :			}
430 :			else { /* triangular with diag='N'): now drop the diagonal */
431 :			/* duplicate, since chx will be modified: */
432 :			CHM_SP chx = AS_CHM_SP(duplicate(x));
433 :			int uploT = (*uplo_P(x) == 'U') ? 1 : -1,
434 :			Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
435 :			R_CheckStack();
436 :
437 :			chm_diagN2U(chx, uploT, /* do_realloc */ FALSE);
438 :
439 :			return chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */,
440 :			uploT, Rkind, "U",
441 :			GET_SLOT(x, Matrix_DimNamesSym));
442 :			}
443 :			}
444 :
445 :	bates	1366	SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j)
446 :			{
447 :	maechler	1960	CHM_SP chx = AS_CHM_SP(x);
448 :	bates	1366	int rsize = (isNull(i)) ? -1 : LENGTH(i),
449 :			csize = (isNull(j)) ? -1 : LENGTH(j);
450 :	maechler	1736	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
451 :	maechler	1960	R_CheckStack();
452 :	bates	1366
453 :			if (rsize >= 0 && !isInteger(i))
454 :			error(_("Index i must be NULL or integer"));
455 :			if (csize >= 0 && !isInteger(j))
456 :			error(_("Index j must be NULL or integer"));
457 :	maechler	1736
458 :	bates	1366	return chm_sparse_to_SEXP(cholmod_submatrix(chx, INTEGER(i), rsize,
459 :	maechler	1375	INTEGER(j), csize,
460 :	bates	1366	TRUE, TRUE, &c),
461 :	maechler	1736	1, 0, Rkind, "",
462 :			/* FIXME: drops dimnames */ R_NilValue);
463 :	bates	1366	}
464 :	bates	2049
465 :			SEXP Csparse_MatrixMarket(SEXP x, SEXP fname)
466 :			{
467 :			FILE *f = fopen(CHAR(asChar(fname)), "w");
468 :
469 :			if (!f)
470 :			error(_("failure to open file \"%s\" for writing"),
471 :			CHAR(asChar(fname)));
472 :	maechler	2120	if (!cholmod_write_sparse(f, AS_CHM_SP(Csparse_diagU2N(x)),
473 :			(CHM_SP)NULL, (char*) NULL, &c))
474 :	bates	2049	error(_("cholmod_write_sparse returned error code"));
475 :			fclose(f);
476 :			return R_NilValue;
477 :			}
478 :	maechler	2137
479 :
480 :			/**
481 :			* Extract the diagonal entries from *triangular* Csparse matrix __or__ a
482 :			* cholmod_sparse factor (LDL = TRUE).
483 :			*
484 :			* @param n dimension of the matrix.
485 :			* @param x_p 'p' (column pointer) slot contents
486 :			* @param x_x 'x' (non-zero entries) slot contents
487 :			* @param perm 'perm' (= permutation vector) slot contents
488 :			* @param resultKind a (SEXP) string indicating which kind of result is desired.
489 :			*
490 :			* @return a SEXP, either a (double) number or a length n-vector of diagonal entries
491 :			*/
492 :			SEXP diag_tC_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind)
493 :			/* ^^^^^^ FIXME[Generalize] to int / ... */
494 :			{
495 :			const char* res_ch = CHAR(STRING_ELT(resultKind,0));
496 :			enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log
497 :			} res_kind = ((!strcmp(res_ch, "trace")) ? trace :
498 :			((!strcmp(res_ch, "sumLog")) ? sum_log :
499 :			((!strcmp(res_ch, "prod")) ? prod :
500 :			((!strcmp(res_ch, "diag")) ? diag :
501 :			((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :
502 :			-1)))));
503 :			int i, n_x, i_from = 0;
504 :			SEXP ans = PROTECT(allocVector(REALSXP,
505 :			/* ^^^^ FIXME[Generalize] */
506 :			(res_kind == diag ||
507 :			res_kind == diag_backpermuted) ? n : 1));
508 :			double *v = REAL(ans);
509 :			/* ^^^^^^ ^^^^ FIXME[Generalize] */
510 :
511 :			#define for_DIAG(v_ASSIGN) \
512 :			for(i = 0; i < n; i++, i_from += n_x) { \
513 :			/* looking at i-th column */ \
514 :			n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \
515 :			v_ASSIGN; \
516 :			}
517 :
518 :			/* NOTA BENE: we assume -- uplo = "L" i.e. lower triangular matrix
519 :			* for uplo = "U" (makes sense with a "dtCMatrix" !),
520 :			* should use x_x[i_from + (nx - 1)] instead of x_x[i_from],
521 :			* where nx = (x_p[i+1] - x_p[i])
522 :			*/
523 :
524 :			switch(res_kind) {
525 :			case trace:
526 :			v[0] = 0.;
527 :			for_DIAG(v[0] += x_x[i_from]);
528 :			break;
529 :
530 :			case sum_log:
531 :			v[0] = 0.;
532 :			for_DIAG(v[0] += log(x_x[i_from]));
533 :			break;
534 :
535 :			case prod:
536 :			v[0] = 1.;
537 :			for_DIAG(v[0] *= x_x[i_from]);
538 :			break;
539 :
540 :			case diag:
541 :			for_DIAG(v[i] = x_x[i_from]);
542 :			break;
543 :
544 :			case diag_backpermuted:
545 :			for_DIAG(v[i] = x_x[i_from]);
546 :			/* now back_permute : */
547 :
548 :			error(_("resultKind = 'diagBack' (back-permuted) is not yet implemented"));
549 :			break;
550 :
551 :			default: /* -1 from above */
552 :			error("diag_tC(): invalid 'resultKind'");
553 :			/* Wall: */ ans = R_NilValue; v = REAL(ans);
554 :			}
555 :
556 :			UNPROTECT(1);
557 :			return ans;
558 :			}
559 :
560 :			/**
561 :			* Extract the diagonal entries from *triangular* Csparse matrix __or__ a
562 :			* cholmod_sparse factor (LDL = TRUE).
563 :			*
564 :			* @param pslot 'p' (column pointer) slot of Csparse matrix/factor
565 :			* @param xslot 'x' (non-zero entries) slot of Csparse matrix/factor
566 :			* @param perm_slot 'perm' (= permutation vector) slot of corresponding CHMfactor
567 :			* @param resultKind a (SEXP) string indicating which kind of result is desired.
568 :			*
569 :			* @return a SEXP, either a (double) number or a length n-vector of diagonal entries
570 :			*/
571 :			SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind)
572 :			{
573 :			int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */
574 :			*x_p = INTEGER(pslot),
575 :			*perm = INTEGER(perm_slot);
576 :			double *x_x = REAL(xslot);
577 :			/* ^^^^^^ ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/
578 :
579 :			return diag_tC_ptr(n, x_p, x_x, perm, resultKind);
580 :			}

---

R-Forge@R-project.org	ViewVC Help
Powered by ViewVC 1.0.0