Annotation of /pkg/src/Csparse.c

Parent Directory | Revision Log

---

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

1 :	bates	1218	/* Sparse matrices in compressed column-oriented form */
2 :	bates	922	#include "Csparse.h"
3 :			#include "chm_common.h"
4 :
5 :			SEXP Csparse_validate(SEXP x)
6 :			{
7 :			SEXP pslot = GET_SLOT(x, Matrix_pSym),
8 :			islot = GET_SLOT(x, Matrix_iSym);
9 :			int j, ncol = length(pslot) - 1,
10 :			*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
11 :			nrow, *xp = INTEGER(pslot),
12 :			*xi = INTEGER(islot);
13 :
14 :			nrow = dims[0];
15 :			if (length(pslot) <= 0)
16 :			return mkString(_("slot p must have length > 0"));
17 :			if (xp[0] != 0)
18 :			return mkString(_("first element of slot p must be zero"));
19 :			if (length(islot) != xp[ncol])
20 :			return mkString(_("last element of slot p must match length of slots i and x"));
21 :			for (j = 0; j < ncol; j++) {
22 :			if (xp[j] > xp[j+1])
23 :			return mkString(_("slot p must be non-decreasing"));
24 :			}
25 :			for (j = 0; j < length(islot); j++) {
26 :			if (xi[j] < 0 || xi[j] >= nrow)
27 :			return mkString(_("all row indices must be between 0 and nrow-1"));
28 :			}
29 :			return ScalarLogical(1);
30 :			}
31 :
32 :	bates	1059	SEXP Csparse_to_dense(SEXP x)
33 :			{
34 :			cholmod_sparse *chxs = as_cholmod_sparse(x);
35 :			cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c);
36 :
37 :	bates	1141	Free(chxs);
38 :	bates	1059	return chm_dense_to_SEXP(chxd, 1);
39 :			}
40 :
41 :	bates	1371	SEXP Csparse_to_logical(SEXP x, SEXP tri)
42 :			{
43 :			cholmod_sparse *chxs = as_cholmod_sparse(x);
44 :			cholmod_sparse
45 :			*chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c);
46 :			int uploT = 0; char *diag = "";
47 :
48 :			Free(chxs);
49 :			if (asLogical(tri)) { /* triangular sparse matrices */
50 :			uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ?
51 :			-1 : 1;
52 :			diag = CHAR(asChar(GET_SLOT(x, Matrix_diagSym)));
53 :			}
54 :			return chm_sparse_to_SEXP(chxcp, 1, uploT, diag,
55 :			GET_SLOT(x, Matrix_DimNamesSym));
56 :			}
57 :
58 :	bates	1366	SEXP Csparse_to_matrix(SEXP x)
59 :	bates	922	{
60 :	maechler	925	cholmod_sparse *chxs = as_cholmod_sparse(x);
61 :	bates	1366	cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c);
62 :
63 :			Free(chxs);
64 :			return chm_dense_to_matrix(chxd, 1,
65 :			GET_SLOT(x, Matrix_DimNamesSym));
66 :			}
67 :
68 :			SEXP Csparse_to_Tsparse(SEXP x, SEXP tri)
69 :			{
70 :			cholmod_sparse *chxs = as_cholmod_sparse(x);
71 :	bates	922	cholmod_triplet *chxt = cholmod_sparse_to_triplet(chxs, &c);
72 :	bates	1366	int uploT = 0; char *diag = "";
73 :	bates	922
74 :	bates	1141	Free(chxs);
75 :	bates	1366	if (asLogical(tri)) { /* triangular sparse matrices */
76 :	bates	1419	uploT = (*uplo_P(x) == 'U') ? -1 : 1;
77 :			diag = diag_P(x);
78 :	bates	1366	}
79 :			return chm_triplet_to_SEXP(chxt, 1, uploT, diag,
80 :	bates	1371	GET_SLOT(x, Matrix_DimNamesSym));
81 :	bates	922	}
82 :
83 :	bates	1448	/* this used to be called sCMatrix_to_gCMatrix(..) [in ./dsCMatrix.c ]: */
84 :	bates	1371	SEXP Csparse_symmetric_to_general(SEXP x)
85 :			{
86 :			cholmod_sparse *chx = as_cholmod_sparse(x), *chgx;
87 :
88 :			if (!(chx->stype))
89 :			error(_("Nonsymmetric matrix in Csparse_symmeteric_to_general"));
90 :	maechler	1375	chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
91 :			/* xtype: pattern, "real", complex or .. */
92 :	bates	1371	Free(chx);
93 :			return chm_sparse_to_SEXP(chgx, 1, 0, "",
94 :			GET_SLOT(x, Matrix_DimNamesSym));
95 :			}
96 :
97 :	bates	1369	SEXP Csparse_transpose(SEXP x, SEXP tri)
98 :	bates	922	{
99 :			cholmod_sparse *chx = as_cholmod_sparse(x);
100 :			cholmod_sparse *chxt = cholmod_transpose(chx, (int) chx->xtype, &c);
101 :	bates	1366	SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp;
102 :	bates	1369	int uploT = 0; char *diag = "";
103 :
104 :	bates	1141	Free(chx);
105 :	bates	1366	tmp = VECTOR_ELT(dn, 0); /* swap the dimnames */
106 :			SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1));
107 :			SET_VECTOR_ELT(dn, 1, tmp);
108 :			UNPROTECT(1);
109 :	bates	1369	if (asLogical(tri)) { /* triangular sparse matrices */
110 :	bates	1419	uploT = (*uplo_P(x) == 'U') ? -1 : 1;
111 :			diag = diag_P(x);
112 :	bates	1369	}
113 :			return chm_sparse_to_SEXP(chxt, 1, uploT, diag, dn);
114 :	bates	922	}
115 :
116 :			SEXP Csparse_Csparse_prod(SEXP a, SEXP b)
117 :			{
118 :	bates	1059	cholmod_sparse *cha = as_cholmod_sparse(a),
119 :			*chb = as_cholmod_sparse(b);
120 :			cholmod_sparse *chc = cholmod_ssmult(cha, chb, 0, cha->xtype, 1, &c);
121 :	bates	1366	SEXP dn = allocVector(VECSXP, 2);
122 :	bates	922
123 :	bates	1141	Free(cha); Free(chb);
124 :	bates	1366	SET_VECTOR_ELT(dn, 0, /* establish dimnames */
125 :			duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
126 :			SET_VECTOR_ELT(dn, 1,
127 :			duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));
128 :			return chm_sparse_to_SEXP(chc, 1, 0, "", dn);
129 :	bates	922	}
130 :
131 :			SEXP Csparse_dense_prod(SEXP a, SEXP b)
132 :			{
133 :			cholmod_sparse *cha = as_cholmod_sparse(a);
134 :	bates	1461	cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b)));
135 :	bates	1448	cholmod_dense *chc =
136 :			cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow, chb->xtype, &c);
137 :			double alpha[] = {1,0}, beta[] = {0,0};
138 :	bates	922
139 :	bates	1448	cholmod_sdmult(cha, 0, alpha, beta, chb, chc, &c);
140 :	bates	1141	Free(cha); Free(chb);
141 :	bates	1461	UNPROTECT(1);
142 :	bates	923	return chm_dense_to_SEXP(chc, 1);
143 :	bates	922	}
144 :	maechler	925
145 :	bates	1067	SEXP Csparse_dense_crossprod(SEXP a, SEXP b)
146 :			{
147 :			cholmod_sparse *cha = as_cholmod_sparse(a);
148 :	bates	1461	cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b)));
149 :	bates	1448	cholmod_dense *chc =
150 :			cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol, chb->xtype, &c);
151 :			double alpha[] = {1,0}, beta[] = {0,0};
152 :	bates	1067
153 :	bates	1448	cholmod_sdmult(cha, 1, alpha, beta, chb, chc, &c);
154 :	bates	1067	Free(cha); Free(chb);
155 :	bates	1461	UNPROTECT(1);
156 :	bates	1067	return chm_dense_to_SEXP(chc, 1);
157 :			}
158 :
159 :	bates	928	SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet)
160 :	bates	922	{
161 :	maechler	957	int trip = asLogical(triplet),
162 :			tr = asLogical(trans); /* gets reversed because _aat is tcrossprod */
163 :	bates	928	cholmod_triplet
164 :			*cht = trip ? as_cholmod_triplet(x) : (cholmod_triplet*) NULL;
165 :			cholmod_sparse *chcp, *chxt,
166 :			*chx = trip ? cholmod_triplet_to_sparse(cht, cht->nnz, &c)
167 :			: as_cholmod_sparse(x);
168 :	bates	1366	SEXP dn = PROTECT(allocVector(VECSXP, 2));
169 :	bates	922
170 :	bates	923	if (!tr)
171 :	bates	1353	chxt = cholmod_transpose(chx, chx->xtype, &c);
172 :	bates	928	chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);
173 :	maechler	957	if(!chcp)
174 :			error("Csparse_crossprod(): error return from cholmod_aat()");
175 :	bates	1360	cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c);
176 :			chcp->stype = 1;
177 :	bates	930	if (trip) {
178 :			cholmod_free_sparse(&chx, &c);
179 :	bates	1141	Free(cht);
180 :	bates	930	} else {
181 :	bates	1141	Free(chx);
182 :	bates	930	}
183 :	bates	923	if (!tr) cholmod_free_sparse(&chxt, &c);
184 :	bates	1366	/* create dimnames */
185 :			SET_VECTOR_ELT(dn, 0,
186 :			duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym),
187 :			(tr) ? 1 : 0)));
188 :			SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0)));
189 :			UNPROTECT(1);
190 :			return chm_sparse_to_SEXP(chcp, 1, 0, "", dn);
191 :	bates	922	}
192 :	bates	923
193 :	bates	1218	SEXP Csparse_horzcat(SEXP x, SEXP y)
194 :			{
195 :			cholmod_sparse *chx = as_cholmod_sparse(x),
196 :			*chy = as_cholmod_sparse(y), *ans;
197 :	maechler	1375
198 :	bates	1218	ans = cholmod_horzcat(chx, chy, 1, &c);
199 :			Free(chx); Free(chy);
200 :	bates	1366	/* FIXME: currently drops dimnames */
201 :			return chm_sparse_to_SEXP(ans, 1, 0, "", R_NilValue);
202 :	bates	1218	}
203 :
204 :			SEXP Csparse_vertcat(SEXP x, SEXP y)
205 :			{
206 :			cholmod_sparse *chx = as_cholmod_sparse(x),
207 :			*chy = as_cholmod_sparse(y), *ans;
208 :	maechler	1375
209 :	bates	1218	ans = cholmod_vertcat(chx, chy, 1, &c);
210 :			Free(chx); Free(chy);
211 :	bates	1366	/* FIXME: currently drops dimnames */
212 :			return chm_sparse_to_SEXP(ans, 1, 0, "", R_NilValue);
213 :	bates	1218	}
214 :	bates	1265
215 :			SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2)
216 :			{
217 :			cholmod_sparse *chx = as_cholmod_sparse(x), *ans;
218 :
219 :			ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c);
220 :			Free(chx);
221 :	bates	1366	return chm_sparse_to_SEXP(ans, 1, 0, "", R_NilValue);
222 :	bates	1265	}
223 :	bates	1366
224 :			SEXP Csparse_diagU2N(SEXP x)
225 :			{
226 :			cholmod_sparse *chx = as_cholmod_sparse(x);
227 :			cholmod_sparse *eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c);
228 :			double one[] = {1, 0};
229 :			cholmod_sparse *ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c);
230 :			int uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ?
231 :			-1 : 1;
232 :
233 :			Free(chx); cholmod_free_sparse(&eye, &c);
234 :			return chm_sparse_to_SEXP(ans, 1, uploT, "N",
235 :			duplicate(GET_SLOT(x, Matrix_DimNamesSym)));
236 :			}
237 :
238 :			SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j)
239 :			{
240 :			cholmod_sparse *chx = as_cholmod_sparse(x);
241 :			int rsize = (isNull(i)) ? -1 : LENGTH(i),
242 :			csize = (isNull(j)) ? -1 : LENGTH(j);
243 :
244 :			if (rsize >= 0 && !isInteger(i))
245 :			error(_("Index i must be NULL or integer"));
246 :			if (csize >= 0 && !isInteger(j))
247 :			error(_("Index j must be NULL or integer"));
248 :			return chm_sparse_to_SEXP(cholmod_submatrix(chx, INTEGER(i), rsize,
249 :	maechler	1375	INTEGER(j), csize,
250 :	bates	1366	TRUE, TRUE, &c),
251 :			1, 0, "", R_NilValue);
252 :			}

---

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