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

Parent Directory | Revision Log | Patch

	revision 1657, Wed Nov 1 16:29:53 2006 UTC	revision 2144, Tue Mar 18 23:08:12 2008 UTC
#	Line 1	Line 1
1	/* Sparse matrices in compressed column-oriented form */	/* Sparse matrices in compressed column-oriented form */
2	#include "Csparse.h"	#include "Csparse.h"
3		#include "Tsparse.h"
4	#include "chm_common.h"	#include "chm_common.h"
5
6	SEXP Csparse_validate(SEXP x)	SEXP Csparse_validate(SEXP x)
#	Line 7	Line 8
8	/* NB: we do *NOT* check a potential 'x' slot here, at all */	/* NB: we do *NOT* check a potential 'x' slot here, at all */
9	SEXP pslot = GET_SLOT(x, Matrix_pSym),	SEXP pslot = GET_SLOT(x, Matrix_pSym),
10	islot = GET_SLOT(x, Matrix_iSym);	islot = GET_SLOT(x, Matrix_iSym);
11	int j, k, ncol, nrow, sorted,	Rboolean sorted, strictly;
12		int j, k,
13	*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),	*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
14		nrow = dims[0],
15		ncol = dims[1],
16	*xp = INTEGER(pslot),	*xp = INTEGER(pslot),
17	*xi = INTEGER(islot);	*xi = INTEGER(islot);
18
	nrow = dims[0];
	ncol = dims[1];
19	if (length(pslot) != dims[1] + 1)	if (length(pslot) != dims[1] + 1)
20	return mkString(_("slot p must have length = ncol(.) + 1"));	return mkString(_("slot p must have length = ncol(.) + 1"));
21	if (xp[0] != 0)	if (xp[0] != 0)
22	return mkString(_("first element of slot p must be zero"));	return mkString(_("first element of slot p must be zero"));
23	if (length(islot) != xp[ncol])	if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/
24	return	return
25	mkString(_("last element of slot p must match length of slots i and x"));	mkString(_("last element of slot p must match length of slots i and x"));
26	for (j = 0; j < length(islot); j++) {	for (j = 0; j < length(islot); j++) {
27	if (xi[j] < 0 || xi[j] >= nrow)	if (xi[j] < 0 || xi[j] >= nrow)
28	return mkString(_("all row indices must be between 0 and nrow-1"));	return mkString(_("all row indices must be between 0 and nrow-1"));
29	}	}
30	sorted = TRUE;	sorted = TRUE; strictly = TRUE;
31	for (j = 0; j < ncol; j++) {	for (j = 0; j < ncol; j++) {
32	if (xp[j] > xp[j+1])	if (xp[j] > xp[j+1])
33	return mkString(_("slot p must be non-decreasing"));	return mkString(_("slot p must be non-decreasing"));
34	for (k = xp[j] + 1; k < xp[j + 1]; k++)	if(sorted)
35	if (xi[k] < xi[k - 1]) sorted = FALSE;	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	}	}
42	if (!sorted) {	if (!sorted) {
43	cholmod_sparse *chx = as_cholmod_sparse(x);	CHM_SP chx = AS_CHM_SP(x);
44		R_CheckStack();
45
46	cholmod_sort(chx, &c);	cholmod_sort(chx, &c);
47	Free(chx);	/* 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	}	}
58	return ScalarLogical(1);	return ScalarLogical(1);
59	}	}
60
61		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		/* Called from ../R/Csparse.R : */
108		/* Can only return [dln]geMatrix (no symm/triang);
109		* FIXME: replace by non-CHOLMOD code ! */
110	SEXP Csparse_to_dense(SEXP x)	SEXP Csparse_to_dense(SEXP x)
111	{	{
112	cholmod_sparse *chxs = as_cholmod_sparse(x);	CHM_SP chxs = AS_CHM_SP(x);
113	cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c);	/* 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		CHM_DN chxd = cholmod_sparse_to_dense(chxs, &c);
117		int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x);
118		R_CheckStack();
119
120	Free(chxs);	return chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym));
	return chm_dense_to_SEXP(chxd, 1, Real_kind(x));
121	}	}
122
123	SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri)	SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri)
124	{	{
125	cholmod_sparse *chxs = as_cholmod_sparse(x);	CHM_SP chxs = AS_CHM_SP(x);
126	cholmod_sparse	CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c);
127	*chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c);	int tr = asLogical(tri);
128	int uploT = 0; char *diag = "";	R_CheckStack();
129
130	Free(chxs);	return chm_sparse_to_SEXP(chxcp, 1/*do_free*/,
131	if (asLogical(tri)) {       /* triangular sparse matrices */	tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
132	uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ?	0, tr ? diag_P(x) : "",
	-1 : 1;
	diag = CHAR(asChar(GET_SLOT(x, Matrix_diagSym)));
	}
	return chm_sparse_to_SEXP(chxcp, 1, uploT, 0, diag,
133	GET_SLOT(x, Matrix_DimNamesSym));	GET_SLOT(x, Matrix_DimNamesSym));
134	}	}
135
136	SEXP Csparse_to_matrix(SEXP x)	SEXP Csparse_to_matrix(SEXP x)
137	{	{
138	cholmod_sparse *chxs = as_cholmod_sparse(x);	return chm_dense_to_matrix(cholmod_sparse_to_dense(AS_CHM_SP(x), &c),
139	cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c);	1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym));
	Free(chxs);
	return chm_dense_to_matrix(chxd, 1,
	GET_SLOT(x, Matrix_DimNamesSym));
140	}	}
141
142	SEXP Csparse_to_Tsparse(SEXP x, SEXP tri)	SEXP Csparse_to_Tsparse(SEXP x, SEXP tri)
143	{	{
144	cholmod_sparse *chxs = as_cholmod_sparse(x);	CHM_SP chxs = AS_CHM_SP(x);
145	cholmod_triplet *chxt = cholmod_sparse_to_triplet(chxs, &c);	CHM_TR chxt = cholmod_sparse_to_triplet(chxs, &c);
146	int uploT = 0;	int tr = asLogical(tri);
147	char *diag = "";	int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
148	int Rkind = (chxs->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;	R_CheckStack();
149
150	Free(chxs);	return chm_triplet_to_SEXP(chxt, 1,
151	if (asLogical(tri)) {       /* triangular sparse matrices */	tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
152	uploT = (*uplo_P(x) == 'U') ? -1 : 1;	Rkind, tr ? diag_P(x) : "",
	diag = diag_P(x);
	}
	return chm_triplet_to_SEXP(chxt, 1, uploT, Rkind, diag,
153	GET_SLOT(x, Matrix_DimNamesSym));	GET_SLOT(x, Matrix_DimNamesSym));
154	}	}
155
156	/* this used to be called  sCMatrix_to_gCMatrix(..)   [in ./dsCMatrix.c ]: */	/* this used to be called  sCMatrix_to_gCMatrix(..)   [in ./dsCMatrix.c ]: */
157	SEXP Csparse_symmetric_to_general(SEXP x)	SEXP Csparse_symmetric_to_general(SEXP x)
158	{	{
159	cholmod_sparse *chx = as_cholmod_sparse(x), *chgx;	CHM_SP chx = AS_CHM_SP(x), chgx;
160	int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
161		R_CheckStack();
162
163	if (!(chx->stype))	if (!(chx->stype))
164	error(_("Nonsymmetric matrix in Csparse_symmetric_to_general"));	error(_("Nonsymmetric matrix in Csparse_symmetric_to_general"));
165	chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);	chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
166	/* xtype: pattern, "real", complex or .. */	/* xtype: pattern, "real", complex or .. */
	Free(chx);
167	return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",	return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
168	GET_SLOT(x, Matrix_DimNamesSym));	GET_SLOT(x, Matrix_DimNamesSym));
169	}	}
170
171	SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo)	SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo)
172	{	{
173	cholmod_sparse *chx = as_cholmod_sparse(x), *chgx;	CHM_SP chx = AS_CHM_SP(x), chgx;
174	int uploT = (*CHAR(asChar(uplo)) == 'U') ? -1 : 1;	int uploT = (*CHAR(STRING_ELT(uplo,0)) == 'U') ? 1 : -1;
175	int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
176		R_CheckStack();
177
178	chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c);	chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c);
179	/* xtype: pattern, "real", complex or .. */	/* xtype: pattern, "real", complex or .. */
	Free(chx);
180	return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",	return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
181	GET_SLOT(x, Matrix_DimNamesSym));	GET_SLOT(x, Matrix_DimNamesSym));
182	}	}
183
184	SEXP Csparse_transpose(SEXP x, SEXP tri)	SEXP Csparse_transpose(SEXP x, SEXP tri)
185	{	{
186	cholmod_sparse *chx = as_cholmod_sparse(x);	/* TODO: lgCMatrix & igC* currently go via double prec. cholmod -
187	int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;	*       since cholmod (& cs) lacks sparse 'int' matrices */
188	cholmod_sparse *chxt = cholmod_transpose(chx, (int) chx->xtype, &c);	CHM_SP chx = AS_CHM_SP(x);
189		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
190		CHM_SP chxt = cholmod_transpose(chx, chx->xtype, &c);
191	SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp;	SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp;
192	int uploT = 0; char *diag = "";	int tr = asLogical(tri);
193		R_CheckStack();
194
	Free(chx);
195	tmp = VECTOR_ELT(dn, 0);    /* swap the dimnames */	tmp = VECTOR_ELT(dn, 0);    /* swap the dimnames */
196	SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1));	SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1));
197	SET_VECTOR_ELT(dn, 1, tmp);	SET_VECTOR_ELT(dn, 1, tmp);
198	UNPROTECT(1);	UNPROTECT(1);
199	if (asLogical(tri)) {       /* triangular sparse matrices */	return chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */
200	uploT = (*uplo_P(x) == 'U') ? -1 : 1;	tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0,
201	diag = diag_P(x);	Rkind, tr ? diag_P(x) : "", dn);
	}
	return chm_sparse_to_SEXP(chxt, 1, uploT, Rkind, diag, dn);
202	}	}
203
204	SEXP Csparse_Csparse_prod(SEXP a, SEXP b)	SEXP Csparse_Csparse_prod(SEXP a, SEXP b)
205	{	{
206	cholmod_sparse *cha = as_cholmod_sparse(a),	CHM_SP
207	*chb = as_cholmod_sparse(b);	cha = AS_CHM_SP(Csparse_diagU2N(a)),
208	cholmod_sparse *chc = cholmod_ssmult(cha, chb, 0, cha->xtype, 1, &c);	chb = AS_CHM_SP(Csparse_diagU2N(b)),
209		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	SEXP dn = allocVector(VECSXP, 2);	SEXP dn = allocVector(VECSXP, 2);
215		R_CheckStack();
216
217	Free(cha); Free(chb);	/* 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	SET_VECTOR_ELT(dn, 0,       /* establish dimnames */	SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
235	duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));	duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
236	SET_VECTOR_ELT(dn, 1,	SET_VECTOR_ELT(dn, 1,
237	duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));	duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));
238	return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn);	return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
239	}	}
240
241	SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b)	SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans)
242	{	{
243	cholmod_sparse *cha = as_cholmod_sparse(a),	int tr = asLogical(trans);
244	*chb = as_cholmod_sparse(b);	CHM_SP
245	cholmod_sparse *chta = cholmod_transpose(cha, 1, &c);	cha = AS_CHM_SP(Csparse_diagU2N(a)),
246	cholmod_sparse *chc = cholmod_ssmult(chta, chb, 0, cha->xtype, 1, &c);	chb = AS_CHM_SP(Csparse_diagU2N(b)),
247		chTr, chc;
248		const char *cl_a = class_P(a), *cl_b = class_P(b);
249		char diag[] = {'\0', '\0'};
250		int uploT = 0;
251	SEXP dn = allocVector(VECSXP, 2);	SEXP dn = allocVector(VECSXP, 2);
252		R_CheckStack();
253
254		chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c);
255		chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb,
256		/*out_stype:*/ 0, cha->xtype, /*out sorted:*/ 1, &c);
257		cholmod_free_sparse(&chTr, &c);
258
259		/* 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
	Free(cha); Free(chb); cholmod_free_sparse(&chta, &c);
271	SET_VECTOR_ELT(dn, 0,       /* establish dimnames */	SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
272	duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1)));	duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), (tr) ? 0 : 1)));
273	SET_VECTOR_ELT(dn, 1,	SET_VECTOR_ELT(dn, 1,
274	duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));	duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), (tr) ? 0 : 1)));
275	return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn);	return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
276	}	}
277
278	SEXP Csparse_dense_prod(SEXP a, SEXP b)	SEXP Csparse_dense_prod(SEXP a, SEXP b)
279	{	{
280	cholmod_sparse *cha = as_cholmod_sparse(a);	CHM_SP cha = AS_CHM_SP(Csparse_diagU2N(a));
281	cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b)));	SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));
282	cholmod_dense *chc =	CHM_DN chb = AS_CHM_DN(b_M);
283	cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow, chb->xtype, &c);	CHM_DN chc = cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow,
284	double alpha[] = {1,0}, beta[] = {0,0};	chb->xtype, &c);
285		SEXP dn = PROTECT(allocVector(VECSXP, 2));
286		double one[] = {1,0}, zero[] = {0,0};
287		R_CheckStack();
288
289	cholmod_sdmult(cha, 0, alpha, beta, chb, chc, &c);	cholmod_sdmult(cha, 0, one, zero, chb, chc, &c);
290	Free(cha); Free(chb);	SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
291	UNPROTECT(1);	duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
292	return chm_dense_to_SEXP(chc, 1, 0);	SET_VECTOR_ELT(dn, 1,
293		duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1)));
294		UNPROTECT(2);
295		return chm_dense_to_SEXP(chc, 1, 0, dn);
296	}	}
297
298	SEXP Csparse_dense_crossprod(SEXP a, SEXP b)	SEXP Csparse_dense_crossprod(SEXP a, SEXP b)
299	{	{
300	cholmod_sparse *cha = as_cholmod_sparse(a);	CHM_SP cha = AS_CHM_SP(Csparse_diagU2N(a));
301	cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b)));	SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));
302	cholmod_dense *chc =	CHM_DN chb = AS_CHM_DN(b_M);
303	cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol, chb->xtype, &c);	CHM_DN chc = cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol,
304	double alpha[] = {1,0}, beta[] = {0,0};	chb->xtype, &c);
305		SEXP dn = PROTECT(allocVector(VECSXP, 2));
306		double one[] = {1,0}, zero[] = {0,0};
307		R_CheckStack();
308
309	cholmod_sdmult(cha, 1, alpha, beta, chb, chc, &c);	cholmod_sdmult(cha, 1, one, zero, chb, chc, &c);
310	Free(cha); Free(chb);	SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
311	UNPROTECT(1);	duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1)));
312	return chm_dense_to_SEXP(chc, 1, 0);	SET_VECTOR_ELT(dn, 1,
313		duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1)));
314		UNPROTECT(2);
315		return chm_dense_to_SEXP(chc, 1, 0, dn);
316	}	}
317
318		/* Computes   x'x  or  x x' -- *also* for Tsparse (triplet = TRUE)
319		see Csparse_Csparse_crossprod above for  x'y and x y' */
320	SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet)	SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet)
321	{	{
322	int trip = asLogical(triplet),	int trip = asLogical(triplet),
323	tr   = asLogical(trans); /* gets reversed because _aat is tcrossprod */	tr   = asLogical(trans); /* gets reversed because _aat is tcrossprod */
324	cholmod_triplet	CHM_TR cht = trip ? AS_CHM_TR(Tsparse_diagU2N(x)) : (CHM_TR) NULL;
325	*cht = trip ? as_cholmod_triplet(x) : (cholmod_triplet*) NULL;	CHM_SP chcp, chxt,
326	cholmod_sparse *chcp, *chxt,	chx = (trip ?
327	*chx = trip ? cholmod_triplet_to_sparse(cht, cht->nnz, &c)	cholmod_triplet_to_sparse(cht, cht->nnz, &c) :
328	: as_cholmod_sparse(x);	AS_CHM_SP(Csparse_diagU2N(x)));
329	SEXP dn = PROTECT(allocVector(VECSXP, 2));	SEXP dn = PROTECT(allocVector(VECSXP, 2));
330		R_CheckStack();
331
332	if (!tr)	if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c);
	chxt = cholmod_transpose(chx, chx->xtype, &c);
333	chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);	chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);
334	if(!chcp)	if(!chcp) {
335		UNPROTECT(1);
336	error(_("Csparse_crossprod(): error return from cholmod_aat()"));	error(_("Csparse_crossprod(): error return from cholmod_aat()"));
337		}
338	cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c);	cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c);
339	chcp->stype = 1;	chcp->stype = 1;
340	if (trip) {	if (trip) cholmod_free_sparse(&chx, &c);
	cholmod_free_sparse(&chx, &c);
	Free(cht);
	} else {
	Free(chx);
	}
341	if (!tr) cholmod_free_sparse(&chxt, &c);	if (!tr) cholmod_free_sparse(&chxt, &c);
342	/* create dimnames */	SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
	SET_VECTOR_ELT(dn, 0,
343	duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym),	duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym),
344	(tr) ? 1 : 0)));	(tr) ? 0 : 1)));
345	SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0)));	SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0)));
346	UNPROTECT(1);	UNPROTECT(1);
347	return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);	return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);
#	Line 236	Line 349
349
350	SEXP Csparse_drop(SEXP x, SEXP tol)	SEXP Csparse_drop(SEXP x, SEXP tol)
351	{	{
352	cholmod_sparse *chx = as_cholmod_sparse(x),	CHM_SP chx = AS_CHM_SP(x);
353	*ans = cholmod_copy(chx, chx->stype, chx->xtype, &c);	CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c);
354	double dtol = asReal(tol);	double dtol = asReal(tol);
355	int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
356		R_CheckStack();
357
358	if(!cholmod_drop(dtol, ans, &c))	if(!cholmod_drop(dtol, ans, &c))
359	error(_("cholmod_drop() failed"));	error(_("cholmod_drop() failed"));
360	Free(chx);	return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "",
361	/* FIXME: currently drops dimnames */	GET_SLOT(x, Matrix_DimNamesSym));
	return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue);
362	}	}
363
364	SEXP Csparse_horzcat(SEXP x, SEXP y)	SEXP Csparse_horzcat(SEXP x, SEXP y)
365	{	{
366	cholmod_sparse *chx = as_cholmod_sparse(x),	CHM_SP chx = AS_CHM_SP(x), chy = AS_CHM_SP(y);
	*chy = as_cholmod_sparse(y), *ans;
367	int Rkind = 0; /* only for "d" - FIXME */	int Rkind = 0; /* only for "d" - FIXME */
368		R_CheckStack();
369
	ans = cholmod_horzcat(chx, chy, 1, &c);
	Free(chx); Free(chy);
370	/* FIXME: currently drops dimnames */	/* FIXME: currently drops dimnames */
371	return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue);	return chm_sparse_to_SEXP(cholmod_horzcat(chx, chy, 1, &c),
372		1, 0, Rkind, "", R_NilValue);
373	}	}
374
375	SEXP Csparse_vertcat(SEXP x, SEXP y)	SEXP Csparse_vertcat(SEXP x, SEXP y)
376	{	{
377	cholmod_sparse *chx = as_cholmod_sparse(x),	CHM_SP chx = AS_CHM_SP(x), chy = AS_CHM_SP(y);
	*chy = as_cholmod_sparse(y), *ans;
378	int Rkind = 0; /* only for "d" - FIXME */	int Rkind = 0; /* only for "d" - FIXME */
379		R_CheckStack();
380
	ans = cholmod_vertcat(chx, chy, 1, &c);
	Free(chx); Free(chy);
381	/* FIXME: currently drops dimnames */	/* FIXME: currently drops dimnames */
382	return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue);	return chm_sparse_to_SEXP(cholmod_vertcat(chx, chy, 1, &c),
383		1, 0, Rkind, "", R_NilValue);
384	}	}
385
386	SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2)	SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2)
387	{	{
388	cholmod_sparse *chx = as_cholmod_sparse(x), *ans;	CHM_SP chx = AS_CHM_SP(x);
389	int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
390		CHM_SP ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c);
391		R_CheckStack();
392
393	ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c);	return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "",
394	Free(chx);	GET_SLOT(x, Matrix_DimNamesSym));
	return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue);
395	}	}
396
397	SEXP Csparse_diagU2N(SEXP x)	SEXP Csparse_diagU2N(SEXP x)
398	{	{
399	cholmod_sparse *chx = as_cholmod_sparse(x);	const char *cl = class_P(x);
400	cholmod_sparse *eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c);	/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
401		if (cl[1] != 't' || *diag_P(x) != 'U') {
402		/* "trivially fast" when not triangular (<==> no 'diag' slot),
403		or not *unit* triangular */
404		return (x);
405		}
406		else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */
407		CHM_SP chx = AS_CHM_SP(x);
408		CHM_SP eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c);
409	double one[] = {1, 0};	double one[] = {1, 0};
410	cholmod_sparse *ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c);	CHM_SP ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c);
411	int uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ?	int uploT = (*uplo_P(x) == 'U') ? 1 : -1;
412	-1 : 1;	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
	int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;
413
414	Free(chx); cholmod_free_sparse(&eye, &c);	R_CheckStack();
415		cholmod_free_sparse(&eye, &c);
416	return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N",	return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N",
417	duplicate(GET_SLOT(x, Matrix_DimNamesSym)));	GET_SLOT(x, Matrix_DimNamesSym));
418		}
419		}
420
421		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	SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j)	SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j)
446	{	{
447	cholmod_sparse *chx = as_cholmod_sparse(x);	CHM_SP chx = AS_CHM_SP(x);
448	int rsize = (isNull(i)) ? -1 : LENGTH(i),	int rsize = (isNull(i)) ? -1 : LENGTH(i),
449	csize = (isNull(j)) ? -1 : LENGTH(j);	csize = (isNull(j)) ? -1 : LENGTH(j);
450	int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;	int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
451		R_CheckStack();
452
453	if (rsize >= 0 && !isInteger(i))	if (rsize >= 0 && !isInteger(i))
454	error(_("Index i must be NULL or integer"));	error(_("Index i must be NULL or integer"));
455	if (csize >= 0 && !isInteger(j))	if (csize >= 0 && !isInteger(j))
456	error(_("Index j must be NULL or integer"));	error(_("Index j must be NULL or integer"));
457
458	return chm_sparse_to_SEXP(cholmod_submatrix(chx, INTEGER(i), rsize,	return chm_sparse_to_SEXP(cholmod_submatrix(chx, INTEGER(i), rsize,
459	INTEGER(j), csize,	INTEGER(j), csize,
460	TRUE, TRUE, &c),	TRUE, TRUE, &c),
461	1, 0, Rkind, "", R_NilValue);	1, 0, Rkind, "",
462		/* FIXME: drops dimnames */ R_NilValue);
463		}
464
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		if (!cholmod_write_sparse(f, AS_CHM_SP(Csparse_diagU2N(x)),
473		(CHM_SP)NULL, (char*) NULL, &c))
474		error(_("cholmod_write_sparse returned error code"));
475		fclose(f);
476		return R_NilValue;
477		}
478
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
547		error(_("resultKind = 'diagBack' (back-permuted) is not yet implemented"));
548		/* now back_permute : */
549		for(i = 0; i < n; i++) {
550		double tmp = v[i]; v[i] = v[perm[i]]; v[perm[i]] = tmp;
551		/*^^^^ FIXME[Generalize] */
552		}
553		break;
554
555		default: /* -1 from above */
556		error("diag_tC(): invalid 'resultKind'");
557		/* Wall: */ ans = R_NilValue; v = REAL(ans);
558		}
559
560		UNPROTECT(1);
561		return ans;
562		}
563
564		/**
565		* Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
566		* cholmod_sparse factor (LDL = TRUE).
567		*
568		* @param pslot  'p' (column pointer)   slot of Csparse matrix/factor
569		* @param xslot  'x' (non-zero entries) slot of Csparse matrix/factor
570		* @param perm_slot  'perm' (= permutation vector) slot of corresponding CHMfactor
571		* @param resultKind a (SEXP) string indicating which kind of result is desired.
572		*
573		* @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
574		*/
575		SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind)
576		{
577		int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */
578		*x_p  = INTEGER(pslot),
579		*perm = INTEGER(perm_slot);
580		double *x_x = REAL(xslot);
581		/*  ^^^^^^        ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/
582
583		return diag_tC_ptr(n, x_p, x_x, perm, resultKind);
584	}	}

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

Legend: Removed from v.1657   changed lines   Added in v.2144

---

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