Annotation of /pkg/Matrix/src/dgCMatrix.c

Parent Directory | Revision Log

---

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

1 :	mmaechler	3186	#ifdef __GLIBC__
2 :	mmaechler	3182	// to get strdup declared in glibc (when strict -std=c11 or -stdc99):
3 :			#define _POSIX_C_SOURCE 200809L
4 :	mmaechler	3186	#endif
5 :
6 :	mmaechler	3182	#include <string.h>
7 :
8 :	bates	478	#include "dgCMatrix.h"
9 :	bates	10
10 :	maechler	1747	/* for Csparse_transpose() : */
11 :			#include "Csparse.h"
12 :	bates	922	#include "chm_common.h"
13 :	mmaechler	2327	/* -> Mutils.h / SPQR ... */
14 :	bates	922
15 :	maechler	956	/* FIXME -- we "forget" about dimnames almost everywhere : */
16 :
17 :	maechler	1661	/* for dgCMatrix _and_ lgCMatrix and others (but *not* ngC...) : */
18 :			SEXP xCMatrix_validate(SEXP x)
19 :	bates	10	{
20 :	maechler	1660	/* Almost everything now in Csparse_validate ( ./Csparse.c )
21 :			* *but* the checking of the 'x' slot : */
22 :	mmaechler	3268	if (xlength(GET_SLOT(x, Matrix_iSym)) !=
23 :			xlength(GET_SLOT(x, Matrix_xSym)))
24 :	maechler	1660	return mkString(_("lengths of slots 'i' and 'x' must match"));
25 :	maechler	534
26 :	bates	10	return ScalarLogical(1);
27 :			}
28 :	maechler	534
29 :	maechler	1968	/* for dgRMatrix _and_ lgRMatrix and others (but *not* ngC...) : */
30 :			SEXP xRMatrix_validate(SEXP x)
31 :			{
32 :			/* Almost everything now in Rsparse_validate ( ./Csparse.c )
33 :			* *but* the checking of the 'x' slot : */
34 :	mmaechler	3268	if (xlength(GET_SLOT(x, Matrix_jSym)) !=
35 :			xlength(GET_SLOT(x, Matrix_xSym)))
36 :	maechler	1968	return mkString(_("lengths of slots 'j' and 'x' must match"));
37 :
38 :			return ScalarLogical(1);
39 :			}
40 :
41 :	maechler	1760	/* This and the following R_to_CMatrix() lead to memory-not-mapped seg.faults
42 :			* only with {32bit + R-devel + enable-R-shlib} -- no idea why */
43 :	maechler	1747	SEXP compressed_to_TMatrix(SEXP x, SEXP colP)
44 :	bates	10	{
45 :	maechler	890	int col = asLogical(colP); /* 1 if "C"olumn compressed; 0 if "R"ow */
46 :	maechler	1760	/* however, for Csparse, we now effectively use the cholmod-based
47 :			* Csparse_to_Tsparse() in ./Csparse.c ; maybe should simply write
48 :			* an as_cholmod_Rsparse() function and then do "as there" ...*/
49 :	mmaechler	3213	SEXP indSym = col ? Matrix_iSym : Matrix_jSym, ans,
50 :			indP = PROTECT(GET_SLOT(x, indSym)),
51 :			pP = PROTECT(GET_SLOT(x, Matrix_pSym));
52 :	bates	679	int npt = length(pP) - 1;
53 :	bates	1867	char *ncl = strdup(class_P(x));
54 :	mmaechler	2713	static const char *valid[] = { MATRIX_VALID_Csparse, MATRIX_VALID_Rsparse, ""};
55 :	mmaechler	3147	int ctype = R_check_class_etc(x, valid);
56 :	maechler	1747
57 :			if (ctype < 0)
58 :	bates	1867	error(_("invalid class(x) '%s' in compressed_to_TMatrix(x)"), ncl);
59 :	maechler	1747
60 :	maechler	1766	/* replace 'C' or 'R' with 'T' :*/
61 :	bates	1763	ncl[2] = 'T';
62 :	mmaechler	3270	ans = PROTECT(NEW_OBJECT_OF_CLASS(ncl));
63 :	maechler	1747
64 :	maechler	1968	slot_dup(ans, x, Matrix_DimSym);
65 :	maechler	1747	if((ctype / 3) % 4 != 2) /* not n..Matrix */
66 :	maechler	1968	slot_dup(ans, x, Matrix_xSym);
67 :	maechler	1747	if(ctype % 3) { /* s(ymmetric) or t(riangular) : */
68 :	maechler	1968	slot_dup(ans, x, Matrix_uploSym);
69 :	maechler	1747	if(ctype % 3 == 2) /* t(riangular) : */
70 :	maechler	1968	slot_dup(ans, x, Matrix_diagSym);
71 :	maechler	1747	}
72 :			SET_DimNames(ans, x);
73 :	mmaechler	3055	// possibly asymmetric for symmetricMatrix is ok
74 :	bates	677	SET_SLOT(ans, indSym, duplicate(indP));
75 :	bates	679	expand_cmprPt(npt, INTEGER(pP),
76 :			INTEGER(ALLOC_SLOT(ans, col ? Matrix_jSym : Matrix_iSym,
77 :			INTSXP, length(indP))));
78 :	bates	1763	free(ncl);
79 :	mmaechler	3213	UNPROTECT(3);
80 :	bates	10	return ans;
81 :			}
82 :
83 :	maechler	1747	SEXP R_to_CMatrix(SEXP x)
84 :			{
85 :			SEXP ans, tri = PROTECT(allocVector(LGLSXP, 1));
86 :	bates	1867	char *ncl = strdup(class_P(x));
87 :	mmaechler	2713	static const char *valid[] = { MATRIX_VALID_Rsparse, ""};
88 :	mmaechler	3147	int ctype = R_check_class_etc(x, valid);
89 :	maechler	1747	int *x_dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), *a_dims;
90 :	mmaechler	2455	PROTECT_INDEX ipx;
91 :	maechler	1747
92 :			if (ctype < 0)
93 :	bates	1867	error(_("invalid class(x) '%s' in R_to_CMatrix(x)"), ncl);
94 :	maechler	1747
95 :	maechler	1766	/* replace 'R' with 'C' : */
96 :			ncl[2] = 'C';
97 :	mmaechler	3270	PROTECT_WITH_INDEX(ans = NEW_OBJECT_OF_CLASS(ncl), &ipx);
98 :	maechler	1766
99 :	maechler	1747	a_dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
100 :			/* reversed dim() since we will transpose: */
101 :			a_dims[0] = x_dims[1];
102 :			a_dims[1] = x_dims[0];
103 :
104 :			/* triangular: */ LOGICAL(tri)[0] = 0;
105 :			if((ctype / 3) != 2) /* not n..Matrix */
106 :	maechler	1968	slot_dup(ans, x, Matrix_xSym);
107 :	maechler	1747	if(ctype % 3) { /* s(ymmetric) or t(riangular) : */
108 :	maechler	1968	SET_SLOT(ans, Matrix_uploSym,
109 :			mkString((*uplo_P(x) == 'U') ? "L" : "U"));
110 :	maechler	1747	if(ctype % 3 == 2) { /* t(riangular) : */
111 :			LOGICAL(tri)[0] = 1;
112 :	maechler	1968	slot_dup(ans, x, Matrix_diagSym);
113 :	maechler	1747	}
114 :			}
115 :			SET_SLOT(ans, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_jSym)));
116 :	maechler	1968	slot_dup(ans, x, Matrix_pSym);
117 :	mmaechler	2455	REPROTECT(ans = Csparse_transpose(ans, tri), ipx);
118 :	maechler	1747	SET_DimNames(ans, x);
119 :	mmaechler	3055	// possibly asymmetric for symmetricMatrix is ok
120 :	maechler	1771	free(ncl);
121 :	maechler	1747	UNPROTECT(2);
122 :			return ans;
123 :			}
124 :
125 :	mmaechler	2814	/** Return a 2 column matrix '' cbind(i, j) '' of 0-origin index vectors (i,j)
126 :			* which entirely correspond to the (i,j) slots of
127 :			* as(x, "TsparseMatrix") :
128 :			*/
129 :	maechler	919	SEXP compressed_non_0_ij(SEXP x, SEXP colP)
130 :			{
131 :			int col = asLogical(colP); /* 1 if "C"olumn compressed; 0 if "R"ow */
132 :			SEXP ans, indSym = col ? Matrix_iSym : Matrix_jSym;
133 :	mmaechler	3213	SEXP indP = PROTECT(GET_SLOT(x, indSym)),
134 :			pP = PROTECT(GET_SLOT(x, Matrix_pSym));
135 :	mmaechler	2236	int i, *ij;
136 :	dmbates	2759	int nouter = INTEGER(GET_SLOT(x, Matrix_DimSym))[col ? 1 : 0],
137 :			n_el = INTEGER(pP)[nouter]; /* is only == length(indP), if the
138 :			inner slot is not over-allocated */
139 :	maechler	919
140 :			ij = INTEGER(ans = PROTECT(allocMatrix(INTSXP, n_el, 2)));
141 :			/* expand the compressed margin to 'i' or 'j' : */
142 :	dmbates	2759	expand_cmprPt(nouter, INTEGER(pP), &ij[col ? n_el : 0]);
143 :	maechler	919	/* and copy the other one: */
144 :			if (col)
145 :			for(i = 0; i < n_el; i++)
146 :			ij[i] = INTEGER(indP)[i];
147 :			else /* row compressed */
148 :			for(i = 0; i < n_el; i++)
149 :			ij[i + n_el] = INTEGER(indP)[i];
150 :
151 :	mmaechler	3213	UNPROTECT(3);
152 :	maechler	919	return ans;
153 :			}
154 :
155 :	dmbates	2401	#if 0 /* unused */
156 :	bates	1380	SEXP dgCMatrix_lusol(SEXP x, SEXP y)
157 :	bates	332	{
158 :	mmaechler	2237	SEXP ycp = PROTECT((TYPEOF(y) == REALSXP) ?
159 :			duplicate(y) : coerceVector(y, REALSXP));
160 :	mmaechler	2227	CSP xc = AS_CSP__(x);
161 :	maechler	1960	R_CheckStack();
162 :	bates	310
163 :	bates	1380	if (xc->m != xc->n || xc->m <= 0)
164 :			error(_("dgCMatrix_lusol requires a square, non-empty matrix"));
165 :	mmaechler	2237	if (LENGTH(ycp) != xc->m)
166 :	bates	1380	error(_("Dimensions of system to be solved are inconsistent"));
167 :	bates	1383	if (!cs_lusol(/*order*/ 1, xc, REAL(ycp), /*tol*/ 1e-7))
168 :	bates	1380	error(_("cs_lusol failed"));
169 :	maechler	1960
170 :	bates	1380	UNPROTECT(1);
171 :			return ycp;
172 :			}
173 :	dmbates	2401	#endif
174 :	bates	332
175 :	mmaechler	2978	// called from package MatrixModels's R code
176 :	mmaechler	2239	SEXP dgCMatrix_qrsol(SEXP x, SEXP y, SEXP ord)
177 :	bates	1380	{
178 :	mmaechler	2529	/* FIXME: extend this to work in multivariate case, i.e. y a matrix with > 1 column ! */
179 :	mmaechler	2237	SEXP ycp = PROTECT((TYPEOF(y) == REALSXP) ?
180 :			duplicate(y) : coerceVector(y, REALSXP));
181 :			CSP xc = AS_CSP(x); /* <--> x may be dgC* or dtC* */
182 :	mmaechler	3033	int order = asInteger(ord);
183 :	mmaechler	2529	#ifdef _not_yet_do_FIXME__
184 :			const char *nms[] = {"L", "coef", "Xty", "resid", ""};
185 :	mmaechler	2645	SEXP ans = PROTECT(Rf_mkNamed(VECSXP, nms));
186 :	mmaechler	2529	#endif
187 :	maechler	1960	R_CheckStack();
188 :	bates	1380
189 :	mmaechler	2239	if (order < 0 || order > 3)
190 :	mmaechler	2529	error(_("dgCMatrix_qrsol(., order) needs order in {0,..,3}"));
191 :	mmaechler	2239	/* --> cs_amd() --- order 0: natural, 1: Chol, 2: LU, 3: QR */
192 :			if (LENGTH(ycp) != xc->m)
193 :			error(_("Dimensions of system to be solved are inconsistent"));
194 :			/* FIXME? Note that qr_sol() would allow *under-determined systems;
195 :	mmaechler	2529	* In general, we'd need LENGTH(ycp) = max(n,m)
196 :			* FIXME also: multivariate y (see above)
197 :	mmaechler	2239	*/
198 :	bates	1380	if (xc->m < xc->n || xc->n <= 0)
199 :	mmaechler	2237	error(_("dgCMatrix_qrsol(<%d x %d>-matrix) requires a 'tall' rectangular matrix"),
200 :			xc->m, xc->n);
201 :	mmaechler	2239
202 :			/* cs_qrsol(): Tim Davis (2006) .. "8.2 Using a QR factorization", p.136f , calling
203 :			* ------- cs_sqr(order, ..), see p.76 */
204 :	mmaechler	2529	/* MM: FIXME: write our *OWN* version of - the first case (m >= n) - of cs_qrsol()
205 :			* --------- which will (1) work with a *multivariate* y
206 :			* (2) compute coefficients properly, not overwriting RHS
207 :			*/
208 :	mmaechler	2239	if (!cs_qrsol(order, xc, REAL(ycp)))
209 :			/* return value really is 0 or 1 - no more info there */
210 :	mmaechler	2529	error(_("cs_qrsol() failed inside dgCMatrix_qrsol()"));
211 :	maechler	1960
212 :	mmaechler	2239	/* Solution is only in the first part of ycp -- cut its length back to n : */
213 :	mmaechler	3161	ycp = lengthgets(ycp, (R_xlen_t) xc->n);
214 :	mmaechler	2879
215 :	bates	332	UNPROTECT(1);
216 :	bates	1380	return ycp;
217 :	bates	332	}
218 :	bates	1380
219 :	mmaechler	2854	// Modified version of Tim Davis's cs_qr_mex.c file for MATLAB (in CSparse)
220 :			// Usage: [V,beta,p,R,q] = cs_qr(A) ;
221 :	mmaechler	3115	SEXP dgCMatrix_QR(SEXP Ap, SEXP order, SEXP keep_dimnames)
222 :	bates	1380	{
223 :	mmaechler	2227	CSP A = AS_CSP__(Ap), D;
224 :	mmaechler	2854	int io = INTEGER(order)[0];
225 :	mmaechler	2933	Rboolean verbose = (io < 0);// verbose=TRUE, encoded with negative 'order'
226 :	mmaechler	3115	int m0 = A->m, m = m0, n = A->n, ord = asLogical(order) ? 3 : 0, *p;
227 :	maechler	1960	R_CheckStack();
228 :	bates	1380
229 :	mmaechler	2387	if (m < n) error(_("A must have #{rows} >= #{columns}")) ;
230 :	mmaechler	3273	SEXP ans = PROTECT(NEW_OBJECT_OF_CLASS("sparseQR"));
231 :	mmaechler	2854	int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
232 :	maechler	2106	dims[0] = m; dims[1] = n;
233 :	mmaechler	2854	css *S = cs_sqr(ord, A, 1); /* symbolic QR ordering & analysis*/
234 :	mmaechler	2387	if (!S) error(_("cs_sqr failed"));
235 :	mmaechler	3115	int keep_dimnms = asLogical(keep_dimnames);
236 :			if(keep_dimnms == NA_LOGICAL) { keep_dimnms = TRUE;
237 :			warning(_("dgcMatrix_QR(*, keep_dimnames = NA): NA taken as TRUE"));
238 :			}
239 :	mmaechler	2854	if(verbose && S->m2 > m) // in ./cs.h , m2 := # of rows for QR, after adding fictitious rows
240 :			Rprintf("Symbolic QR(): Matrix structurally rank deficient (m2-m = %d)\n",
241 :			S->m2 - m);
242 :			csn *N = cs_qr(A, S); /* numeric QR factorization */
243 :	mmaechler	2387	if (!N) error(_("cs_qr failed")) ;
244 :	bates	1381	cs_dropzeros(N->L); /* drop zeros from V and sort */
245 :	mmaechler	2854	D = cs_transpose(N->L, 1); cs_spfree(N->L);
246 :			N->L = cs_transpose(D, 1); cs_spfree(D);
247 :	bates	1381	cs_dropzeros(N->U); /* drop zeros from R and sort */
248 :	mmaechler	2854	D = cs_transpose(N->U, 1); cs_spfree(N->U) ;
249 :			N->U = cs_transpose(D, 1); cs_spfree(D);
250 :	bates	1381	m = N->L->m; /* m may be larger now */
251 :	mmaechler	2854	// MM: m := S->m2 also counting the ficticious rows (Tim Davis, p.72, 74f)
252 :	bates	1383	p = cs_pinv(S->pinv, m); /* p = pinv' */
253 :	mmaechler	3115	SEXP dn = R_NilValue; Rboolean do_dn = FALSE;
254 :			if(keep_dimnms) {
255 :			dn = GET_SLOT(Ap, Matrix_DimNamesSym);
256 :			do_dn = !isNull(VECTOR_ELT(dn, 0)) && m == m0;
257 :			// FIXME? also deal with case m > m0 ?
258 :			if(do_dn) { // keep rownames
259 :			dn = PROTECT(duplicate(dn));
260 :			SET_VECTOR_ELT(dn, 1, R_NilValue);
261 :			} else dn = R_NilValue;
262 :			}
263 :			SET_SLOT(ans, Matrix_VSym, Matrix_cs_to_SEXP(N->L, "dgCMatrix", 0, dn)); // "V"
264 :			Memcpy(REAL(ALLOC_SLOT(ans, Matrix_betaSym, REALSXP, n)), N->B, n);
265 :			Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, m)), p, m);
266 :			if(do_dn) {
267 :			UNPROTECT(1); // dn
268 :			dn = R_NilValue; do_dn = FALSE;
269 :			}
270 :			if (ord) {
271 :			Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"), INTSXP, n)), S->q, n);
272 :			if(keep_dimnms) {
273 :			dn = GET_SLOT(Ap, Matrix_DimNamesSym);
274 :			do_dn = !isNull(VECTOR_ELT(dn, 1));
275 :			if(do_dn) {
276 :			dn = PROTECT(duplicate(dn));
277 :			// permute colnames by S->q : cn <- cn[ S->q ] :
278 :			SEXP cns = PROTECT(duplicate(VECTOR_ELT(dn, 1)));
279 :			for(int j=0; j < n; j++)
280 :			SET_STRING_ELT(VECTOR_ELT(dn, 1), j, STRING_ELT(cns, S->q[j]));
281 :			UNPROTECT(1);
282 :			SET_VECTOR_ELT(dn, 0, R_NilValue);
283 :			} else dn = R_NilValue;
284 :			}
285 :			} else
286 :	bates	1397	ALLOC_SLOT(ans, install("q"), INTSXP, 0);
287 :	mmaechler	3115	SET_SLOT(ans, install("R"), Matrix_cs_to_SEXP(N->U, "dgCMatrix", 0, dn));
288 :			if(do_dn) UNPROTECT(1); // dn
289 :	bates	1380	cs_nfree(N);
290 :			cs_sfree(S);
291 :			cs_free(p);
292 :			UNPROTECT(1);
293 :			return ans;
294 :			}
295 :
296 :	mmaechler	2336	#ifdef Matrix_with_SPQR
297 :	dmbates	2276	/**
298 :			* Return a SuiteSparse QR factorization of the sparse matrix A
299 :			*
300 :	mmaechler	2322	* @param Ap (pointer to) a [m x n] dgCMatrix
301 :			* @param ordering integer SEXP specifying the ordering strategy to be used
302 :			* see SPQR/Include/SuiteSparseQR_definitions.h
303 :			* @param econ integer SEXP ("economy"): number of rows of R and columns of Q
304 :			* to return. The default is m. Using n gives the standard economy form.
305 :			* A value less than the estimated rank r is set to r, so econ=0 gives the
306 :			* "rank-sized" factorization, where nrow(R)==nnz(diag(R))==r.
307 :	mmaechler	2327	* @param tol double SEXP: if tol <= -2 use SPQR's default,
308 :			* if -2 < tol < 0, then no tol is used; otherwise,
309 :			* tol > 0, use as tolerance: columns with 2-norm <= tol treated as 0
310 :	dmbates	2276	*
311 :	mmaechler	2327	*
312 :			* @return SEXP "SPQR" object with slots (Q, R, p, rank, Dim):
313 :	mmaechler	2322	* Q: dgCMatrix; R: dgCMatrix [subject to change to dtCMatrix FIXME ?]
314 :	mmaechler	2327	* p: integer: 0-based permutation (or length 0 <=> identity);
315 :			* rank: integer, the "revealed" rank Dim: integer, original matrix dim.
316 :	dmbates	2276	*/
317 :	mmaechler	2322	SEXP dgCMatrix_SPQR(SEXP Ap, SEXP ordering, SEXP econ, SEXP tol)
318 :	dmbates	2276	{
319 :	mmaechler	2327	/* SEXP ans = PROTECT(allocVector(VECSXP, 4)); */
320 :	mmaechler	3273	SEXP ans = PROTECT(NEW_OBJECT_OF_CLASS("SPQR"));
321 :	mmaechler	2327
322 :	dmbates	2299	CHM_SP A = AS_CHM_SP(Ap), Q, R;
323 :	mmaechler	3001	SuiteSparse_long *E, rank;/* not always = int FIXME (Windows_64 ?) */
324 :	mmaechler	2322
325 :			if ((rank = SuiteSparseQR_C_QR(asInteger(ordering),
326 :			asReal(tol),/* originally had SPQR_DEFAULT_TOL */
327 :	mmaechler	3001	(SuiteSparse_long)asInteger(econ),/* originally had 0 */
328 :	mmaechler	2661	A, &Q, &R, &E, &cl)) == -1)
329 :	dmbates	2276	error(_("SuiteSparseQR_C_QR returned an error code"));
330 :	mmaechler	2327
331 :	mmaechler	3011	slot_dup(ans, Ap, Matrix_DimSym);
332 :	mmaechler	2327	/* SET_VECTOR_ELT(ans, 0, */
333 :			/* chm_sparse_to_SEXP(Q, 0, 0, 0, "", R_NilValue)); */
334 :			SET_SLOT(ans, install("Q"),
335 :			chm_sparse_to_SEXP(Q, 0, 0, 0, "", R_NilValue));
336 :
337 :	mmaechler	2322	/* Also gives a dgCMatrix (not a dtC* *triangular*) :
338 :			* may make sense if to be used in the "spqr_solve" routines .. ?? */
339 :	mmaechler	2327	/* SET_VECTOR_ELT(ans, 1, */
340 :			/* chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue)); */
341 :			SET_SLOT(ans, install("R"),
342 :			chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue));
343 :	mmaechler	2661	cholmod_free_sparse(&Al, &cl);
344 :			cholmod_free_sparse(&R, &cl);
345 :			cholmod_free_sparse(&Q, &cl);
346 :	dmbates	2276	if (E) {
347 :	mmaechler	2661	int *Er;
348 :			SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, A->ncol));
349 :			Er = INTEGER(VECTOR_ELT(ans, 2));
350 :			for (int i = 0; i < A->ncol; i++) Er[i] = (int) E[i];
351 :	dmbates	2276	Free(E);
352 :	mmaechler	2661	} else SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, 0));
353 :			SET_VECTOR_ELT(ans, 3, ScalarInteger((int)rank));
354 :	dmbates	2276	UNPROTECT(1);
355 :			return ans;
356 :			}
357 :	mmaechler	2338	#endif
358 :			/* Matrix_with_SPQR */
359 :	dmbates	2276
360 :	bates	1463	/* Modified version of Tim Davis's cs_lu_mex.c file for MATLAB */
361 :	mmaechler	3115	void install_lu(SEXP Ap, int order, double tol, Rboolean err_sing, Rboolean keep_dimnms)
362 :	bates	1383	{
363 :	mmaechler	2713	// (order, tol) == (1, 1) by default, when called from R.
364 :	dmbates	2401	SEXP ans;
365 :	bates	1383	css *S;
366 :			csn *N;
367 :	dmbates	2401	int n, *p, *dims;
368 :			CSP A = AS_CSP__(Ap), D;
369 :	maechler	1960	R_CheckStack();
370 :	bates	1383
371 :	bates	1385	n = A->n;
372 :	bates	1383	if (A->m != n)
373 :	mmaechler	2387	error(_("LU decomposition applies only to square matrices"));
374 :	bates	1383	if (order) { /* not using natural order */
375 :			order = (tol == 1) ? 2 /* amd(S'*S) w/dense rows or I */
376 :			: 1; /* amd (A+A'), or natural */
377 :			}
378 :	mmaechler	2713	S = cs_sqr(order, A, /*qr = */ 0); /* symbolic ordering */
379 :	dmbates	2401	N = cs_lu(A, S, tol); /* numeric factorization */
380 :	mmaechler	2488	if (!N) {
381 :			if(err_sing)
382 :			error(_("cs_lu(A) failed: near-singular A (or out of memory)"));
383 :			else {
384 :			/* No warning: The useR should be careful :
385 :	mmaechler	2692	* Put NA into "LU" factor */
386 :	mmaechler	2488	set_factors(Ap, ScalarLogical(NA_LOGICAL), "LU");
387 :			return;
388 :			}
389 :			}
390 :	dmbates	2401	cs_dropzeros(N->L); /* drop zeros from L and sort it */
391 :			D = cs_transpose(N->L, 1);
392 :			cs_spfree(N->L);
393 :			N->L = cs_transpose(D, 1);
394 :			cs_spfree(D);
395 :			cs_dropzeros(N->U); /* drop zeros from U and sort it */
396 :			D = cs_transpose(N->U, 1);
397 :			cs_spfree(N->U);
398 :			N->U = cs_transpose(D, 1);
399 :			cs_spfree(D);
400 :			p = cs_pinv(N->pinv, n); /* p=pinv' */
401 :	mmaechler	3273	ans = PROTECT(NEW_OBJECT_OF_CLASS("sparseLU"));
402 :	mmaechler	2268	dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
403 :			dims[0] = n; dims[1] = n;
404 :	mmaechler	3115	SEXP dn; Rboolean do_dn = FALSE;
405 :			if(keep_dimnms) {
406 :			dn = GET_SLOT(Ap, Matrix_DimNamesSym);
407 :			do_dn = !isNull(VECTOR_ELT(dn, 0));
408 :			if(do_dn) {
409 :			dn = PROTECT(duplicate(dn));
410 :			// permute rownames by p : rn <- rn[ p ] :
411 :			SEXP rn = PROTECT(duplicate(VECTOR_ELT(dn, 0)));
412 :			for(int i=0; i < n; i++)
413 :			SET_STRING_ELT(VECTOR_ELT(dn, 0), i, STRING_ELT(rn, p[i]));
414 :			UNPROTECT(1); // rn
415 :			SET_VECTOR_ELT(dn, 1, R_NilValue); // colnames(.) := NULL
416 :			}
417 :			}
418 :	bates	1383	SET_SLOT(ans, install("L"),
419 :	mmaechler	3115	Matrix_cs_to_SEXP(N->L, "dtCMatrix", 0, do_dn ? dn : R_NilValue));
420 :
421 :			if(keep_dimnms) {
422 :			if(do_dn) {
423 :			UNPROTECT(1); // dn
424 :			dn = GET_SLOT(Ap, Matrix_DimNamesSym);
425 :			}
426 :			do_dn = !isNull(VECTOR_ELT(dn, 1));
427 :			if(do_dn) {
428 :			dn = PROTECT(duplicate(dn));
429 :			if(order) { // permute colnames by S->q : cn <- cn[ S->q ] :
430 :			SEXP cn = PROTECT(duplicate(VECTOR_ELT(dn, 1)));
431 :			for(int j=0; j < n; j++)
432 :			SET_STRING_ELT(VECTOR_ELT(dn, 1), j, STRING_ELT(cn, S->q[j]));
433 :			UNPROTECT(1); // cn
434 :			}
435 :			SET_VECTOR_ELT(dn, 0, R_NilValue); // rownames(.) := NULL
436 :			}
437 :			}
438 :	bates	1383	SET_SLOT(ans, install("U"),
439 :	mmaechler	3115	Matrix_cs_to_SEXP(N->U, "dtCMatrix", 0, do_dn ? dn : R_NilValue));
440 :			if(do_dn) UNPROTECT(1); // dn
441 :	mmaechler	2247	Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, /* "p" */
442 :	bates	1383	INTSXP, n)), p, n);
443 :			if (order)
444 :			Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"),
445 :			INTSXP, n)), S->q, n);
446 :			cs_nfree(N);
447 :			cs_sfree(S);
448 :			cs_free(p);
449 :			UNPROTECT(1);
450 :	dmbates	2401	set_factors(Ap, ans, "LU");
451 :	bates	1385	}
452 :
453 :	mmaechler	3115	SEXP dgCMatrix_LU(SEXP Ap, SEXP orderp, SEXP tolp, SEXP error_on_sing, SEXP keep_dimnames)
454 :	dmbates	2401	{
455 :			SEXP ans;
456 :	mmaechler	2488	Rboolean err_sing = asLogical(error_on_sing);
457 :	dmbates	2401	/* FIXME: dgCMatrix_LU should check ans for consistency in
458 :			* permutation type with the requested value - Should have two
459 :			* classes or two different names in the factors list for LU with
460 :			* permuted columns or not.
461 :			* OTOH, currently (order, tol) === (1, 1) always.
462 :	dmbates	2403	* It is true that length(LU@q) does flag the order argument.
463 :	dmbates	2401	*/
464 :			if (!isNull(ans = get_factors(Ap, "LU")))
465 :			return ans;
466 :	mmaechler	3115	int keep_dimnms = asLogical(keep_dimnames);
467 :			if(keep_dimnms == NA_LOGICAL) { keep_dimnms = TRUE;
468 :			warning(_("dgcMatrix_LU(*, keep_dimnames = NA): NA taken as TRUE"));
469 :			}
470 :			install_lu(Ap, asInteger(orderp), asReal(tolp), err_sing, keep_dimnms);
471 :	dmbates	2401	return get_factors(Ap, "LU");
472 :			}
473 :
474 :	mmaechler	2889	SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b, SEXP give_sparse)
475 :	mmaechler	3115	// FIXME: add 'keep_dimnames' as argument
476 :	bates	1385	{
477 :	mmaechler	2889	Rboolean sparse = asLogical(give_sparse);
478 :			if(sparse) {
479 :			// FIXME: implement this
480 :			error(_("dgCMatrix_matrix_solve(.., sparse=TRUE) not yet implemented"));
481 :
482 :			/* Idea: in the for(j = 0; j < nrhs ..) loop below, build the *sparse* result matrix
483 :			* ----- *column* wise -- which is perfect for dgCMatrix
484 :			* --> build (i,p,x) slots "increasingly" [well, allocate in batches ..]
485 :			*
486 :			* --> maybe first a protoype in R
487 :			*/
488 :
489 :			}
490 :	dmbates	2401	SEXP ans = PROTECT(dup_mMatrix_as_dgeMatrix(b)),
491 :			lu, qslot;
492 :			CSP L, U;
493 :			int *bdims = INTEGER(GET_SLOT(ans, Matrix_DimSym)), *p, *q;
494 :	bates	1385	int j, n = bdims[0], nrhs = bdims[1];
495 :	mmaechler	3105	double *x, *ax = REAL(GET_SLOT(ans, Matrix_xSym));
496 :			C_or_Alloca_TO(x, n, double);
497 :	bates	1385
498 :	dmbates	2401	if (isNull(lu = get_factors(Ap, "LU"))) {
499 :	mmaechler	3115	install_lu(Ap, /* order = */ 1, /* tol = */ 1.0, /* err_sing = */ TRUE,
500 :			/* keep_dimnames = */ TRUE);
501 :	dmbates	2401	lu = get_factors(Ap, "LU");
502 :			}
503 :			qslot = GET_SLOT(lu, install("q"));
504 :			L = AS_CSP__(GET_SLOT(lu, install("L")));
505 :			U = AS_CSP__(GET_SLOT(lu, install("U")));
506 :			R_CheckStack();
507 :	mmaechler	3011	if (U->n != n)
508 :			error(_("Dimensions of system to be solved are inconsistent"));
509 :			if(nrhs >= 1 && n >= 1) {
510 :			p = INTEGER(GET_SLOT(lu, Matrix_pSym));
511 :			q = LENGTH(qslot) ? INTEGER(qslot) : (int *) NULL;
512 :	dmbates	2401
513 :	mmaechler	3011	for (j = 0; j < nrhs; j++) {
514 :			cs_pvec(p, ax + j * n, x, n); /* x = b(p) */
515 :			cs_lsolve(L, x); /* x = L\x */
516 :			cs_usolve(U, x); /* x = U\x */
517 :	mmaechler	3023	if (q) /* r(q) = x , hence
518 :			r = Q' U{^-1} L{^-1} P b = A^{-1} b */
519 :	mmaechler	3011	cs_ipvec(q, x, ax + j * n, n);
520 :			else
521 :			Memcpy(ax + j * n, x, n);
522 :			}
523 :	bates	1385	}
524 :	mmaechler	3105	if(n >= SMALL_4_Alloca) Free(x);
525 :	dmbates	2401	UNPROTECT(1);
526 :	bates	1383	return ans;
527 :			}
528 :	bates	1385
529 :	mmaechler	2978	// called from package MatrixModels's R code:
530 :	bates	1895	SEXP dgCMatrix_cholsol(SEXP x, SEXP y)
531 :			{
532 :	mmaechler	2239	/* Solve Sparse Least Squares X %*% beta ~= y with dense RHS y,
533 :			* where X = t(x) i.e. we pass x = t(X) as argument,
534 :			* via "Cholesky(X'X)" .. well not really:
535 :			* cholmod_factorize("x", ..) finds L in X'X = L'L directly */
536 :	mmaechler	2237	CHM_SP cx = AS_CHM_SP(x);
537 :	mmaechler	2529	/* FIXME: extend this to work in multivariate case, i.e. y a matrix with > 1 column ! */
538 :	mmaechler	3276	SEXP y_ = PROTECT(coerceVector(y, REALSXP));
539 :			CHM_DN cy = AS_CHM_DN(y_), rhs, cAns, resid;
540 :	maechler	1960	CHM_FR L;
541 :	mmaechler	2529	int n = cx->ncol;/* #{obs.} {x = t(X) !} */
542 :			double one[] = {1,0}, zero[] = {0,0}, neg1[] = {-1,0};
543 :			const char *nms[] = {"L", "coef", "Xty", "resid", ""};
544 :	mmaechler	2645	SEXP ans = PROTECT(Rf_mkNamed(VECSXP, nms));
545 :	maechler	1960	R_CheckStack();
546 :	bates	1895
547 :	mmaechler	2529	if (n < cx->nrow || n <= 0)
548 :	bates	1895	error(_("dgCMatrix_cholsol requires a 'short, wide' rectangular matrix"));
549 :	mmaechler	2529	if (cy->nrow != n)
550 :	bates	1895	error(_("Dimensions of system to be solved are inconsistent"));
551 :	mmaechler	2661	rhs = cholmod_allocate_dense(cx->nrow, 1, cx->nrow, CHOLMOD_REAL, &c);
552 :	mmaechler	2529	/* cholmod_sdmult(A, transp, alpha, beta, X, Y, &c):
553 :			* Y := alpha*(A*X) + beta*Y or alpha*(A'*X) + beta*Y ;
554 :			* here: rhs := 1 * x %*% y + 0 = x %*% y = X'y */
555 :	mmaechler	2661	if (!(cholmod_sdmult(cx, 0 /* trans */, one, zero, cy, rhs, &c)))
556 :			error(_("cholmod_sdmult error (rhs)"));
557 :			L = cholmod_analyze(cx, &c);
558 :			if (!cholmod_factorize(cx, L, &c))
559 :			error(_("cholmod_factorize failed: status %d, minor %d from ncol %d"),
560 :	bates	1895	c.status, L->minor, L->n);
561 :			/* FIXME: Do this in stages so an "effects" vector can be calculated */
562 :	mmaechler	2661	if (!(cAns = cholmod_solve(CHOLMOD_A, L, rhs, &c)))
563 :			error(_("cholmod_solve (CHOLMOD_A) failed: status %d, minor %d from ncol %d"),
564 :	bates	1895	c.status, L->minor, L->n);
565 :	mmaechler	2529	/* L : */
566 :	bates	1895	SET_VECTOR_ELT(ans, 0, chm_factor_to_SEXP(L, 0));
567 :	mmaechler	2529	/* coef : */
568 :	bates	1895	SET_VECTOR_ELT(ans, 1, allocVector(REALSXP, cx->nrow));
569 :			Memcpy(REAL(VECTOR_ELT(ans, 1)), (double*)(cAns->x), cx->nrow);
570 :	mmaechler	2529	/* X'y : */
571 :	bates	1895	/* FIXME: Change this when the "effects" vector is available */
572 :			SET_VECTOR_ELT(ans, 2, allocVector(REALSXP, cx->nrow));
573 :	mmaechler	2527	Memcpy(REAL(VECTOR_ELT(ans, 2)), (double*)(rhs->x), cx->nrow);
574 :	mmaechler	2529	/* resid := y */
575 :	mmaechler	2661	resid = cholmod_copy_dense(cy, &c);
576 :	mmaechler	2529	/* cholmod_sdmult(A, transp, alp, bet, X, Y, &c):
577 :			* Y := alp*(A*X) + bet*Y or alp*(A'*X) + beta*Y ;
578 :			* here: resid := -1 * x' %*% coef + 1 * y = y - X %*% coef */
579 :	mmaechler	2661	if (!(cholmod_sdmult(cx, 1/* trans */, neg1, one, cAns, resid, &c)))
580 :			error(_("cholmod_sdmult error (resid)"));
581 :	mmaechler	2529	/* FIXME: for multivariate case, i.e. resid *matrix* with > 1 column ! */
582 :			SET_VECTOR_ELT(ans, 3, allocVector(REALSXP, n));
583 :			Memcpy(REAL(VECTOR_ELT(ans, 3)), (double*)(resid->x), n);
584 :	bates	1895
585 :	mmaechler	2661	cholmod_free_factor(&L, &c);
586 :			cholmod_free_dense(&rhs, &c);
587 :			cholmod_free_dense(&cAns, &c);
588 :	mmaechler	3276	UNPROTECT(2);
589 :	bates	1895	return ans;
590 :			}
591 :	bates	1901
592 :
593 :	maechler	2015	/* Define all of
594 :	maechler	1911	* dgCMatrix_colSums(....)
595 :			* igCMatrix_colSums(....)
596 :	maechler	2015	* lgCMatrix_colSums_d(....)
597 :			* lgCMatrix_colSums_i(....)
598 :			* ngCMatrix_colSums_d(....)
599 :			* ngCMatrix_colSums_i(....)
600 :	maechler	1911	*/
601 :	maechler	1920	#define _dgC_
602 :			#include "t_gCMatrix_colSums.c"
603 :	bates	1901
604 :	maechler	1920	#define _igC_
605 :			#include "t_gCMatrix_colSums.c"
606 :	bates	1901
607 :	maechler	1920	#define _lgC_
608 :			#include "t_gCMatrix_colSums.c"
609 :	maechler	1911
610 :	maechler	1920	#define _ngC_
611 :			#include "t_gCMatrix_colSums.c"
612 :	maechler	1911
613 :	maechler	1920	#define _lgC_mn
614 :			#include "t_gCMatrix_colSums.c"
615 :	maechler	1911
616 :	maechler	1920	#define _ngC_mn
617 :			#include "t_gCMatrix_colSums.c"
618 :	maechler	1911
619 :
620 :	maechler	1920	SEXP lgCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means)
621 :			{
622 :			if(asLogical(means)) /* ==> result will be "double" / "dsparseVector" */
623 :			return lgCMatrix_colSums_d(x, NArm, spRes, trans, means);
624 :			else
625 :			return lgCMatrix_colSums_i(x, NArm, spRes, trans, means);
626 :			}
627 :	maechler	1911
628 :	maechler	1920	SEXP ngCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means)
629 :			{
630 :			if(asLogical(means)) /* ==> result will be "double" / "dsparseVector" */
631 :			return ngCMatrix_colSums_d(x, NArm, spRes, trans, means);
632 :			else
633 :			return ngCMatrix_colSums_i(x, NArm, spRes, trans, means);
634 :			}

---

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