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

Parent Directory | Revision Log | Patch

	pkg/src/Csparse.c revision 1267, Tue May 16 19:12:13 2006 UTC	pkg/Matrix/src/Csparse.c revision 3072, Fri Mar 27 15:10:48 2015 UTC
#	Line 1	Line 1
1	/* Sparse matrices in compressed column-oriented form */	/* Sparse matrices in compressed column-oriented form */
2
3	#include "Csparse.h"	#include "Csparse.h"
4		#include "Tsparse.h"
5	#include "chm_common.h"	#include "chm_common.h"
6
7	SEXP Csparse_validate(SEXP x)	/** "Cheap" C version of  Csparse_validate() - *not* sorting : */
8		Rboolean isValid_Csparse(SEXP x)
9	{	{
10		/* NB: we do *NOT* check a potential 'x' slot here, at all */
11	SEXP pslot = GET_SLOT(x, Matrix_pSym),	SEXP pslot = GET_SLOT(x, Matrix_pSym),
12	islot = GET_SLOT(x, Matrix_iSym);	islot = GET_SLOT(x, Matrix_iSym);
13	int j, ncol = length(pslot) - 1,	int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), j,
14	*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),	nrow = dims[0],
15	nrow, *xp = INTEGER(pslot),	ncol = dims[1],
16		*xp = INTEGER(pslot),
17	*xi = INTEGER(islot);	*xi = INTEGER(islot);
18
19	nrow = dims[0];	if (length(pslot) != dims[1] + 1)
20	if (length(pslot) <= 0)	return FALSE;
	return mkString(_("slot p must have length > 0"));
21	if (xp[0] != 0)	if (xp[0] != 0)
22	return mkString(_("first element of slot p must be zero"));	return FALSE;
23	if (length(islot) != xp[ncol])	if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/
24	return mkString(_("last element of slot p must match length of slots i and x"));	return FALSE;
25		for (j = 0; j < xp[ncol]; j++) {
26		if (xi[j] < 0 || xi[j] >= nrow)
27		return FALSE;
28		}
29	for (j = 0; j < ncol; j++) {	for (j = 0; j < ncol; j++) {
30	if (xp[j] > xp[j+1])	if (xp[j] > xp[j+1])
31		return FALSE;
32		}
33		return TRUE;
34		}
35
36		SEXP Csparse_validate(SEXP x) {
37		return Csparse_validate_(x, FALSE);
38		}
39
40
41		#define _t_Csparse_validate
42		#include "t_Csparse_validate.c"
43
44		#define _t_Csparse_sort
45		#include "t_Csparse_validate.c"
46
47		// R: .validateCsparse(x, sort.if.needed = FALSE) :
48		SEXP Csparse_validate2(SEXP x, SEXP maybe_modify) {
49		return Csparse_validate_(x, asLogical(maybe_modify));
50		}
51
52		// R: Matrix:::.sortCsparse(x) :
53		SEXP Csparse_sort (SEXP x) {
54		int ok = Csparse_sort_2(x, TRUE); // modifying x directly
55		if(!ok) warning(_("Csparse_sort(x): x is not a valid (apart from sorting) CsparseMatrix"));
56		return x;
57		}
58
59		SEXP Rsparse_validate(SEXP x)
60		{
61		/* NB: we do *NOT* check a potential 'x' slot here, at all */
62		SEXP pslot = GET_SLOT(x, Matrix_pSym),
63		jslot = GET_SLOT(x, Matrix_jSym);
64		Rboolean sorted, strictly;
65		int i, k,
66		*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
67		nrow = dims[0],
68		ncol = dims[1],
69		*xp = INTEGER(pslot),
70		*xj = INTEGER(jslot);
71
72		if (length(pslot) != dims[0] + 1)
73		return mkString(_("slot p must have length = nrow(.) + 1"));
74		if (xp[0] != 0)
75		return mkString(_("first element of slot p must be zero"));
76		if (length(jslot) < xp[nrow]) /* allow larger slots from over-allocation!*/
77		return
78		mkString(_("last element of slot p must match length of slots j and x"));
79		for (i = 0; i < length(jslot); i++) {
80		if (xj[i] < 0 || xj[i] >= ncol)
81		return mkString(_("all column indices must be between 0 and ncol-1"));
82		}
83		sorted = TRUE; strictly = TRUE;
84		for (i = 0; i < nrow; i++) {
85		if (xp[i] > xp[i+1])
86	return mkString(_("slot p must be non-decreasing"));	return mkString(_("slot p must be non-decreasing"));
87		if(sorted)
88		for (k = xp[i] + 1; k < xp[i + 1]; k++) {
89		if (xj[k] < xj[k - 1])
90		sorted = FALSE;
91		else if (xj[k] == xj[k - 1])
92		strictly = FALSE;
93	}	}
	for (j = 0; j < length(islot); j++) {
	if (xi[j] < 0 || xi[j] >= nrow)
	return mkString(_("all row indices must be between 0 and nrow-1"));
94	}	}
95		if (!sorted)
96		/* cannot easily use cholmod_sort(.) ... -> "error out" :*/
97		return mkString(_("slot j is not increasing inside a column"));
98		else if(!strictly) /* sorted, but not strictly */
99		return mkString(_("slot j is not *strictly* increasing inside a column"));
100
101	return ScalarLogical(1);	return ScalarLogical(1);
102	}	}
103
104	SEXP Csparse_to_dense(SEXP x)	/**
105		* From a CsparseMatrix, produce a dense one.
106		* Directly deals with symmetric, triangular and general.
107		* Called from ../R/Csparse.R's  C2dense()
108		*
109		* @param x a CsparseMatrix: currently all 9 of  "[dln][gst]CMatrix"
110		* @param symm_or_tri integer (NA, < 0, > 0, = 0) specifying the knowledge of the caller about x:
111		*      NA  : unknown => will be determined
112		*      = 0 : "generalMatrix" (not symm or tri);
113		*      < 0 : "triangularMatrix"
114		*      > 0 : "symmetricMatrix"
115		*
116		* @return a "denseMatrix"
117		*/
118		SEXP Csparse_to_dense(SEXP x, SEXP symm_or_tri)
119		{
120		Rboolean is_sym, is_tri;
121		int is_sym_or_tri = asInteger(symm_or_tri),
122		ctype = 0; // <- default = "dgC"
123		static const char *valid[] = { MATRIX_VALID_Csparse, ""};
124		if(is_sym_or_tri == NA_INTEGER) { // find if  is(x, "symmetricMatrix") :
125		ctype = Matrix_check_class_etc(x, valid);
126		is_sym = (ctype % 3 == 1);
127		is_tri = (ctype % 3 == 2);
128		} else {
129		is_sym = is_sym_or_tri > 0;
130		is_tri = is_sym_or_tri < 0;
131		// => both are FALSE  iff  is_.. == 0
132		if(is_sym || is_tri)
133		ctype = Matrix_check_class_etc(x, valid);
134		}
135		CHM_SP chxs = AS_CHM_SP__(x);// -> chxs->stype = +- 1 <==> symmetric
136		R_CheckStack();
137		if(is_tri && *diag_P(x) == 'U') { // ==>  x := diagU2N(x), directly for chxs
138		CHM_SP eye = cholmod_speye(chxs->nrow, chxs->ncol, chxs->xtype, &c);
139		double one[] = {1, 0};
140		CHM_SP ans = cholmod_add(chxs, eye, one, one,
141		/* values: */ ((ctype / 3) != 2), // TRUE iff not "nMatrix"
142		TRUE, &c);
143		cholmod_free_sparse(&eye, &c);
144		chxs = cholmod_copy_sparse(ans, &c);
145		cholmod_free_sparse(&ans, &c);
146		}
147		/* The following loses the symmetry property, since cholmod_dense has none,
148		* BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices
149		* to numeric (CHOLMOD_REAL) ones {and we "revert" via chm_dense_to_SEXP()}: */
150		CHM_DN chxd = cholmod_sparse_to_dense(chxs, &c);
151		int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x);
152
153		SEXP ans = chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym),
154		/* transp: */ FALSE);
155		// -> a [dln]geMatrix
156		if(is_sym) { // ==> want  [dln]syMatrix
157		const char cl1 = class_P(ans)[0];
158		PROTECT(ans);
159		SEXP aa = PROTECT(NEW_OBJECT(MAKE_CLASS((cl1 == 'd') ? "dsyMatrix" :
160		((cl1 == 'l') ? "lsyMatrix" : "nsyMatrix"))));
161		// No need to duplicate() as slots of ans are freshly allocated and ans will not be used
162		SET_SLOT(aa, Matrix_xSym,       GET_SLOT(ans, Matrix_xSym));
163		SET_SLOT(aa, Matrix_DimSym,     GET_SLOT(ans, Matrix_DimSym));
164		SET_SLOT(aa, Matrix_DimNamesSym,GET_SLOT(ans, Matrix_DimNamesSym));
165		SET_SLOT(aa, Matrix_uploSym, mkString((chxs->stype > 0) ? "U" : "L"));
166		UNPROTECT(2);
167		return aa;
168		}
169		else if(is_tri) { // ==> want  [dln]trMatrix
170		const char cl1 = class_P(ans)[0];
171		PROTECT(ans);
172		SEXP aa = PROTECT(NEW_OBJECT(MAKE_CLASS((cl1 == 'd') ? "dtrMatrix" :
173		((cl1 == 'l') ? "ltrMatrix" : "ntrMatrix"))));
174		// No need to duplicate() as slots of ans are freshly allocated and ans will not be used
175		SET_SLOT(aa, Matrix_xSym,       GET_SLOT(ans, Matrix_xSym));
176		SET_SLOT(aa, Matrix_DimSym,     GET_SLOT(ans, Matrix_DimSym));
177		SET_SLOT(aa, Matrix_DimNamesSym,GET_SLOT(ans, Matrix_DimNamesSym));
178		slot_dup(aa, x, Matrix_uploSym);
179		/* already by NEW_OBJECT(..) above:
180		SET_SLOT(aa, Matrix_diagSym, mkString("N")); */
181		UNPROTECT(2);
182		return aa;
183		}
184		else
185		return ans;
186		}
187
188		// FIXME: do not go via CHM (should not be too hard, to just *drop* the x-slot, right?
189		SEXP Csparse2nz(SEXP x, Rboolean tri)
190		{
191		CHM_SP chxs = AS_CHM_SP__(x);
192		CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c);
193		R_CheckStack();
194
195		return chm_sparse_to_SEXP(chxcp, 1/*do_free*/,
196		tri ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
197		/* Rkind: pattern */ 0,
198		/* diag = */ tri ? diag_P(x) : "",
199		GET_SLOT(x, Matrix_DimNamesSym));
200		}
201		SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri)
202		{
203		int tr_ = asLogical(tri);
204		if(tr_ == NA_LOGICAL) {
205		warning(_("Csparse_to_nz_pattern(x, tri = NA): 'tri' is taken as TRUE"));
206		tr_ = TRUE;
207		}
208		return Csparse2nz(x, (Rboolean) tr_);
209		}
210
211		// n.CMatrix --> [dli].CMatrix  (not going through CHM!)
212		SEXP nz_pattern_to_Csparse(SEXP x, SEXP res_kind)
213	{	{
214	cholmod_sparse *chxs = as_cholmod_sparse(x);	return nz2Csparse(x, asInteger(res_kind));
215	cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c);	}
216
217		// n.CMatrix --> [dli].CMatrix  (not going through CHM!)
218		// NOTE: use chm_MOD_xtype(() to change type of  'cholmod_sparse' matrix
219		SEXP nz2Csparse(SEXP x, enum x_slot_kind r_kind)
220		{
221		const char *cl_x = class_P(x);
222		if(cl_x[0] != 'n') error(_("not a 'n.CMatrix'"));
223		if(cl_x[2] != 'C') error(_("not a CsparseMatrix"));
224		int nnz = LENGTH(GET_SLOT(x, Matrix_iSym));
225		SEXP ans;
226		char *ncl = alloca(strlen(cl_x) + 1); /* not much memory required */
227		strcpy(ncl, cl_x);
228		double *dx_x; int *ix_x;
229		ncl[0] = (r_kind == x_double ? 'd' :
230		(r_kind == x_logical ? 'l' :
231		/* else (for now):  r_kind == x_integer : */ 'i'));
232		PROTECT(ans = NEW_OBJECT(MAKE_CLASS(ncl)));
233		// create a correct 'x' slot:
234		switch(r_kind) {
235		int i;
236		case x_double: // 'd'
237		dx_x = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz));
238		for (i=0; i < nnz; i++) dx_x[i] = 1.;
239		break;
240		case x_logical: // 'l'
241		ix_x = LOGICAL(ALLOC_SLOT(ans, Matrix_xSym, LGLSXP, nnz));
242		for (i=0; i < nnz; i++) ix_x[i] = TRUE;
243		break;
244		case x_integer: // 'i'
245		ix_x = INTEGER(ALLOC_SLOT(ans, Matrix_xSym, INTSXP, nnz));
246		for (i=0; i < nnz; i++) ix_x[i] = 1;
247		break;
248
249		default:
250		error(_("nz2Csparse(): invalid/non-implemented r_kind = %d"),
251		r_kind);
252		}
253
254		// now copy all other slots :
255		slot_dup(ans, x, Matrix_iSym);
256		slot_dup(ans, x, Matrix_pSym);
257		slot_dup(ans, x, Matrix_DimSym);
258		slot_dup(ans, x, Matrix_DimNamesSym);
259		if(ncl[1] != 'g') { // symmetric or triangular ...
260		slot_dup_if_has(ans, x, Matrix_uploSym);
261		slot_dup_if_has(ans, x, Matrix_diagSym);
262		}
263		UNPROTECT(1);
264		return ans;
265		}
266
267		SEXP Csparse_to_matrix(SEXP x, SEXP chk, SEXP symm)
268		{
269		int is_sym = asLogical(symm);
270		if(is_sym == NA_LOGICAL) { // find if  is(x, "symmetricMatrix") :
271		static const char *valid[] = { MATRIX_VALID_Csparse, ""};
272		int ctype = Matrix_check_class_etc(x, valid);
273		is_sym = (ctype % 3 == 1);
274		}
275		return chm_dense_to_matrix(
276		cholmod_sparse_to_dense(AS_CHM_SP2(x, asLogical(chk)), &c),
277		1 /*do_free*/,
278		(is_sym
279		? symmetric_DimNames(GET_SLOT(x, Matrix_DimNamesSym))
280		:                    GET_SLOT(x, Matrix_DimNamesSym)));
281		}
282
283	Free(chxs);	SEXP Csparse_to_vector(SEXP x)
284	return chm_dense_to_SEXP(chxd, 1);	{
285		return chm_dense_to_vector(cholmod_sparse_to_dense(AS_CHM_SP__(x), &c), 1);
286	}	}
287
288	SEXP Csparse_to_Tsparse(SEXP x)	SEXP Csparse_to_Tsparse(SEXP x, SEXP tri)
289	{	{
290	cholmod_sparse *chxs = as_cholmod_sparse(x);	CHM_SP chxs = AS_CHM_SP__(x);
291	cholmod_triplet *chxt = cholmod_sparse_to_triplet(chxs, &c);	CHM_TR chxt = cholmod_sparse_to_triplet(chxs, &c);
292		int tr = asLogical(tri);
293		int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
294		R_CheckStack();
295
296		return chm_triplet_to_SEXP(chxt, 1,
297		tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
298		Rkind, tr ? diag_P(x) : "",
299		GET_SLOT(x, Matrix_DimNamesSym));
300		}
301
302	Free(chxs);	SEXP Csparse_to_tCsparse(SEXP x, SEXP uplo, SEXP diag)
303	return chm_triplet_to_SEXP(chxt, 1);	{
304		CHM_SP chxs = AS_CHM_SP__(x);
305		int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
306		R_CheckStack();
307		return chm_sparse_to_SEXP(chxs, /* dofree = */ 0,
308		/* uploT = */ (*CHAR(asChar(uplo)) == 'U')? 1: -1,
309		Rkind, /* diag = */ CHAR(STRING_ELT(diag, 0)),
310		GET_SLOT(x, Matrix_DimNamesSym));
311	}	}
312
313	SEXP Csparse_transpose(SEXP x)	SEXP Csparse_to_tTsparse(SEXP x, SEXP uplo, SEXP diag)
314	{	{
315	cholmod_sparse *chx = as_cholmod_sparse(x);	CHM_SP chxs = AS_CHM_SP__(x);
316	cholmod_sparse *chxt = cholmod_transpose(chx, (int) chx->xtype, &c);	CHM_TR chxt = cholmod_sparse_to_triplet(chxs, &c);
317		int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
318		R_CheckStack();
319		return chm_triplet_to_SEXP(chxt, 1,
320		/* uploT = */ (*CHAR(asChar(uplo)) == 'U')? 1: -1,
321		Rkind, /* diag = */ CHAR(STRING_ELT(diag, 0)),
322		GET_SLOT(x, Matrix_DimNamesSym));
323		}
324
325	Free(chx);
326	return chm_sparse_to_SEXP(chxt, 1);	SEXP Csparse_symmetric_to_general(SEXP x)
327		{
328		CHM_SP chx = AS_CHM_SP__(x), chgx;
329		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
330		R_CheckStack();
331
332		if (!(chx->stype))
333		error(_("Nonsymmetric matrix in Csparse_symmetric_to_general"));
334		chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
335		/* xtype: pattern, "real", complex or .. */
336		return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
337		symmetric_DimNames(GET_SLOT(x, Matrix_DimNamesSym)));
338	}	}
339
340	SEXP Csparse_Csparse_prod(SEXP a, SEXP b)	SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo, SEXP sym_dmns)
341	{	{
342	cholmod_sparse *cha = as_cholmod_sparse(a),	int *adims = INTEGER(GET_SLOT(x, Matrix_DimSym)), n = adims[0];
343	*chb = as_cholmod_sparse(b);	if(n != adims[1]) {
344	cholmod_sparse *chc = cholmod_ssmult(cha, chb, 0, cha->xtype, 1, &c);	error(_("Csparse_general_to_symmetric(): matrix is not square!"));
345		return R_NilValue; /* -Wall */
346		}
347		CHM_SP chx = AS_CHM_SP__(x), chgx;
348		int uploT = (*CHAR(asChar(uplo)) == 'U') ? 1 : -1;
349		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
350		R_CheckStack();
351		chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c);
352
353		SEXP dns = GET_SLOT(x, Matrix_DimNamesSym);
354		if(asLogical(sym_dmns))
355		dns = symmetric_DimNames(dns);
356		else if((!isNull(VECTOR_ELT(dns, 0)) &&
357		!isNull(VECTOR_ELT(dns, 1))) ||
358		!isNull(getAttrib(dns, R_NamesSymbol))) {
359		/* symmetrize them if both are not NULL
360		* or names(dimnames(.)) is asymmetric : */
361		dns = PROTECT(duplicate(dns));
362		if(!equal_string_vectors(VECTOR_ELT(dns, 0),
363		VECTOR_ELT(dns, 1))) {
364		if(uploT == 1)
365		SET_VECTOR_ELT(dns, 0, VECTOR_ELT(dns,1));
366		else
367		SET_VECTOR_ELT(dns, 1, VECTOR_ELT(dns,0));
368		}
369		SEXP nms_dns = getAttrib(dns, R_NamesSymbol);
370		if(!isNull(nms_dns) &&  // names(dimnames(.)) :
371		!R_compute_identical(STRING_ELT(nms_dns, 0),
372		STRING_ELT(nms_dns, 1), 16)) {
373		if(uploT == 1)
374		SET_STRING_ELT(nms_dns, 0, STRING_ELT(nms_dns,1));
375		else
376		SET_STRING_ELT(nms_dns, 1, STRING_ELT(nms_dns,0));
377		setAttrib(dns, R_NamesSymbol, nms_dns);
378		}
379		UNPROTECT(1);
380		}
381		/* xtype: pattern, "real", complex or .. */
382		return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", dns);
383		}
384
385	Free(cha); Free(chb);	SEXP Csparse_transpose(SEXP x, SEXP tri)
386	return chm_sparse_to_SEXP(chc, 1);	{
387		/* TODO: lgCMatrix & igC* currently go via double prec. cholmod -
388		*       since cholmod (& cs) lacks sparse 'int' matrices */
389		CHM_SP chx = AS_CHM_SP__(x);
390		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
391		CHM_SP chxt = cholmod_transpose(chx, chx->xtype, &c);
392		SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp;
393		int tr = asLogical(tri);
394		R_CheckStack();
395
396		tmp = VECTOR_ELT(dn, 0);    /* swap the dimnames */
397		SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1));
398		SET_VECTOR_ELT(dn, 1, tmp);
399		if(!isNull(tmp = getAttrib(dn, R_NamesSymbol))) { // swap names(dimnames(.)):
400		SEXP nms_dns = PROTECT(allocVector(VECSXP, 2));
401		SET_VECTOR_ELT(nms_dns, 1, STRING_ELT(tmp, 0));
402		SET_VECTOR_ELT(nms_dns, 0, STRING_ELT(tmp, 1));
403		setAttrib(dn, R_NamesSymbol, nms_dns);
404		UNPROTECT(1);
405	}	}
406		UNPROTECT(1);
407		return chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */
408		tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0,
409		Rkind, tr ? diag_P(x) : "", dn);
410		}
411
412		/* NOTA BENE:  cholmod_ssmult(A,B, ...) ->  ./CHOLMOD/MatrixOps/cholmod_ssmult.c
413		* ---------  computes a patter*n* matrix __always_ when
414		* *one* of A or B is pattern*n*, because of this (line 73-74):
415		---------------------------------------------------------------------------
416		values = values &&
417		(A->xtype != CHOLMOD_PATTERN) && (B->xtype != CHOLMOD_PATTERN) ;
418		---------------------------------------------------------------------------
419		* ==> Often need to copy the patter*n* to a *l*ogical matrix first !!!
420		*/
421
422	SEXP Csparse_dense_prod(SEXP a, SEXP b)	SEXP Csparse_Csparse_prod(SEXP a, SEXP b, SEXP bool_arith)
423	{	{
424	cholmod_sparse *cha = as_cholmod_sparse(a);	CHM_SP
425	cholmod_dense *chb = as_cholmod_dense(b);	cha = AS_CHM_SP(a),
426	cholmod_dense *chc = cholmod_allocate_dense(cha->nrow, chb->ncol,	chb = AS_CHM_SP(b), chc;
427	cha->nrow, chb->xtype, &c);	R_CheckStack();
428	double alpha = 1, beta = 0;	// const char *cl_a = class_P(a), *cl_b = class_P(b);
429		static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };
430		char diag[] = {'\0', '\0'};
431		int uploT = 0, nprot = 1,
432		do_bool = asLogical(bool_arith); // TRUE / NA / FALSE
433		Rboolean
434		a_is_n = (cha->xtype == CHOLMOD_PATTERN),
435		b_is_n = (chb->xtype == CHOLMOD_PATTERN),
436		force_num = (do_bool == FALSE),
437		maybe_bool= (do_bool == NA_LOGICAL);
438
439		#ifdef DEBUG_Matrix_verbose
440		Rprintf("DBG Csparse_C*_prod(%s, %s)\n", class_P(a), class_P(b));
441		#endif
442
443		if(a_is_n && (force_num || (maybe_bool && !b_is_n))) {
444		/* coerce 'a' to  double;
445		* have no CHOLMOD function (pattern -> logical) --> use "our" code */
446		SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;
447		cha = AS_CHM_SP(da);
448		R_CheckStack();
449		a_is_n = FALSE;
450		}
451		else if(b_is_n && (force_num || (maybe_bool && !a_is_n))) {
452		// coerce 'b' to  double
453		SEXP db = PROTECT(nz2Csparse(b, x_double)); nprot++;
454		chb = AS_CHM_SP(db);
455		R_CheckStack();
456		b_is_n = FALSE;
457		}
458		chc = cholmod_ssmult(cha, chb, /*out_stype:*/ 0,
459		/* values : */ do_bool != TRUE,
460		/* sorted = TRUE: */ 1, &c);
461
462		/* Preserve triangularity and even unit-triangularity if appropriate.
463		* Note that in that case, the multiplication itself should happen
464		* faster.  But there's no support for that in CHOLMOD */
465
466		if(Matrix_check_class_etc(a, valid_tri) >= 0 &&
467		Matrix_check_class_etc(b, valid_tri) >= 0)
468		if(*uplo_P(a) == *uplo_P(b)) { /* both upper, or both lower tri. */
469		uploT = (*uplo_P(a) == 'U') ? 1 : -1;
470		if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */
471		/* "remove the diagonal entries": */
472		chm_diagN2U(chc, uploT, /* do_realloc */ FALSE);
473		diag[0]= 'U';
474		}
475		else diag[0]= 'N';
476		}
477
478	cholmod_sdmult(cha, 0, &alpha, &beta, chb, chc, &c);	SEXP dn = PROTECT(allocVector(VECSXP, 2));
479	Free(cha); Free(chb);	SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
480	return chm_dense_to_SEXP(chc, 1);	duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
481		SET_VECTOR_ELT(dn, 1,
482		duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));
483		UNPROTECT(nprot);
484		return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
485	}	}
486
487	SEXP Csparse_dense_crossprod(SEXP a, SEXP b)	/* trans = FALSE:  crossprod(a,b)
488		* trans = TRUE : tcrossprod(a,b) */
489		SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans, SEXP bool_arith)
490	{	{
491	cholmod_sparse *cha = as_cholmod_sparse(a);	int tr = asLogical(trans), nprot = 1,
492	cholmod_dense *chb = as_cholmod_dense(b);	do_bool = asLogical(bool_arith); // TRUE / NA / FALSE
493	cholmod_dense *chc = cholmod_allocate_dense(cha->ncol, chb->ncol,	CHM_SP
494	cha->ncol, chb->xtype, &c);	cha = AS_CHM_SP(a),
495	double alpha = 1, beta = 0;	chb = AS_CHM_SP(b),
496		chTr, chc;
497		R_CheckStack();
498		static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };
499		char diag[] = {'\0', '\0'};
500		int uploT = 0;
501		Rboolean
502		a_is_n = (cha->xtype == CHOLMOD_PATTERN),
503		b_is_n = (chb->xtype == CHOLMOD_PATTERN),
504		force_num = (do_bool == FALSE),
505		maybe_bool= (do_bool == NA_LOGICAL);
506
507		if(a_is_n && (force_num || (maybe_bool && !b_is_n))) {
508		// coerce 'a' to  double
509		SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;
510		cha = AS_CHM_SP(da);
511		R_CheckStack();
512		// a_is_n = FALSE;
513		}
514		else if(b_is_n && (force_num || (maybe_bool && !a_is_n))) {
515		// coerce 'b' to  double
516		SEXP db = PROTECT(nz2Csparse(b, x_double)); nprot++;
517		chb = AS_CHM_SP(db);
518		R_CheckStack();
519		// b_is_n = FALSE;
520		}
521		else if(do_bool == TRUE) { // Want boolean arithmetic: sufficient if *one* is pattern:
522		if(!a_is_n && !b_is_n) {
523		// coerce 'a' to pattern
524		SEXP da = PROTECT(Csparse2nz(a, /* tri = */
525		Matrix_check_class_etc(a, valid_tri) >= 0)); nprot++;
526		cha = AS_CHM_SP(da);
527		R_CheckStack();
528		// a_is_n = TRUE;
529		}
530		}
531		chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c);
532		chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb,
533		/*out_stype:*/ 0, /* values : */ do_bool != TRUE,
534		/* sorted = TRUE: */ 1, &c);
535		cholmod_free_sparse(&chTr, &c);
536
537		/* Preserve triangularity and unit-triangularity if appropriate;
538		* see Csparse_Csparse_prod() for comments */
539		if(Matrix_check_class_etc(a, valid_tri) >= 0 &&
540		Matrix_check_class_etc(b, valid_tri) >= 0)
541		if(*uplo_P(a) != *uplo_P(b)) { /* one 'U', the other 'L' */
542		uploT = (*uplo_P(b) == 'U') ? 1 : -1;
543		if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */
544		chm_diagN2U(chc, uploT, /* do_realloc */ FALSE);
545		diag[0]= 'U';
546		}
547		else diag[0]= 'N';
548		}
549
550	cholmod_sdmult(cha, 1, &alpha, &beta, chb, chc, &c);	SEXP dn = PROTECT(allocVector(VECSXP, 2));
551	Free(cha); Free(chb);	SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
552	return chm_dense_to_SEXP(chc, 1);	duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym),
553		(tr) ? 0 : 1)));
554		SET_VECTOR_ELT(dn, 1,
555		duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym),
556		(tr) ? 0 : 1)));
557		UNPROTECT(nprot);
558		return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
559	}	}
560
561	SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet)	/**
562		* All (dense * sparse)  Matrix products and cross products
563		*
564		*   f( f(<Csparse>)  %*%  f(<dense>) )   for  f in {t(), identity()}
565		*
566		* @param a CsparseMatrix  (n x m)
567		* @param b numeric vector, matrix, or denseMatrix (m x k) or (k x m)  if `transp` is '2' or 'B'
568		* @param transp character.
569		*        = " " : nothing transposed {apart from a}
570		*        = "2" : "transpose 2nd arg": use  t(b) instead of b (= 2nd argument)
571		*        = "c" : "transpose c":       Return  t(c) instead of c
572		*        = "B" : "transpose both":    use t(b) and return t(c) instead of c
573		* NB: For "2", "c", "B", need to transpose a *dense* matrix, B or C --> chm_transpose_dense()
574		*
575		* @return a dense matrix, the matrix product c = g(a,b) :
576		*
577		*                                                Condition (R)   Condition (C)
578		*   R notation            Math notation          cross  transp   t.a t.b t.ans
579		*   ~~~~~~~~~~~~~~~~~     ~~~~~~~~~~~~~~~~~~     ~~~~~~~~~~~~~   ~~~~~~~~~~~~~
580		*   c <-   a %*%   b      C :=      A B            .       " "    .   .   .
581		*   c <-   a %*% t(b)     C :=      A B'           .       "2"    .   |   .
582		*   c <- t(a %*%   b)     C := (A B)'  = B'A'      .       "c"    .   .   |
583		*   c <- t(a %*% t(b))    C := (A B')' = B A'      .       "B"    .   |   |
584		*
585		*   c <-   t(a) %*%   b   C :=      A'B           TRUE     " "    |   .   .
586		*   c <-   t(a) %*% t(b)  C :=      A'B'          TRUE     "2"    |   |   .
587		*   c <- t(t(a) %*%   b)  C := (A'B)'  = B'A      TRUE     "c"    |   .   |
588		*   c <- t(t(a) %*% t(b)) C := (A'B')' = B A      TRUE     "B"    |   |   |
589		*/
590		SEXP Csp_dense_products(SEXP a, SEXP b,
591		Rboolean transp_a, Rboolean transp_b, Rboolean transp_ans)
592	{	{
593	int trip = asLogical(triplet),	CHM_SP cha = AS_CHM_SP(a);
594	tr   = asLogical(trans); /* gets reversed because _aat is tcrossprod */	int a_nc = transp_a ? cha->nrow : cha->ncol,
595	cholmod_triplet	a_nr = transp_a ? cha->ncol : cha->nrow;
596	*cht = trip ? as_cholmod_triplet(x) : (cholmod_triplet*) NULL;	Rboolean
597	cholmod_sparse *chcp, *chxt,	maybe_transp_b = (a_nc == 1),
598	*chx = trip ? cholmod_triplet_to_sparse(cht, cht->nnz, &c)	b_is_vector = FALSE;
599	: as_cholmod_sparse(x);	/* NOTE: trans_b {<--> "use t(b) instead of b" }
600		----  "interferes" with the  case automatic treatment of *vector* b.
601		In that case,  t(b) or b is used "whatever make more sense",
602		according to the general R philosophy of treating vectors in matrix products.
603		*/
604
605		/* repeating a "cheap part" of  mMatrix_as_dgeMatrix2(b, .)  to see if
606		* we have a vector that we might 'transpose_if_vector' : */
607		static const char *valid[] = {"_NOT_A_CLASS_", MATRIX_VALID_ddense, ""};
608		/* int ctype = Matrix_check_class_etc(b, valid);
609		* if (ctype > 0)   /.* a ddenseMatrix object */
610		if (Matrix_check_class_etc(b, valid) < 0) {
611		// not a ddenseM*:  is.matrix() or vector:
612		b_is_vector = !isMatrix(b);
613		}
614
615	if (!tr)	if(b_is_vector) {
616	chxt = cholmod_transpose(chx, (int) chx->xtype, &c);	/* determine *if* we want/need to transpose at all:
617	chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);	* if (length(b) == ncol(A)) have match: use dim = c(n, 1) (<=> do *not* transp);
618	if(!chcp)	*  otherwise, try to transpose: ok  if (ncol(A) == 1) [see also above]:  */
619	error("Csparse_crossprod(): error return from cholmod_aat()");	maybe_transp_b = (LENGTH(b) != a_nc);
620	cholmod_band_inplace((tr) ? -(chcp->nrow) : 0, (tr) ? 0 : chcp->ncol,	// Here, we transpose already in mMatrix_as_dge*()  ==> don't do it later:
621	chcp->xtype, chcp, &c);	transp_b = FALSE;
622	if (trip) {	}
623	cholmod_free_sparse(&chx, &c);	SEXP b_M = PROTECT(mMatrix_as_dgeMatrix2(b, maybe_transp_b));
624	Free(cht);
625	} else {	CHM_DN chb = AS_CHM_DN(b_M), b_t;
626	Free(chx);	R_CheckStack();
627		int ncol_b;
628		if(transp_b) { // transpose b:
629		b_t = cholmod_allocate_dense(chb->ncol, chb->nrow, chb->ncol, chb->xtype, &c);
630		chm_transpose_dense(b_t, chb);
631		ncol_b = b_t->ncol;
632		} else
633		ncol_b = chb->ncol;
634		// Result C {with dim() before it may be transposed}:
635		CHM_DN chc = cholmod_allocate_dense(a_nr, ncol_b, a_nr, chb->xtype, &c);
636		double one[] = {1,0}, zero[] = {0,0};
637		int nprot = 2;
638
639		/* Tim Davis, please FIXME:  currently (2010-11) *fails* when  a  is a pattern matrix:*/
640		if(cha->xtype == CHOLMOD_PATTERN) {
641		/* warning(_("Csparse_dense_prod(): cholmod_sdmult() not yet implemented for pattern./ ngCMatrix" */
642		/*        " --> slightly inefficient coercion")); */
643
644		// This *fails* to produce a CHOLMOD_REAL ..
645		// CHM_SP chd = cholmod_l_copy(cha, cha->stype, CHOLMOD_REAL, &c);
646		// --> use our Matrix-classes
647		SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;
648		cha = AS_CHM_SP(da);
649		}
650
651		/* cholmod_sdmult(A, transp, alpha, beta, X,  Y,  &c): depending on transp == 0 / != 0:
652		*  Y := alpha*(A*X) + beta*Y or alpha*(A'*X) + beta*Y;  here, alpha = 1, beta = 0:
653		*  Y := A*X  or  A'*X
654		*                       NB: always  <sparse> %*% <dense> !
655		*/
656		cholmod_sdmult(cha, transp_a, one, zero, (transp_b ? b_t : chb), /* -> */ chc, &c);
657
658		SEXP dn = PROTECT(allocVector(VECSXP, 2));  /* establish dimnames */
659		SET_VECTOR_ELT(dn, transp_ans ? 1 : 0,
660		duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), transp_a ? 1 : 0)));
661		SET_VECTOR_ELT(dn, transp_ans ? 0 : 1,
662		duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym),
663		transp_b ? 0 : 1)));
664		if(transp_b) cholmod_free_dense(&b_t, &c);
665		UNPROTECT(nprot);
666		return chm_dense_to_SEXP(chc, 1, 0, dn, transp_ans);
667		}
668
669
670		SEXP Csparse_dense_prod(SEXP a, SEXP b, SEXP transp)
671		{
672		return
673		Csp_dense_products(a, b,
674		/* transp_a = */ FALSE,
675		/* transp_b   = */ (*CHAR(asChar(transp)) == '2' || *CHAR(asChar(transp)) == 'B'),
676		/* transp_ans = */ (*CHAR(asChar(transp)) == 'c' || *CHAR(asChar(transp)) == 'B'));
677		}
678
679		SEXP Csparse_dense_crossprod(SEXP a, SEXP b, SEXP transp)
680		{
681		return
682		Csp_dense_products(a, b,
683		/* transp_a = */ TRUE,
684		/* transp_b   = */ (*CHAR(asChar(transp)) == '2' || *CHAR(asChar(transp)) == 'B'),
685		/* transp_ans = */ (*CHAR(asChar(transp)) == 'c' || *CHAR(asChar(transp)) == 'B'));
686		}
687
688
689		/* Computes   x'x  or  x x' -- *also* for Tsparse (triplet = TRUE)
690		see Csparse_Csparse_crossprod above for  x'y and x y' */
691		SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet, SEXP bool_arith)
692		{
693		int tripl = asLogical(triplet),
694		tr   = asLogical(trans), /* gets reversed because _aat is tcrossprod */
695		do_bool = asLogical(bool_arith); // TRUE / NA / FALSE
696		#ifdef AS_CHM_DIAGU2N_FIXED_FINALLY
697		CHM_TR cht = tripl ? AS_CHM_TR(x) : (CHM_TR) NULL;  int nprot = 1;
698		#else /* workaround needed:*/
699		SEXP xx = PROTECT(Tsparse_diagU2N(x));
700		CHM_TR cht = tripl ? AS_CHM_TR__(xx) : (CHM_TR) NULL; int nprot = 2;
701		#endif
702		CHM_SP chcp, chxt,
703		chx = (tripl ?
704		cholmod_triplet_to_sparse(cht, cht->nnz, &c) :
705		AS_CHM_SP(x));
706		SEXP dn = PROTECT(allocVector(VECSXP, 2));
707		R_CheckStack();
708		Rboolean
709		x_is_n = (chx->xtype == CHOLMOD_PATTERN),
710		force_num = (do_bool == FALSE);
711
712		if(x_is_n && force_num) {
713		// coerce 'x' to  double
714		SEXP dx = PROTECT(nz2Csparse(x, x_double)); nprot++;
715		chx = AS_CHM_SP(dx);
716		R_CheckStack();
717		}
718		else if(do_bool == TRUE && !x_is_n) { // Want boolean arithmetic; need patter[n]
719		// coerce 'x' to pattern
720		static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };
721		SEXP dx = PROTECT(Csparse2nz(x, /* tri = */
722		Matrix_check_class_etc(x, valid_tri) >= 0)); nprot++;
723		chx = AS_CHM_SP(dx);
724		R_CheckStack();
725		}
726
727		if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c);
728
729		chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0,
730		/* mode: */ chx->xtype, &c);
731		if(!chcp) {
732		UNPROTECT(1);
733		error(_("Csparse_crossprod(): error return from cholmod_aat()"));
734	}	}
735		cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c);
736		chcp->stype = 1; // symmetric
737		if (tripl) cholmod_free_sparse(&chx, &c);
738	if (!tr) cholmod_free_sparse(&chxt, &c);	if (!tr) cholmod_free_sparse(&chxt, &c);
739	return chm_sparse_to_SEXP(chcp, 1);	SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
740		duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym),
741		(tr) ? 0 : 1)));
742		SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0)));
743		UNPROTECT(nprot);
744		// FIXME: uploT for symmetric ?
745		return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);
746		}
747
748		/* Csparse_drop(x, tol):  drop entries with absolute value < tol, i.e,
749		*  at least all "explicit" zeros */
750		SEXP Csparse_drop(SEXP x, SEXP tol)
751		{
752		const char *cl = class_P(x);
753		/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
754		int tr = (cl[1] == 't');
755		CHM_SP chx = AS_CHM_SP__(x);
756		CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c);
757		double dtol = asReal(tol);
758		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
759		R_CheckStack();
760
761		if(!cholmod_drop(dtol, ans, &c))
762		error(_("cholmod_drop() failed"));
763		return chm_sparse_to_SEXP(ans, 1,
764		tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
765		Rkind, tr ? diag_P(x) : "",
766		GET_SLOT(x, Matrix_DimNamesSym));
767	}	}
768
769	SEXP Csparse_horzcat(SEXP x, SEXP y)	SEXP Csparse_horzcat(SEXP x, SEXP y)
770	{	{
771	cholmod_sparse *chx = as_cholmod_sparse(x),	#define CSPARSE_CAT(_KIND_)                                             \
772	*chy = as_cholmod_sparse(y), *ans;	CHM_SP chx = AS_CHM_SP__(x), chy = AS_CHM_SP__(y);                  \
773		R_CheckStack();                                                     \
774		int Rk_x = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : -3,     \
775		Rk_y = (chy->xtype != CHOLMOD_PATTERN) ? Real_kind(y) : -3, Rkind; \
776		if(Rk_x == -3 || Rk_y == -3) { /* at least one of them is patter"n" */ \
777		if(Rk_x == -3 && Rk_y == -3) { /* fine */                       \
778		} else { /* only one is a patter"n"                             \
779		* "Bug" in cholmod_horzcat()/vertcat(): returns patter"n" matrix if one of them is */ \
780		Rboolean ok;                                                \
781		if(Rk_x == -3) {                                            \
782		ok = chm_MOD_xtype(CHOLMOD_REAL, chx, &c); Rk_x = 0;    \
783		} else if(Rk_y == -3) {                                     \
784		ok = chm_MOD_xtype(CHOLMOD_REAL, chy, &c); Rk_y = 0;    \
785		} else                                                      \
786		error(_("Impossible Rk_x/Rk_y in Csparse_%s(), please report"), _KIND_); \
787		if(!ok)                                                     \
788		error(_("chm_MOD_xtype() was not successful in Csparse_%s(), please report"), \
789		_KIND_);                                          \
790		}                                                               \
791		}                                                                   \
792		Rkind = /* logical if both x and y are */ (Rk_x == 1 && Rk_y == 1) ? 1 : 0
793
794		CSPARSE_CAT("horzcat");
795		// TODO: currently drops dimnames - and we fix at R level;
796
797	ans = cholmod_horzcat(chx, chy, 1, &c);	return chm_sparse_to_SEXP(cholmod_horzcat(chx, chy, 1, &c),
798	Free(chx); Free(chy);	1, 0, Rkind, "", R_NilValue);
	return chm_sparse_to_SEXP(ans, 1);
799	}	}
800
801	SEXP Csparse_vertcat(SEXP x, SEXP y)	SEXP Csparse_vertcat(SEXP x, SEXP y)
802	{	{
803	cholmod_sparse *chx = as_cholmod_sparse(x),	CSPARSE_CAT("vertcat");
804	*chy = as_cholmod_sparse(y), *ans;	// TODO: currently drops dimnames - and we fix at R level;
805
806	ans = cholmod_vertcat(chx, chy, 1, &c);	return chm_sparse_to_SEXP(cholmod_vertcat(chx, chy, 1, &c),
807	Free(chx); Free(chy);	1, 0, Rkind, "", R_NilValue);
	return chm_sparse_to_SEXP(ans, 1);
808	}	}
809
810	SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2)	SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2)
811	{	{
812	cholmod_sparse *chx = as_cholmod_sparse(x), *ans;	CHM_SP chx = AS_CHM_SP__(x);
813		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
814		CHM_SP ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c);
815		R_CheckStack();
816
817		return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "",
818		GET_SLOT(x, Matrix_DimNamesSym));
819		}
820
821		SEXP Csparse_diagU2N(SEXP x)
822		{
823		const char *cl = class_P(x);
824		/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
825		if (cl[1] != 't' || *diag_P(x) != 'U') {
826		/* "trivially fast" when not triangular (<==> no 'diag' slot),
827		or not *unit* triangular */
828		return (x);
829		}
830		else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */
831		CHM_SP chx = AS_CHM_SP__(x);
832		CHM_SP eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c);
833		double one[] = {1, 0};
834		CHM_SP ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c);
835		int uploT = (*uplo_P(x) == 'U') ? 1 : -1;
836		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
837
838		R_CheckStack();
839		cholmod_free_sparse(&eye, &c);
840		return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N",
841		GET_SLOT(x, Matrix_DimNamesSym));
842		}
843		}
844
845		SEXP Csparse_diagN2U(SEXP x)
846		{
847		const char *cl = class_P(x);
848		/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
849		if (cl[1] != 't' || *diag_P(x) != 'N') {
850		/* "trivially fast" when not triangular (<==> no 'diag' slot),
851		or already *unit* triangular */
852		return (x);
853		}
854		else { /* triangular with diag='N'): now drop the diagonal */
855		/* duplicate, since chx will be modified: */
856		SEXP xx = PROTECT(duplicate(x));
857		CHM_SP chx = AS_CHM_SP__(xx);
858		int uploT = (*uplo_P(x) == 'U') ? 1 : -1,
859		Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
860		R_CheckStack();
861
862		chm_diagN2U(chx, uploT, /* do_realloc */ FALSE);
863
864		SEXP ans = chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */,
865		uploT, Rkind, "U",
866		GET_SLOT(x, Matrix_DimNamesSym));
867		UNPROTECT(1);// only now !
868		return ans;
869		}
870		}
871
872		/**
873		* "Indexing" aka subsetting : Compute  x[i,j], also for vectors i and j
874		* Working via CHOLMOD_submatrix, see ./CHOLMOD/MatrixOps/cholmod_submatrix.c
875		* @param x CsparseMatrix
876		* @param i row     indices (0-origin), or NULL (R's)
877		* @param j columns indices (0-origin), or NULL
878		*
879		* @return x[i,j]  still CsparseMatrix --- currently, this loses dimnames
880		*/
881		SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j)
882		{
883		CHM_SP chx = AS_CHM_SP(x); /* << does diagU2N() when needed */
884		int rsize = (isNull(i)) ? -1 : LENGTH(i),
885		csize = (isNull(j)) ? -1 : LENGTH(j);
886		int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
887		R_CheckStack();
888
889		if (rsize >= 0 && !isInteger(i))
890		error(_("Index i must be NULL or integer"));
891		if (csize >= 0 && !isInteger(j))
892		error(_("Index j must be NULL or integer"));
893
894		#define CHM_SUB(_M_, _i_, _j_)                                  \
895		cholmod_submatrix(_M_,                                      \
896		(rsize < 0) ? NULL : INTEGER(_i_), rsize, \
897		(csize < 0) ? NULL : INTEGER(_j_), csize, \
898		TRUE, TRUE, &c)
899		CHM_SP ans;
900		if (!chx->stype) {/* non-symmetric Matrix */
901		ans = CHM_SUB(chx, i, j);
902		}
903		else {
904		/* for now, cholmod_submatrix() only accepts "generalMatrix" */
905		CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
906		ans = CHM_SUB(tmp, i, j);
907		cholmod_free_sparse(&tmp, &c);
908		}
909
910		// "FIXME": currently dropping dimnames, and adding them afterwards in R :
911		/* // dimnames: */
912		/* SEXP x_dns = GET_SLOT(x, Matrix_DimNamesSym), */
913		/*  dn = PROTECT(allocVector(VECSXP, 2)); */
914		return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", /* dimnames: */ R_NilValue);
915		}
916
917		#define _d_Csp_
918		#include "t_Csparse_subassign.c"
919
920		#define _l_Csp_
921		#include "t_Csparse_subassign.c"
922
923		#define _i_Csp_
924		#include "t_Csparse_subassign.c"
925
926		#define _n_Csp_
927		#include "t_Csparse_subassign.c"
928
929		#define _z_Csp_
930		#include "t_Csparse_subassign.c"
931
932
933
934		SEXP Csparse_MatrixMarket(SEXP x, SEXP fname)
935		{
936		FILE *f = fopen(CHAR(asChar(fname)), "w");
937
938		if (!f)
939		error(_("failure to open file \"%s\" for writing"),
940		CHAR(asChar(fname)));
941		if (!cholmod_write_sparse(f, AS_CHM_SP(x),
942		(CHM_SP)NULL, (char*) NULL, &c))
943		error(_("cholmod_write_sparse returned error code"));
944		fclose(f);
945		return R_NilValue;
946		}
947
948
949		/**
950		* Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
951		* cholmod_sparse factor (LDL = TRUE).
952		*
953		* @param n  dimension of the matrix.
954		* @param x_p  'p' (column pointer) slot contents
955		* @param x_x  'x' (non-zero entries) slot contents
956		* @param perm 'perm' (= permutation vector) slot contents; only used for "diagBack"
957		* @param resultKind a (SEXP) string indicating which kind of result is desired.
958		*
959		* @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
960		*/
961		SEXP diag_tC_ptr(int n, int *x_p, double *x_x, Rboolean is_U, int *perm,
962		/*                                ^^^^^^ FIXME[Generalize] to int / ... */
963		SEXP resultKind)
964		{
965		const char* res_ch = CHAR(STRING_ELT(resultKind,0));
966		enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log, min, max, range
967		} res_kind = ((!strcmp(res_ch, "trace")) ? trace :
968		((!strcmp(res_ch, "sumLog")) ? sum_log :
969		((!strcmp(res_ch, "prod")) ? prod :
970		((!strcmp(res_ch, "min")) ? min :
971		((!strcmp(res_ch, "max")) ? max :
972		((!strcmp(res_ch, "range")) ? range :
973		((!strcmp(res_ch, "diag")) ? diag :
974		((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :
975		-1))))))));
976		int i, n_x, i_from;
977		SEXP ans = PROTECT(allocVector(REALSXP,
978		/*                                 ^^^^  FIXME[Generalize] */
979		(res_kind == diag ||
980		res_kind == diag_backpermuted) ? n :
981		(res_kind == range ? 2 : 1)));
982		double *v = REAL(ans);
983		/*  ^^^^^^      ^^^^  FIXME[Generalize] */
984
985		i_from = (is_U ? -1 : 0);
986
987		#define for_DIAG(v_ASSIGN)                                      \
988		for(i = 0; i < n; i++) {                                    \
989		/* looking at i-th column */                            \
990		n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \
991		if( is_U) i_from += n_x;                                \
992		v_ASSIGN;                                               \
993		if(!is_U) i_from += n_x;                                \
994		}
995
996		/* NOTA BENE: we assume  -- uplo = "L" i.e. lower triangular matrix
997		*            for uplo = "U" (makes sense with a "dtCMatrix" !),
998		*            should use  x_x[i_from + (n_x - 1)] instead of x_x[i_from],
999		*            where n_x = (x_p[i+1] - x_p[i])
1000		*/
1001
1002		switch(res_kind) {
1003		case trace: // = sum
1004		v[0] = 0.;
1005		for_DIAG(v[0] += x_x[i_from]);
1006		break;
1007
1008		case sum_log:
1009		v[0] = 0.;
1010		for_DIAG(v[0] += log(x_x[i_from]));
1011		break;
1012
1013		case prod:
1014		v[0] = 1.;
1015		for_DIAG(v[0] *= x_x[i_from]);
1016		break;
1017
1018		case min:
1019		v[0] = R_PosInf;
1020		for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from]);
1021		break;
1022
1023		case max:
1024		v[0] = R_NegInf;
1025		for_DIAG(if(v[0] < x_x[i_from]) v[0] = x_x[i_from]);
1026		break;
1027
1028		case range:
1029		v[0] = R_PosInf;
1030		v[1] = R_NegInf;
1031		for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from];
1032		if(v[1] < x_x[i_from]) v[1] = x_x[i_from]);
1033		break;
1034
1035		case diag:
1036		for_DIAG(v[i] = x_x[i_from]);
1037		break;
1038
1039		case diag_backpermuted:
1040		for_DIAG(v[i] = x_x[i_from]);
1041
1042		warning(_("%s = '%s' (back-permuted) is experimental"),
1043		"resultKind", "diagBack");
1044		/* now back_permute : */
1045		for(i = 0; i < n; i++) {
1046		double tmp = v[i]; v[i] = v[perm[i]]; v[perm[i]] = tmp;
1047		/*^^^^ FIXME[Generalize] */
1048		}
1049		break;
1050
1051		default: /* -1 from above */
1052		error(_("diag_tC(): invalid 'resultKind'"));
1053		/* Wall: */ ans = R_NilValue; v = REAL(ans);
1054		}
1055
1056		UNPROTECT(1);
1057		return ans;
1058		}
1059
1060		/**
1061		* Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
1062		* cholmod_sparse factor (LDL = TRUE).
1063		*
1064		* @param obj -- now a cholmod_sparse factor or a dtCMatrix
1065		* @param pslot  'p' (column pointer)   slot of Csparse matrix/factor
1066		* @param xslot  'x' (non-zero entries) slot of Csparse matrix/factor
1067		* @param perm_slot  'perm' (= permutation vector) slot of corresponding CHMfactor;
1068		*                   only used for "diagBack"
1069		* @param resultKind a (SEXP) string indicating which kind of result is desired.
1070		*
1071		* @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
1072		*/
1073		SEXP diag_tC(SEXP obj, SEXP resultKind)
1074		{
1075
1076		SEXP
1077		pslot = GET_SLOT(obj, Matrix_pSym),
1078		xslot = GET_SLOT(obj, Matrix_xSym);
1079		Rboolean is_U = (R_has_slot(obj, Matrix_uploSym) &&
1080		*CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))) == 'U');
1081		int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */
1082		*x_p  = INTEGER(pslot), pp = -1, *perm;
1083		double *x_x = REAL(xslot);
1084		/*  ^^^^^^        ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/
1085
1086		if(R_has_slot(obj, Matrix_permSym))
1087		perm = INTEGER(GET_SLOT(obj, Matrix_permSym));
1088		else perm = &pp;
1089
1090		return diag_tC_ptr(n, x_p, x_x, is_U, perm, resultKind);
1091		}
1092
1093
1094	ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c);	/**
1095	Free(chx);	* Create a Csparse matrix object from indices and/or pointers.
1096	return chm_sparse_to_SEXP(ans, 1);	*
1097		* @param cls name of actual class of object to create
1098		* @param i optional integer vector of length nnz of row indices
1099		* @param j optional integer vector of length nnz of column indices
1100		* @param p optional integer vector of length np of row or column pointers
1101		* @param np length of integer vector p.  Must be zero if p == (int*)NULL
1102		* @param x optional vector of values
1103		* @param nnz length of vectors i, j and/or x, whichever is to be used
1104		* @param dims optional integer vector of length 2 to be used as
1105		*     dimensions.  If dims == (int*)NULL then the maximum row and column
1106		*     index are used as the dimensions.
1107		* @param dimnames optional list of length 2 to be used as dimnames
1108		* @param index1 indicator of 1-based indices
1109		*
1110		* @return an SEXP of class cls inheriting from CsparseMatrix.
1111		*/
1112		SEXP create_Csparse(char* cls, int* i, int* j, int* p, int np,
1113		void* x, int nnz, int* dims, SEXP dimnames,
1114		int index1)
1115		{
1116		SEXP ans;
1117		int *ij = (int*)NULL, *tri, *trj,
1118		mi, mj, mp, nrow = -1, ncol = -1;
1119		int xtype = -1;             /* -Wall */
1120		CHM_TR T;
1121		CHM_SP A;
1122
1123		if (np < 0 || nnz < 0)
1124		error(_("negative vector lengths not allowed: np = %d, nnz = %d"),
1125		np, nnz);
1126		if (1 != ((mi = (i == (int*)NULL)) +
1127		(mj = (j == (int*)NULL)) +
1128		(mp = (p == (int*)NULL))))
1129		error(_("exactly 1 of 'i', 'j' or 'p' must be NULL"));
1130		if (mp) {
1131		if (np) error(_("np = %d, must be zero when p is NULL"), np);
1132		} else {
1133		if (np) {               /* Expand p to form i or j */
1134		if (!(p[0])) error(_("p[0] = %d, should be zero"), p[0]);
1135		for (int ii = 0; ii < np; ii++)
1136		if (p[ii] > p[ii + 1])
1137		error(_("p must be non-decreasing"));
1138		if (p[np] != nnz)
1139		error("p[np] = %d != nnz = %d", p[np], nnz);
1140		ij = Calloc(nnz, int);
1141		if (mi) {
1142		i = ij;
1143		nrow = np;
1144		} else {
1145		j = ij;
1146		ncol = np;
1147		}
1148		/* Expand p to 0-based indices */
1149		for (int ii = 0; ii < np; ii++)
1150		for (int jj = p[ii]; jj < p[ii + 1]; jj++) ij[jj] = ii;
1151		} else {
1152		if (nnz)
1153		error(_("Inconsistent dimensions: np = 0 and nnz = %d"),
1154		nnz);
1155		}
1156		}
1157		/* calculate nrow and ncol */
1158		if (nrow < 0) {
1159		for (int ii = 0; ii < nnz; ii++) {
1160		int i1 = i[ii] + (index1 ? 0 : 1); /* 1-based index */
1161		if (i1 < 1) error(_("invalid row index at position %d"), ii);
1162		if (i1 > nrow) nrow = i1;
1163		}
1164		}
1165		if (ncol < 0) {
1166		for (int jj = 0; jj < nnz; jj++) {
1167		int j1 = j[jj] + (index1 ? 0 : 1);
1168		if (j1 < 1) error(_("invalid column index at position %d"), jj);
1169		if (j1 > ncol) ncol = j1;
1170		}
1171		}
1172		if (dims != (int*)NULL) {
1173		if (dims[0] > nrow) nrow = dims[0];
1174		if (dims[1] > ncol) ncol = dims[1];
1175		}
1176		/* check the class name */
1177		if (strlen(cls) != 8)
1178		error(_("strlen of cls argument = %d, should be 8"), strlen(cls));
1179		if (!strcmp(cls + 2, "CMatrix"))
1180		error(_("cls = \"%s\" does not end in \"CMatrix\""), cls);
1181		switch(cls[0]) {
1182		case 'd':
1183		case 'l':
1184		xtype = CHOLMOD_REAL;
1185		break;
1186		case 'n':
1187		xtype = CHOLMOD_PATTERN;
1188		break;
1189		default:
1190		error(_("cls = \"%s\" must begin with 'd', 'l' or 'n'"), cls);
1191		}
1192		if (cls[1] != 'g')
1193		error(_("Only 'g'eneral sparse matrix types allowed"));
1194		/* allocate and populate the triplet */
1195		T = cholmod_allocate_triplet((size_t)nrow, (size_t)ncol, (size_t)nnz, 0,
1196		xtype, &c);
1197		T->x = x;
1198		tri = (int*)T->i;
1199		trj = (int*)T->j;
1200		for (int ii = 0; ii < nnz; ii++) {
1201		tri[ii] = i[ii] - ((!mi && index1) ? 1 : 0);
1202		trj[ii] = j[ii] - ((!mj && index1) ? 1 : 0);
1203		}
1204		/* create the cholmod_sparse structure */
1205		A = cholmod_triplet_to_sparse(T, nnz, &c);
1206		cholmod_free_triplet(&T, &c);
1207		/* copy the information to the SEXP */
1208		ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls)));
1209		// FIXME: This has been copied from chm_sparse_to_SEXP in  chm_common.c
1210		/* allocate and copy common slots */
1211		nnz = cholmod_nnz(A, &c);
1212		dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
1213		dims[0] = A->nrow; dims[1] = A->ncol;
1214		Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, A->ncol + 1)), (int*)A->p, A->ncol + 1);
1215		Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)), (int*)A->i, nnz);
1216		switch(cls[1]) {
1217		case 'd':
1218		Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)), (double*)A->x, nnz);
1219		break;
1220		case 'l':
1221		error(_("code not yet written for cls = \"lgCMatrix\""));
1222		}
1223		/* FIXME: dimnames are *NOT* put there yet (if non-NULL) */
1224		cholmod_free_sparse(&A, &c);
1225		UNPROTECT(1);
1226		return ans;
1227	}	}

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

Legend: Removed from v.1267   changed lines   Added in v.3072

---

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