The R Project SVN – R-packages – Diff – /branches/Matrix-APIchange/src/dgCMatrix.c



}
 
SEXP dgCMatrix_lusol(SEXP x, SEXP y)
{
    SEXP ycp = PROTECT(duplicate(y));
    CSP xc = AS_CSP(x);
 
    if (xc->m != xc->n || xc->m <= 0)
	error(_("dgCMatrix_lusol requires a square, non-empty matrix"));
    if (!isReal(ycp) || LENGTH(ycp) != xc->m)
	error(_("Dimensions of system to be solved are inconsistent"));
    if (!cs_lusol(/*order*/ 1, xc, REAL(ycp), /*tol*/ 1e-7))
	error(_("cs_lusol failed"));
 
    UNPROTECT(1);
    return ycp;
}
 
SEXP dgCMatrix_qrsol(SEXP x, SEXP y)
{
    SEXP ycp = PROTECT(duplicate(y));
    CSP xc = AS_CSP(x);
    R_CheckStack();
 
    if (xc->m < xc->n || xc->n <= 0)
	error(_("dgCMatrix_qrsol requires a 'tall' rectangular matrix"));
    if (!isReal(ycp) || LENGTH(ycp) != xc->m)
	error(_("Dimensions of system to be solved are inconsistent"));
    if (!cs_qrsol(/*order*/ 1, xc, REAL(ycp)))
	error(_("cs_qrsol failed"));
 
    UNPROTECT(1);
    return ycp;
}
 
/* Modified version of Tim Davis's cs_qr_mex.c file for MATLAB */
SEXP dgCMatrix_QR(SEXP Ap, SEXP order)
{
    SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseQR")));
    CSP A = AS_CSP(Ap), D;
    css *S;
    csn *N;
    int m = A->m, n = A->n, ord = asLogical(order) ? 3 : 0, *p;
    R_CheckStack();
 
    if (m < n) error("A must have # rows >= # columns") ;
    S = cs_sqr(ord, A, 1);	/* symbolic QR ordering & analysis*/
    if (!S) error("cs_sqr failed");
    N = cs_qr(A, S);		/* numeric QR factorization */

 
/* Modified version of Tim Davis's cs_lu_mex.c file for MATLAB */
SEXP dgCMatrix_LU(SEXP Ap, SEXP orderp, SEXP tolp)
{
    SEXP ans = get_factors(Ap, "LU");
    CSP A = AS_CSP(Ap), D;
    css *S;
    csn *N;
    int n, order = asInteger(orderp), *p;
    double tol = asReal(tolp);
    R_CheckStack();
 
    /* FIXME: dgCMatrix_LU should check ans for consistency in
     * permutation type with the requested value - Should have two
     * classes or two different names in the factors list for LU with
     * permuted columns or not. */
 
    if (ans != R_NilValue) return ans;
    ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseLU")));
 
    n = A->n;
    if (A->m != n)
	error("LU decomposition applies only to square matrices");
    if (order) {		/* not using natural order */
	order = (tol == 1) ? 2	/* amd(S'*S) w/dense rows or I */

	Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"),
				  INTSXP, n)), S->q, n);
    cs_nfree(N);
    cs_sfree(S);
    cs_free(p);
 
    UNPROTECT(1);
    return set_factors(Ap, ans, "LU");
}
 
SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b)
{
    /* b is dense or NULL [ <--> solve(A) */
    SEXP lu = dgCMatrix_LU(Ap, ScalarLogical(1), ScalarReal(1));
    SEXP qslot = GET_SLOT(lu, install("q"));
    CSP L = AS_CSP(GET_SLOT(lu, install("L"))),
	U = AS_CSP(GET_SLOT(lu, install("U")));
    SEXP ans = PROTECT( !isNull(b) ? dup_mMatrix_as_dgeMatrix(b)
			: new_dgeMatrix(U->n, U->n));
    int *bdims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
    int j, n = bdims[0], nrhs = bdims[1];
    int *p = INTEGER(GET_SLOT(lu, Matrix_pSym)),
	*q = LENGTH(qslot) ? INTEGER(qslot) : (int *) NULL;
    double *ax = REAL(GET_SLOT(ans, Matrix_xSym)),
	*x = Calloc(n, double);
    R_CheckStack();
 
    if (U->n != n || nrhs < 1 || n < 1)
	error(_("Dimensions of system to be solved are inconsistent"));
    for (j = 0; j < nrhs; j++) {
	if(!isNull(b))

	if (q)			       /* b(q) = x */
	    cs_ipvec(q, x, ax + j * n, n);
	else
	    Memcpy(ax + j * n, x, n);
    }
    Free(x);
    UNPROTECT(1);
    return ans;
}
 
SEXP dgCMatrix_cholsol(SEXP x, SEXP y)

Rev 4530	Rev 4543
Line 137...	Line 137...
137	`}`	137	`}`
138		138
139	`SEXP dgCMatrix_lusol(SEXP x, SEXP y)`	139	`SEXP dgCMatrix_lusol(SEXP x, SEXP y)`
140	`{`	140	`{`
141	`SEXP ycp = PROTECT(duplicate(y));`	141	`SEXP ycp = PROTECT(duplicate(y));`
142	`cs *xc = Matrix_as_cs(x);`	142	`CSP xc = AS_CSP(x);`
143		143
144	`if (xc->m != xc->n \|\| xc->m <= 0)`	144	`if (xc->m != xc->n \|\| xc->m <= 0)`
145	`error(_("dgCMatrix_lusol requires a square, non-empty matrix"));`	145	`error(_("dgCMatrix_lusol requires a square, non-empty matrix"));`
146	`if (!isReal(ycp) \|\| LENGTH(ycp) != xc->m)`	146	`if (!isReal(ycp) \|\| LENGTH(ycp) != xc->m)`
147	`error(_("Dimensions of system to be solved are inconsistent"));`	147	`error(_("Dimensions of system to be solved are inconsistent"));`
148	`if (!cs_lusol(/order/ 1, xc, REAL(ycp), /tol/ 1e-7))`	148	`if (!cs_lusol(/order/ 1, xc, REAL(ycp), /tol/ 1e-7))`
149	`error(_("cs_lusol failed"));`	149	`error(_("cs_lusol failed"));`
150	`Free(xc);`	150
151	`UNPROTECT(1);`	151	`UNPROTECT(1);`
152	`return ycp;`	152	`return ycp;`
153	`}`	153	`}`
154		154
155	`SEXP dgCMatrix_qrsol(SEXP x, SEXP y)`	155	`SEXP dgCMatrix_qrsol(SEXP x, SEXP y)`
156	`{`	156	`{`
157	`SEXP ycp = PROTECT(duplicate(y));`	157	`SEXP ycp = PROTECT(duplicate(y));`
158	`cs *xc = Matrix_as_cs(x);`	158	`CSP xc = AS_CSP(x);`
-		159	`R_CheckStack();`
159		160
160	`if (xc->m < xc->n \|\| xc->n <= 0)`	161	`if (xc->m < xc->n \|\| xc->n <= 0)`
161	`error(_("dgCMatrix_qrsol requires a 'tall' rectangular matrix"));`	162	`error(_("dgCMatrix_qrsol requires a 'tall' rectangular matrix"));`
162	`if (!isReal(ycp) \|\| LENGTH(ycp) != xc->m)`	163	`if (!isReal(ycp) \|\| LENGTH(ycp) != xc->m)`
163	`error(_("Dimensions of system to be solved are inconsistent"));`	164	`error(_("Dimensions of system to be solved are inconsistent"));`
164	`if (!cs_qrsol(/order/ 1, xc, REAL(ycp)))`	165	`if (!cs_qrsol(/order/ 1, xc, REAL(ycp)))`
165	`error(_("cs_qrsol failed"));`	166	`error(_("cs_qrsol failed"));`
166	`Free(xc);`	167
167	`UNPROTECT(1);`	168	`UNPROTECT(1);`
168	`return ycp;`	169	`return ycp;`
169	`}`	170	`}`
170		171
171	`/* Modified version of Tim Davis's cs_qr_mex.c file for MATLAB */`	172	`/* Modified version of Tim Davis's cs_qr_mex.c file for MATLAB */`
172	`SEXP dgCMatrix_QR(SEXP Ap, SEXP order)`	173	`SEXP dgCMatrix_QR(SEXP Ap, SEXP order)`
173	`{`	174	`{`
174	`SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseQR")));`	175	`SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseQR")));`
175	`cs A = Matrix_as_cs(Ap), D;`	176	`CSP A = AS_CSP(Ap), D;`
176	`css *S;`	177	`css *S;`
177	`csn *N;`	178	`csn *N;`
178	`int m = A->m, n = A->n, ord = asLogical(order) ? 3 : 0, *p;`	179	`int m = A->m, n = A->n, ord = asLogical(order) ? 3 : 0, *p;`
-		180	`R_CheckStack();`
179		181
180	`if (m < n) error("A must have # rows >= # columns") ;`	182	`if (m < n) error("A must have # rows >= # columns") ;`
181	`S = cs_sqr(ord, A, 1); /* symbolic QR ordering & analysis*/`	183	`S = cs_sqr(ord, A, 1); /* symbolic QR ordering & analysis*/`
182	`if (!S) error("cs_sqr failed");`	184	`if (!S) error("cs_sqr failed");`
183	`N = cs_qr(A, S); /* numeric QR factorization */`	185	`N = cs_qr(A, S); /* numeric QR factorization */`
Line 216...	Line 218...
216		218
217	`/* Modified version of Tim Davis's cs_lu_mex.c file for MATLAB */`	219	`/* Modified version of Tim Davis's cs_lu_mex.c file for MATLAB */`
218	`SEXP dgCMatrix_LU(SEXP Ap, SEXP orderp, SEXP tolp)`	220	`SEXP dgCMatrix_LU(SEXP Ap, SEXP orderp, SEXP tolp)`
219	`{`	221	`{`
220	`SEXP ans = get_factors(Ap, "LU");`	222	`SEXP ans = get_factors(Ap, "LU");`
221	`cs A, D;`	223	`CSP A = AS_CSP(Ap), D;`
222	`css *S;`	224	`css *S;`
223	`csn *N;`	225	`csn *N;`
224	`int n, order = asInteger(orderp), *p;`	226	`int n, order = asInteger(orderp), *p;`
225	`double tol = asReal(tolp);`	227	`double tol = asReal(tolp);`
-		228	`R_CheckStack();`
226		229
227	`/* FIXME: dgCMatrix_LU should check ans for consistency in`	230	`/* FIXME: dgCMatrix_LU should check ans for consistency in`
228	`* permutation type with the requested value - Should have two`	231	`* permutation type with the requested value - Should have two`
229	`* classes or two different names in the factors list for LU with`	232	`* classes or two different names in the factors list for LU with`
230	`* permuted columns or not. */`	233	`* permuted columns or not. */`
231		234
232	`if (ans != R_NilValue) return ans;`	235	`if (ans != R_NilValue) return ans;`
233	`ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseLU")));`	236	`ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseLU")));`
234	`A = Matrix_as_cs(Ap);`	-
235	`n = A->n;`	237	`n = A->n;`
236	`if (A->m != n)`	238	`if (A->m != n)`
237	`error("LU decomposition applies only to square matrices");`	239	`error("LU decomposition applies only to square matrices");`
238	`if (order) { /* not using natural order */`	240	`if (order) { /* not using natural order */`
239	`order = (tol == 1) ? 2 /* amd(S'S) w/dense rows or I /`	241	`order = (tol == 1) ? 2 /* amd(S'S) w/dense rows or I /`
Line 263...	Line 265...
263	`Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"),`	265	`Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"),`
264	`INTSXP, n)), S->q, n);`	266	`INTSXP, n)), S->q, n);`
265	`cs_nfree(N);`	267	`cs_nfree(N);`
266	`cs_sfree(S);`	268	`cs_sfree(S);`
267	`cs_free(p);`	269	`cs_free(p);`
268	`Free(A);`	-
269	`UNPROTECT(1);`	270	`UNPROTECT(1);`
270	`return set_factors(Ap, ans, "LU");`	271	`return set_factors(Ap, ans, "LU");`
271	`}`	272	`}`
272		273
273	`SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b)`	274	`SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b)`
274	`{`	275	`{`
275	`/* b is dense or NULL [ <--> solve(A) */`	276	`/* b is dense or NULL [ <--> solve(A) */`
276	`SEXP lu = dgCMatrix_LU(Ap, ScalarLogical(1), ScalarReal(1));`	277	`SEXP lu = dgCMatrix_LU(Ap, ScalarLogical(1), ScalarReal(1));`
277	`SEXP qslot = GET_SLOT(lu, install("q"));`	278	`SEXP qslot = GET_SLOT(lu, install("q"));`
278	`cs *L = Matrix_as_cs(GET_SLOT(lu, install("L"))),`	279	`CSP L = AS_CSP(GET_SLOT(lu, install("L"))),`
279	`*U = Matrix_as_cs(GET_SLOT(lu, install("U")));`	280	`U = AS_CSP(GET_SLOT(lu, install("U")));`
280	`SEXP ans = PROTECT( !isNull(b) ? dup_mMatrix_as_dgeMatrix(b)`	281	`SEXP ans = PROTECT( !isNull(b) ? dup_mMatrix_as_dgeMatrix(b)`
281	`: new_dgeMatrix(U->n, U->n));`	282	`: new_dgeMatrix(U->n, U->n));`
282	`int *bdims = INTEGER(GET_SLOT(ans, Matrix_DimSym));`	283	`int *bdims = INTEGER(GET_SLOT(ans, Matrix_DimSym));`
283	`int j, n = bdims[0], nrhs = bdims[1];`	284	`int j, n = bdims[0], nrhs = bdims[1];`
284	`int *p = INTEGER(GET_SLOT(lu, Matrix_pSym)),`	285	`int *p = INTEGER(GET_SLOT(lu, Matrix_pSym)),`
285	`q = LENGTH(qslot) ? INTEGER(qslot) : (int ) NULL;`	286	`q = LENGTH(qslot) ? INTEGER(qslot) : (int ) NULL;`
286	`double *ax = REAL(GET_SLOT(ans, Matrix_xSym)),`	287	`double *ax = REAL(GET_SLOT(ans, Matrix_xSym)),`
287	`*x = Calloc(n, double);`	288	`*x = Calloc(n, double);`
-		289	`R_CheckStack();`
288		290
289	`if (U->n != n \|\| nrhs < 1 \|\| n < 1)`	291	`if (U->n != n \|\| nrhs < 1 \|\| n < 1)`
290	`error(_("Dimensions of system to be solved are inconsistent"));`	292	`error(_("Dimensions of system to be solved are inconsistent"));`
291	`for (j = 0; j < nrhs; j++) {`	293	`for (j = 0; j < nrhs; j++) {`
292	`if(!isNull(b))`	294	`if(!isNull(b))`
Line 300...	Line 302...
300	`if (q) /* b(q) = x */`	302	`if (q) /* b(q) = x */`
301	`cs_ipvec(q, x, ax + j * n, n);`	303	`cs_ipvec(q, x, ax + j * n, n);`
302	`else`	304	`else`
303	`Memcpy(ax + j * n, x, n);`	305	`Memcpy(ax + j * n, x, n);`
304	`}`	306	`}`
305	`Free(L); Free(U); Free(x);`	307	`Free(x);`
306	`UNPROTECT(1);`	308	`UNPROTECT(1);`
307	`return ans;`	309	`return ans;`
308	`}`	310	`}`
309		311
310	`SEXP dgCMatrix_cholsol(SEXP x, SEXP y)`	312	`SEXP dgCMatrix_cholsol(SEXP x, SEXP y)`

The R Project SVN R-packages

(root)/branches/Matrix-APIchange/src/dgCMatrix.c – Rev 4530 → 4543