Diff of /pkg/Matrix/src/Csparse.c

Parent Directory | Revision Log | Patch

	pkg/src/Csparse.c revision 930, Mon Sep 19 21:29:37 2005 UTC	pkg/Matrix/src/Csparse.c revision 2586, Sun Jul 25 02:32:06 2010 UTC
#	Line 1	Line 1
1	/* Sparse matrices in compress column-oriented form */	/* Sparse matrices in compressed column-oriented form */
2	#include "Csparse.h"	#include "Csparse.h"
3	#ifdef USE_CHOLMOD	#include "Tsparse.h"
4	#include "chm_common.h"	#include "chm_common.h"
	#endif  /* USE_CHOLMOD */
5
6	SEXP Csparse_validate(SEXP x)	/** "Cheap" C version of  Csparse_validate() - *not* sorting : */
7		Rboolean isValid_Csparse(SEXP x)
8	{	{
9		/* NB: we do *NOT* check a potential 'x' slot here, at all */
10	SEXP pslot = GET_SLOT(x, Matrix_pSym),	SEXP pslot = GET_SLOT(x, Matrix_pSym),
11	islot = GET_SLOT(x, Matrix_iSym);	islot = GET_SLOT(x, Matrix_iSym);
12	int j, ncol = length(pslot) - 1,	int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), j,
13		nrow = dims[0],
14		ncol = dims[1],
15		*xp = INTEGER(pslot),
16		*xi = INTEGER(islot);
17
18		if (length(pslot) != dims[1] + 1)
19		return FALSE;
20		if (xp[0] != 0)
21		return FALSE;
22		if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/
23		return FALSE;
24		for (j = 0; j < xp[ncol]; j++) {
25		if (xi[j] < 0 || xi[j] >= nrow)
26		return FALSE;
27		}
28		for (j = 0; j < ncol; j++) {
29		if (xp[j] > xp[j + 1])
30		return FALSE;
31		}
32		return TRUE;
33		}
34
35		SEXP Csparse_validate(SEXP x) {
36		return Csparse_validate_(x, FALSE);
37		}
38
39		SEXP Csparse_validate2(SEXP x, SEXP maybe_modify) {
40		return Csparse_validate_(x, asLogical(maybe_modify));
41		}
42
43		SEXP Csparse_validate_(SEXP x, Rboolean maybe_modify)
44		{
45		/* NB: we do *NOT* check a potential 'x' slot here, at all */
46		SEXP pslot = GET_SLOT(x, Matrix_pSym),
47		islot = GET_SLOT(x, Matrix_iSym);
48		Rboolean sorted, strictly;
49		int j, k,
50	*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),	*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
51	nrow, *xp = INTEGER(pslot),	nrow = dims[0],
52		ncol = dims[1],
53		*xp = INTEGER(pslot),
54	*xi = INTEGER(islot);	*xi = INTEGER(islot);
55
56	nrow = dims[0];	if (length(pslot) != dims[1] + 1)
57	if (length(pslot) <= 0)	return mkString(_("slot p must have length = ncol(.) + 1"));
	return mkString(_("slot p must have length > 0"));
58	if (xp[0] != 0)	if (xp[0] != 0)
59	return mkString(_("first element of slot p must be zero"));	return mkString(_("first element of slot p must be zero"));
60	if (length(islot) != xp[ncol])	if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/
61	return mkString(_("last element of slot p must match length of slots i and x"));	return
62		mkString(_("last element of slot p must match length of slots i and x"));
63		for (j = 0; j < xp[ncol]; j++) {
64		if (xi[j] < 0 || xi[j] >= nrow)
65		return mkString(_("all row indices must be between 0 and nrow-1"));
66		}
67		sorted = TRUE; strictly = TRUE;
68	for (j = 0; j < ncol; j++) {	for (j = 0; j < ncol; j++) {
69	if (xp[j] > xp[j+1])	if (xp[j] > xp[j+1])
70	return mkString(_("slot p must be non-decreasing"));	return mkString(_("slot p must be non-decreasing"));
71		if(sorted) /* only act if >= 2 entries in column j : */
72		for (k = xp[j] + 1; k < xp[j + 1]; k++) {
73		if (xi[k] < xi[k - 1])
74		sorted = FALSE;
75		else if (xi[k] == xi[k - 1])
76		strictly = FALSE;
77	}	}
78	for (j = 0; j < length(islot); j++) {	}
79	if (xi[j] < 0 || xi[j] >= nrow)	if (!sorted) {
80	return mkString(_("all row indices must be between 0 and nrow-1"));	if(maybe_modify) {
81		CHM_SP chx = (CHM_SP) alloca(sizeof(cholmod_sparse));
82		R_CheckStack();
83		as_cholmod_sparse(chx, x, FALSE, TRUE);/*-> cholmod_l_sort() ! */
84		/* as chx = AS_CHM_SP__(x)  but  ^^^^ sorting x in_place !!! */
85
86		/* Now re-check that row indices are *strictly* increasing
87		* (and not just increasing) within each column : */
88		for (j = 0; j < ncol; j++) {
89		for (k = xp[j] + 1; k < xp[j + 1]; k++)
90		if (xi[k] == xi[k - 1])
91		return mkString(_("slot i is not *strictly* increasing inside a column (even after cholmod_l_sort)"));
92		}
93		} else { /* no modifying sorting : */
94		return mkString(_("row indices are not sorted within columns"));
95		}
96		} else if(!strictly) {  /* sorted, but not strictly */
97		return mkString(_("slot i is not *strictly* increasing inside a column"));
98	}	}
99	return ScalarLogical(1);	return ScalarLogical(1);
100	}	}
101
102	SEXP Csparse_to_Tsparse(SEXP x)	SEXP Rsparse_validate(SEXP x)
103	{	{
104	#ifdef USE_CHOLMOD	/* NB: we do *NOT* check a potential 'x' slot here, at all */
105	cholmod_sparse *chxs = as_cholmod_sparse(x);	SEXP pslot = GET_SLOT(x, Matrix_pSym),
106	cholmod_triplet *chxt = cholmod_sparse_to_triplet(chxs, &c);	jslot = GET_SLOT(x, Matrix_jSym);
107		Rboolean sorted, strictly;
108		int i, k,
109		*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
110		nrow = dims[0],
111		ncol = dims[1],
112		*xp = INTEGER(pslot),
113		*xj = INTEGER(jslot);
114
115	Free(chxs);	if (length(pslot) != dims[0] + 1)
116	return chm_triplet_to_SEXP(chxt, 1);	return mkString(_("slot p must have length = nrow(.) + 1"));
117	#else	if (xp[0] != 0)
118	error("General conversion requires CHOLMOD");	return mkString(_("first element of slot p must be zero"));
119	return R_NilValue;          /* -Wall */	if (length(jslot) < xp[nrow]) /* allow larger slots from over-allocation!*/
120	#endif  /* USE_CHOLMOD */	return
121		mkString(_("last element of slot p must match length of slots j and x"));
122		for (i = 0; i < length(jslot); i++) {
123		if (xj[i] < 0 || xj[i] >= ncol)
124		return mkString(_("all column indices must be between 0 and ncol-1"));
125	}	}
126		sorted = TRUE; strictly = TRUE;
127		for (i = 0; i < nrow; i++) {
128		if (xp[i] > xp[i+1])
129		return mkString(_("slot p must be non-decreasing"));
130		if(sorted)
131		for (k = xp[i] + 1; k < xp[i + 1]; k++) {
132		if (xj[k] < xj[k - 1])
133		sorted = FALSE;
134		else if (xj[k] == xj[k - 1])
135		strictly = FALSE;
136		}
137		}
138		if (!sorted)
139		/* cannot easily use cholmod_l_sort(.) ... -> "error out" :*/
140		return mkString(_("slot j is not increasing inside a column"));
141		else if(!strictly) /* sorted, but not strictly */
142		return mkString(_("slot j is not *strictly* increasing inside a column"));
143
144		return ScalarLogical(1);
145		}
146
147
148	SEXP Csparse_transpose(SEXP x)	/* Called from ../R/Csparse.R : */
149		/* Can only return [dln]geMatrix (no symm/triang);
150		* FIXME: replace by non-CHOLMOD code ! */
151		SEXP Csparse_to_dense(SEXP x)
152	{	{
153	#ifdef USE_CHOLMOD	CHM_SP chxs = AS_CHM_SP__(x);
154	cholmod_sparse *chx = as_cholmod_sparse(x);	/* This loses the symmetry property, since cholmod_dense has none,
155	cholmod_sparse *chxt = cholmod_transpose(chx, (int) chx->xtype, &c);	* BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices
156		* to numeric (CHOLMOD_REAL) ones : */
157		CHM_DN chxd = cholmod_l_sparse_to_dense(chxs, &c);
158		int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x);
159		R_CheckStack();
160
161	Free(chx);	return chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym));
162	return chm_sparse_to_SEXP(chxt, 1);	}
163	#else
164	error("General conversion requires CHOLMOD");	SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri)
165	return R_NilValue;          /* -Wall */	{
166	#endif  /* USE_CHOLMOD */	CHM_SP chxs = AS_CHM_SP__(x);
167		CHM_SP chxcp = cholmod_l_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c);
168		int tr = asLogical(tri);
169		R_CheckStack();
170
171		return chm_sparse_to_SEXP(chxcp, 1/*do_free*/,
172		tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
173		0, tr ? diag_P(x) : "",
174		GET_SLOT(x, Matrix_DimNamesSym));
175		}
176
177		SEXP Csparse_to_matrix(SEXP x)
178		{
179		return chm_dense_to_matrix(cholmod_l_sparse_to_dense(AS_CHM_SP__(x), &c),
180		1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym));
181		}
182
183		SEXP Csparse_to_Tsparse(SEXP x, SEXP tri)
184		{
185		CHM_SP chxs = AS_CHM_SP__(x);
186		CHM_TR chxt = cholmod_l_sparse_to_triplet(chxs, &c);
187		int tr = asLogical(tri);
188		int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
189		R_CheckStack();
190
191		return chm_triplet_to_SEXP(chxt, 1,
192		tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
193		Rkind, tr ? diag_P(x) : "",
194		GET_SLOT(x, Matrix_DimNamesSym));
195	}	}
196
197		/* this used to be called  sCMatrix_to_gCMatrix(..)   [in ./dsCMatrix.c ]: */
198		SEXP Csparse_symmetric_to_general(SEXP x)
199		{
200		CHM_SP chx = AS_CHM_SP__(x), chgx;
201		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
202		R_CheckStack();
203
204		if (!(chx->stype))
205		error(_("Nonsymmetric matrix in Csparse_symmetric_to_general"));
206		chgx = cholmod_l_copy(chx, /* stype: */ 0, chx->xtype, &c);
207		/* xtype: pattern, "real", complex or .. */
208		return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
209		GET_SLOT(x, Matrix_DimNamesSym));
210		}
211
212		SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo)
213		{
214		CHM_SP chx = AS_CHM_SP__(x), chgx;
215		int uploT = (*CHAR(STRING_ELT(uplo,0)) == 'U') ? 1 : -1;
216		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
217		R_CheckStack();
218
219		chgx = cholmod_l_copy(chx, /* stype: */ uploT, chx->xtype, &c);
220		/* xtype: pattern, "real", complex or .. */
221		return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
222		GET_SLOT(x, Matrix_DimNamesSym));
223		}
224
225		SEXP Csparse_transpose(SEXP x, SEXP tri)
226		{
227		/* TODO: lgCMatrix & igC* currently go via double prec. cholmod -
228		*       since cholmod (& cs) lacks sparse 'int' matrices */
229		CHM_SP chx = AS_CHM_SP__(x);
230		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
231		CHM_SP chxt = cholmod_l_transpose(chx, chx->xtype, &c);
232		SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp;
233		int tr = asLogical(tri);
234		R_CheckStack();
235
236		tmp = VECTOR_ELT(dn, 0);    /* swap the dimnames */
237		SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1));
238		SET_VECTOR_ELT(dn, 1, tmp);
239		UNPROTECT(1);
240		return chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */
241		tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0,
242		Rkind, tr ? diag_P(x) : "", dn);
243		}
244
245	SEXP Csparse_Csparse_prod(SEXP a, SEXP b)	SEXP Csparse_Csparse_prod(SEXP a, SEXP b)
246	{	{
247	#ifdef USE_CHOLMOD	CHM_SP
248	cholmod_sparse *cha = as_cholmod_sparse(a), *chb = as_cholmod_sparse(b);	cha = AS_CHM_SP(a),
249	cholmod_sparse *chc = cholmod_ssmult(cha, chb, 0, (int) cha->xtype, 1, &c);	chb = AS_CHM_SP(b),
250		chc = cholmod_l_ssmult(cha, chb, /*out_stype:*/ 0,
251		/* values:= is_numeric (T/F) */ cha->xtype > 0,
252		/*out sorted:*/ 1, &c);
253		const char *cl_a = class_P(a), *cl_b = class_P(b);
254		char diag[] = {'\0', '\0'};
255		int uploT = 0;
256		SEXP dn = PROTECT(allocVector(VECSXP, 2));
257		R_CheckStack();
258
259		#ifdef DEBUG_Matrix_verbose
260		Rprintf("DBG Csparse_C*_prod(%s, %s)\n", cl_a, cl_b);
261		#endif
262
263		/* Preserve triangularity and even unit-triangularity if appropriate.
264		* Note that in that case, the multiplication itself should happen
265		* faster.  But there's no support for that in CHOLMOD */
266
267		/* UGLY hack -- rather should have (fast!) C-level version of
268		*       is(a, "triangularMatrix") etc */
269		if (cl_a[1] == 't' && cl_b[1] == 't')
270		/* FIXME: fails for "Cholesky","BunchKaufmann"..*/
271		if(*uplo_P(a) == *uplo_P(b)) { /* both upper, or both lower tri. */
272		uploT = (*uplo_P(a) == 'U') ? 1 : -1;
273		if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */
274		/* "remove the diagonal entries": */
275		chm_diagN2U(chc, uploT, /* do_realloc */ FALSE);
276		diag[0]= 'U';
277		}
278		else diag[0]= 'N';
279		}
280		SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
281		duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
282		SET_VECTOR_ELT(dn, 1,
283		duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));
284		UNPROTECT(1);
285		return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
286		}
287
288	Free(cha); Free(chb);	SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans)
289	return chm_sparse_to_SEXP(chc, 1);	{
290	#else	int tr = asLogical(trans);
291	error("General multiplication requires CHOLMOD");	CHM_SP
292	return R_NilValue;          /* -Wall */	cha = AS_CHM_SP(a),
293	#endif  /* USE_CHOLMOD */	chb = AS_CHM_SP(b),
294		chTr, chc;
295		const char *cl_a = class_P(a), *cl_b = class_P(b);
296		char diag[] = {'\0', '\0'};
297		int uploT = 0;
298		SEXP dn = PROTECT(allocVector(VECSXP, 2));
299		R_CheckStack();
300
301		chTr = cholmod_l_transpose((tr) ? chb : cha, chb->xtype, &c);
302		chc = cholmod_l_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb,
303		/*out_stype:*/ 0, cha->xtype, /*out sorted:*/ 1, &c);
304		cholmod_l_free_sparse(&chTr, &c);
305
306		/* Preserve triangularity and unit-triangularity if appropriate;
307		* see Csparse_Csparse_prod() for comments */
308		if (cl_a[1] == 't' && cl_b[1] == 't')
309		if(*uplo_P(a) != *uplo_P(b)) { /* one 'U', the other 'L' */
310		uploT = (*uplo_P(b) == 'U') ? 1 : -1;
311		if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */
312		chm_diagN2U(chc, uploT, /* do_realloc */ FALSE);
313		diag[0]= 'U';
314		}
315		else diag[0]= 'N';
316		}
317		SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
318		duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), (tr) ? 0 : 1)));
319		SET_VECTOR_ELT(dn, 1,
320		duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), (tr) ? 0 : 1)));
321		UNPROTECT(1);
322		return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
323	}	}
324
325	SEXP Csparse_dense_prod(SEXP a, SEXP b)	SEXP Csparse_dense_prod(SEXP a, SEXP b)
326	{	{
327	#ifdef USE_CHOLMOD	CHM_SP cha = AS_CHM_SP(a);
328	cholmod_sparse *cha = as_cholmod_sparse(a);	SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));
329	cholmod_dense *chb = as_cholmod_dense(b);	CHM_DN chb = AS_CHM_DN(b_M);
330	cholmod_dense *chc = cholmod_allocate_dense(cha->nrow, chb->ncol,	CHM_DN chc = cholmod_l_allocate_dense(cha->nrow, chb->ncol, cha->nrow,
331	cha->nrow, chb->xtype, &c);	chb->xtype, &c);
332	double alpha = 1, beta = 0;	SEXP dn = PROTECT(allocVector(VECSXP, 2));
333		double one[] = {1,0}, zero[] = {0,0};
334	cholmod_sdmult(cha, 0, &alpha, &beta, chb, chc, &c);	R_CheckStack();
335	Free(cha); Free(chb);
336	return chm_dense_to_SEXP(chc, 1);	cholmod_l_sdmult(cha, 0, one, zero, chb, chc, &c);
337	#else	SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
338	error("General multiplication requires CHOLMOD");	duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
339	return R_NilValue;          /* -Wall */	SET_VECTOR_ELT(dn, 1,
340	#endif  /* USE_CHOLMOD */	duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1)));
341		UNPROTECT(2);
342		return chm_dense_to_SEXP(chc, 1, 0, dn);
343	}	}
344
345		SEXP Csparse_dense_crossprod(SEXP a, SEXP b)
346		{
347		CHM_SP cha = AS_CHM_SP(a);
348		SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));
349		CHM_DN chb = AS_CHM_DN(b_M);
350		CHM_DN chc = cholmod_l_allocate_dense(cha->ncol, chb->ncol, cha->ncol,
351		chb->xtype, &c);
352		SEXP dn = PROTECT(allocVector(VECSXP, 2));
353		double one[] = {1,0}, zero[] = {0,0};
354		R_CheckStack();
355
356		cholmod_l_sdmult(cha, 1, one, zero, chb, chc, &c);
357		SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
358		duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1)));
359		SET_VECTOR_ELT(dn, 1,
360		duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1)));
361		UNPROTECT(2);
362		return chm_dense_to_SEXP(chc, 1, 0, dn);
363		}
364
365		/* Computes   x'x  or  x x' -- *also* for Tsparse (triplet = TRUE)
366		see Csparse_Csparse_crossprod above for  x'y and x y' */
367	SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet)	SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet)
368	{	{
369	#ifdef USE_CHOLMOD	int trip = asLogical(triplet),
370	int trip = asLogical(triplet);	tr   = asLogical(trans); /* gets reversed because _aat is tcrossprod */
371	cholmod_triplet	#ifdef AS_CHM_DIAGU2N_FIXED_FINALLY
372	*cht = trip ? as_cholmod_triplet(x) : (cholmod_triplet*) NULL;	CHM_TR cht = trip ? AS_CHM_TR(x) : (CHM_TR) NULL;
373	cholmod_sparse *chcp, *chxt,	#else /* workaround needed:*/
374	*chx = trip ? cholmod_triplet_to_sparse(cht, cht->nnz, &c)	SEXP xx = PROTECT(Tsparse_diagU2N(x));
375	: as_cholmod_sparse(x);	CHM_TR cht = trip ? AS_CHM_TR__(xx) : (CHM_TR) NULL;
376	int tr = asLogical(trans);  /* gets reversed because _aat is tcrossprod */	#endif
377		CHM_SP chcp, chxt,
378	if (!tr)	chx = (trip ?
379	chxt = cholmod_transpose(chx, (int) chx->xtype, &c);	cholmod_l_triplet_to_sparse(cht, cht->nnz, &c) :
380	chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);	AS_CHM_SP(x));
381		SEXP dn = PROTECT(allocVector(VECSXP, 2));
382	if (trip) {	R_CheckStack();
383	cholmod_free_sparse(&chx, &c);
384	Free(cht);	if (!tr) chxt = cholmod_l_transpose(chx, chx->xtype, &c);
385	} else {	chcp = cholmod_l_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);
386	Free(chx);	if(!chcp) {
387		UNPROTECT(1);
388		error(_("Csparse_crossprod(): error return from cholmod_l_aat()"));
389	}	}
390	if (!tr) cholmod_free_sparse(&chxt, &c);	cholmod_l_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c);
391	return chm_sparse_to_SEXP(chcp, 1);	chcp->stype = 1;
392		if (trip) cholmod_l_free_sparse(&chx, &c);
393		if (!tr) cholmod_l_free_sparse(&chxt, &c);
394		SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
395		duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym),
396		(tr) ? 0 : 1)));
397		SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0)));
398		#ifdef AS_CHM_DIAGU2N_FIXED_FINALLY
399		UNPROTECT(1);
400	#else	#else
401	error("General crossproduct requires CHOLMOD");	UNPROTECT(2);
402	return R_NilValue;          /* -Wall */	#endif
403	#endif  /* USE_CHOLMOD */	return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);
404		}
405
406		SEXP Csparse_drop(SEXP x, SEXP tol)
407		{
408		const char *cl = class_P(x);
409		/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
410		int tr = (cl[1] == 't');
411		CHM_SP chx = AS_CHM_SP__(x);
412		CHM_SP ans = cholmod_l_copy(chx, chx->stype, chx->xtype, &c);
413		double dtol = asReal(tol);
414		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
415		R_CheckStack();
416
417		if(!cholmod_l_drop(dtol, ans, &c))
418		error(_("cholmod_l_drop() failed"));
419		return chm_sparse_to_SEXP(ans, 1,
420		tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
421		Rkind, tr ? diag_P(x) : "",
422		GET_SLOT(x, Matrix_DimNamesSym));
423		}
424
425		SEXP Csparse_horzcat(SEXP x, SEXP y)
426		{
427		CHM_SP chx = AS_CHM_SP__(x), chy = AS_CHM_SP__(y);
428		int Rk_x = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0,
429		Rk_y = (chy->xtype != CHOLMOD_PATTERN) ? Real_kind(y) : 0,
430		Rkind = /* logical if both x and y are */ (Rk_x == 1 && Rk_y == 1) ? 1 : 0;
431		R_CheckStack();
432
433		/* TODO: currently drops dimnames - and we fix at R level */
434		return chm_sparse_to_SEXP(cholmod_l_horzcat(chx, chy, 1, &c),
435		1, 0, Rkind, "", R_NilValue);
436		}
437
438		SEXP Csparse_vertcat(SEXP x, SEXP y)
439		{
440		CHM_SP chx = AS_CHM_SP__(x), chy = AS_CHM_SP__(y);
441		int Rk_x = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0,
442		Rk_y = (chy->xtype != CHOLMOD_PATTERN) ? Real_kind(y) : 0,
443		Rkind = /* logical if both x and y are */ (Rk_x == 1 && Rk_y == 1) ? 1 : 0;
444		R_CheckStack();
445
446		/* TODO: currently drops dimnames - and we fix at R level */
447		return chm_sparse_to_SEXP(cholmod_l_vertcat(chx, chy, 1, &c),
448		1, 0, Rkind, "", R_NilValue);
449		}
450
451		SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2)
452		{
453		CHM_SP chx = AS_CHM_SP__(x);
454		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
455		CHM_SP ans = cholmod_l_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c);
456		R_CheckStack();
457
458		return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "",
459		GET_SLOT(x, Matrix_DimNamesSym));
460		}
461
462		SEXP Csparse_diagU2N(SEXP x)
463		{
464		const char *cl = class_P(x);
465		/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
466		if (cl[1] != 't' || *diag_P(x) != 'U') {
467		/* "trivially fast" when not triangular (<==> no 'diag' slot),
468		or not *unit* triangular */
469		return (x);
470		}
471		else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */
472		CHM_SP chx = AS_CHM_SP__(x);
473		CHM_SP eye = cholmod_l_speye(chx->nrow, chx->ncol, chx->xtype, &c);
474		double one[] = {1, 0};
475		CHM_SP ans = cholmod_l_add(chx, eye, one, one, TRUE, TRUE, &c);
476		int uploT = (*uplo_P(x) == 'U') ? 1 : -1;
477		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
478
479		R_CheckStack();
480		cholmod_l_free_sparse(&eye, &c);
481		return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N",
482		GET_SLOT(x, Matrix_DimNamesSym));
483		}
484		}
485
486		SEXP Csparse_diagN2U(SEXP x)
487		{
488		const char *cl = class_P(x);
489		/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
490		if (cl[1] != 't' || *diag_P(x) != 'N') {
491		/* "trivially fast" when not triangular (<==> no 'diag' slot),
492		or already *unit* triangular */
493		return (x);
494		}
495		else { /* triangular with diag='N'): now drop the diagonal */
496		/* duplicate, since chx will be modified: */
497		CHM_SP chx = AS_CHM_SP__(duplicate(x));
498		int uploT = (*uplo_P(x) == 'U') ? 1 : -1,
499		Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
500		R_CheckStack();
501
502		chm_diagN2U(chx, uploT, /* do_realloc */ FALSE);
503
504		return chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */,
505		uploT, Rkind, "U",
506		GET_SLOT(x, Matrix_DimNamesSym));
507		}
508		}
509
510		/**
511		* "Indexing" aka subsetting : Compute  x[i,j], also for vectors i and j
512		* Working via CHOLMOD_submatrix, see ./CHOLMOD/MatrixOps/cholmod_submatrix.c
513		* @param x CsparseMatrix
514		* @param i row     indices (0-origin), or NULL (R's)
515		* @param j columns indices (0-origin), or NULL
516		*
517		* @return x[i,j]  still CsparseMatrix --- currently, this loses dimnames
518		*/
519		SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j)
520		{
521		CHM_SP chx = AS_CHM_SP(x); /* << does diagU2N() when needed */
522		int rsize = (isNull(i)) ? -1 : LENGTH(i),
523		csize = (isNull(j)) ? -1 : LENGTH(j);
524		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
525		R_CheckStack();
526
527		if (rsize >= 0 && !isInteger(i))
528		error(_("Index i must be NULL or integer"));
529		if (csize >= 0 && !isInteger(j))
530		error(_("Index j must be NULL or integer"));
531
532		if (chx->stype) /* symmetricMatrix */
533		/* for now, cholmod_submatrix() only accepts "generalMatrix" */
534		chx = cholmod_l_copy(chx, /* stype: */ 0, chx->xtype, &c);
535
536		return chm_sparse_to_SEXP(cholmod_l_submatrix(chx,
537		(rsize < 0) ? NULL : INTEGER(i), rsize,
538		(csize < 0) ? NULL : INTEGER(j), csize,
539		TRUE, TRUE, &c),
540		1, 0, Rkind, "",
541		/* FIXME: drops dimnames */ R_NilValue);
542		}
543
544		SEXP Csparse_MatrixMarket(SEXP x, SEXP fname)
545		{
546		FILE *f = fopen(CHAR(asChar(fname)), "w");
547
548		if (!f)
549		error(_("failure to open file \"%s\" for writing"),
550		CHAR(asChar(fname)));
551		if (!cholmod_l_write_sparse(f, AS_CHM_SP(x),
552		(CHM_SP)NULL, (char*) NULL, &c))
553		error(_("cholmod_l_write_sparse returned error code"));
554		fclose(f);
555		return R_NilValue;
556		}
557
558
559		/**
560		* Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
561		* cholmod_sparse factor (LDL = TRUE).
562		*
563		* @param n  dimension of the matrix.
564		* @param x_p  'p' (column pointer) slot contents
565		* @param x_x  'x' (non-zero entries) slot contents
566		* @param perm 'perm' (= permutation vector) slot contents; only used for "diagBack"
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_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind)
572		/*                                ^^^^^^ FIXME[Generalize] to int / ... */
573		{
574		const char* res_ch = CHAR(STRING_ELT(resultKind,0));
575		enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log
576		} res_kind = ((!strcmp(res_ch, "trace")) ? trace :
577		((!strcmp(res_ch, "sumLog")) ? sum_log :
578		((!strcmp(res_ch, "prod")) ? prod :
579		((!strcmp(res_ch, "diag")) ? diag :
580		((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :
581		-1)))));
582		int i, n_x, i_from = 0;
583		SEXP ans = PROTECT(allocVector(REALSXP,
584		/*                                 ^^^^  FIXME[Generalize] */
585		(res_kind == diag ||
586		res_kind == diag_backpermuted) ? n : 1));
587		double *v = REAL(ans);
588		/*  ^^^^^^      ^^^^  FIXME[Generalize] */
589
590		#define for_DIAG(v_ASSIGN)                                              \
591		for(i = 0; i < n; i++, i_from += n_x) {                             \
592		/* looking at i-th column */                                    \
593		n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \
594		v_ASSIGN;                                                       \
595		}
596
597		/* NOTA BENE: we assume  -- uplo = "L" i.e. lower triangular matrix
598		*            for uplo = "U" (makes sense with a "dtCMatrix" !),
599		*            should use  x_x[i_from + (nx - 1)] instead of x_x[i_from],
600		*            where nx = (x_p[i+1] - x_p[i])
601		*/
602
603		switch(res_kind) {
604		case trace:
605		v[0] = 0.;
606		for_DIAG(v[0] += x_x[i_from]);
607		break;
608
609		case sum_log:
610		v[0] = 0.;
611		for_DIAG(v[0] += log(x_x[i_from]));
612		break;
613
614		case prod:
615		v[0] = 1.;
616		for_DIAG(v[0] *= x_x[i_from]);
617		break;
618
619		case diag:
620		for_DIAG(v[i] = x_x[i_from]);
621		break;
622
623		case diag_backpermuted:
624		for_DIAG(v[i] = x_x[i_from]);
625
626		warning(_("resultKind = 'diagBack' (back-permuted) is experimental"));
627		/* now back_permute : */
628		for(i = 0; i < n; i++) {
629		double tmp = v[i]; v[i] = v[perm[i]]; v[perm[i]] = tmp;
630		/*^^^^ FIXME[Generalize] */
631		}
632		break;
633
634		default: /* -1 from above */
635		error(_("diag_tC(): invalid 'resultKind'"));
636		/* Wall: */ ans = R_NilValue; v = REAL(ans);
637		}
638
639		UNPROTECT(1);
640		return ans;
641		}
642
643		/**
644		* Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
645		* cholmod_sparse factor (LDL = TRUE).
646		*
647		* @param pslot  'p' (column pointer)   slot of Csparse matrix/factor
648		* @param xslot  'x' (non-zero entries) slot of Csparse matrix/factor
649		* @param perm_slot  'perm' (= permutation vector) slot of corresponding CHMfactor;
650		*                   only used for "diagBack"
651		* @param resultKind a (SEXP) string indicating which kind of result is desired.
652		*
653		* @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
654		*/
655		SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind)
656		{
657		int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */
658		*x_p  = INTEGER(pslot),
659		*perm = INTEGER(perm_slot);
660		double *x_x = REAL(xslot);
661		/*  ^^^^^^        ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/
662
663		return diag_tC_ptr(n, x_p, x_x, perm, resultKind);
664	}	}
665
666		/**
667		* Create a Csparse matrix object from indices and/or pointers.
668		*
669		* @param cls name of actual class of object to create
670		* @param i optional integer vector of length nnz of row indices
671		* @param j optional integer vector of length nnz of column indices
672		* @param p optional integer vector of length np of row or column pointers
673		* @param np length of integer vector p.  Must be zero if p == (int*)NULL
674		* @param x optional vector of values
675		* @param nnz length of vectors i, j and/or x, whichever is to be used
676		* @param dims optional integer vector of length 2 to be used as
677		*     dimensions.  If dims == (int*)NULL then the maximum row and column
678		*     index are used as the dimensions.
679		* @param dimnames optional list of length 2 to be used as dimnames
680		* @param index1 indicator of 1-based indices
681		*
682		* @return an SEXP of class cls inheriting from CsparseMatrix.
683		*/
684		SEXP create_Csparse(char* cls, int* i, int* j, int* p, int np,
685		void* x, int nnz, int* dims, SEXP dimnames,
686		int index1)
687		{
688		SEXP ans;
689		int *ij = (int*)NULL, *tri, *trj,
690		mi, mj, mp, nrow = -1, ncol = -1;
691		int xtype = -1;             /* -Wall */
692		CHM_TR T;
693		CHM_SP A;
694
695		if (np < 0 || nnz < 0)
696		error(_("negative vector lengths not allowed: np = %d, nnz = %d"),
697		np, nnz);
698		if (1 != ((mi = (i == (int*)NULL)) +
699		(mj = (j == (int*)NULL)) +
700		(mp = (p == (int*)NULL))))
701		error(_("exactly 1 of 'i', 'j' or 'p' must be NULL"));
702		if (mp) {
703		if (np) error(_("np = %d, must be zero when p is NULL"), np);
704		} else {
705		if (np) {               /* Expand p to form i or j */
706		if (!(p[0])) error(_("p[0] = %d, should be zero"), p[0]);
707		for (int ii = 0; ii < np; ii++)
708		if (p[ii] > p[ii + 1])
709		error(_("p must be non-decreasing"));
710		if (p[np] != nnz)
711		error("p[np] = %d != nnz = %d", p[np], nnz);
712		ij = Calloc(nnz, int);
713		if (mi) {
714		i = ij;
715		nrow = np;
716		} else {
717		j = ij;
718		ncol = np;
719		}
720		/* Expand p to 0-based indices */
721		for (int ii = 0; ii < np; ii++)
722		for (int jj = p[ii]; jj < p[ii + 1]; jj++) ij[jj] = ii;
723		} else {
724		if (nnz)
725		error(_("Inconsistent dimensions: np = 0 and nnz = %d"),
726		nnz);
727		}
728		}
729		/* calculate nrow and ncol */
730		if (nrow < 0) {
731		for (int ii = 0; ii < nnz; ii++) {
732		int i1 = i[ii] + (index1 ? 0 : 1); /* 1-based index */
733		if (i1 < 1) error(_("invalid row index at position %d"), ii);
734		if (i1 > nrow) nrow = i1;
735		}
736		}
737		if (ncol < 0) {
738		for (int jj = 0; jj < nnz; jj++) {
739		int j1 = j[jj] + (index1 ? 0 : 1);
740		if (j1 < 1) error(_("invalid column index at position %d"), jj);
741		if (j1 > ncol) ncol = j1;
742		}
743		}
744		if (dims != (int*)NULL) {
745		if (dims[0] > nrow) nrow = dims[0];
746		if (dims[1] > ncol) ncol = dims[1];
747		}
748		/* check the class name */
749		if (strlen(cls) != 8)
750		error(_("strlen of cls argument = %d, should be 8"), strlen(cls));
751		if (!strcmp(cls + 2, "CMatrix"))
752		error(_("cls = \"%s\" does not end in \"CMatrix\""), cls);
753		switch(cls[0]) {
754		case 'd':
755		case 'l':
756		xtype = CHOLMOD_REAL;
757		break;
758		case 'n':
759		xtype = CHOLMOD_PATTERN;
760		break;
761		default:
762		error(_("cls = \"%s\" must begin with 'd', 'l' or 'n'"), cls);
763		}
764		if (cls[1] != 'g')
765		error(_("Only 'g'eneral sparse matrix types allowed"));
766		/* allocate and populate the triplet */
767		T = cholmod_l_allocate_triplet((size_t)nrow, (size_t)ncol, (size_t)nnz, 0,
768		xtype, &c);
769		T->x = x;
770		tri = (int*)T->i;
771		trj = (int*)T->j;
772		for (int ii = 0; ii < nnz; ii++) {
773		tri[ii] = i[ii] - ((!mi && index1) ? 1 : 0);
774		trj[ii] = j[ii] - ((!mj && index1) ? 1 : 0);
775		}
776		/* create the cholmod_sparse structure */
777		A = cholmod_l_triplet_to_sparse(T, nnz, &c);
778		cholmod_l_free_triplet(&T, &c);
779		/* copy the information to the SEXP */
780		ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls)));
781		/* FIXME: This has been copied from chm_sparse_to_SEXP in chm_common.c */
782		/* allocate and copy common slots */
783		nnz = cholmod_l_nnz(A, &c);
784		dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
785		dims[0] = A->nrow; dims[1] = A->ncol;
786		Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, A->ncol + 1)), (int*)A->p, A->ncol + 1);
787		Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)), (int*)A->i, nnz);
788		switch(cls[1]) {
789		case 'd':
790		Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)), (double*)A->x, nnz);
791		break;
792		case 'l':
793		error(_("code not yet written for cls = \"lgCMatrix\""));
794		}
795		cholmod_l_free_sparse(&A, &c);
796		UNPROTECT(1);
797		return ans;
798		}

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

Legend: Removed from v.930   changed lines   Added in v.2586

---

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