The R Project SVN – R-packages – Diff – /trunk/Matrix/src/lsCMatrix.c


#include "lsCMatrix.h"
 
/** 
 * Check the validity of the slots of an lsCMatrix object
 * 
 * @param x Pointer to an lsCMatrix object
 * 
 * @return an SEXP that is either TRUE or a character string
 * describing the way in which the object failed the validity check
 */
SEXP lsCMatrix_validate(SEXP x)
{
    SEXP val = check_scalar_string(GET_SLOT(x, Matrix_uploSym),
				   "LU", "uplo");
    int *Dim = INTEGER(GET_SLOT(x, Matrix_DimSym));
 
    if (isString(val)) return val;
    if (Dim[0] != Dim[1])
	return mkString(_("Symmetric matrix must be square"));
    csc_check_column_sorting(x);
    return ScalarLogical(1);
}
 
/** 
 * Transpose an lsCMatrix
 * 
 * @param x Pointer to an lsCMatrix object
 * 
 * @return the transpose of x.  It represents the same matrix but is
 * stored in the opposite triangle.
 */
SEXP lsCMatrix_trans(SEXP x)
{
    SEXP Xi = GET_SLOT(x, Matrix_iSym),
	xDim = GET_SLOT(x, Matrix_DimSym),
	ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lsCMatrix")));
    int n = INTEGER(xDim)[0], nz = length(Xi);
    int *xj = expand_cmprPt(n, INTEGER(GET_SLOT(x, Matrix_pSym)),
			    Calloc(nz, int)),
	*xi = Memcpy(Calloc(nz, int), Xi, nz);
    int up = CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U';
 
    SET_SLOT(ans, Matrix_DimSym, duplicate(xDim));
    SET_SLOT(ans, Matrix_DimNamesSym,
	     duplicate(GET_SLOT(x, Matrix_DimNamesSym)));
    SET_SLOT(ans, Matrix_uploSym, mkString(up ? "L" : "U"));
    make_upper_triangular(up ? xj : xi, up ? xi : xj, nz);
    triplet_to_col(n, n, nz, xi, xj, (double *) NULL,
		   INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP,  n + 1)),
		   INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)),
		   (double *) NULL);
    Free(xj); Free(xi);
    UNPROTECT(1);
    return ans;
}
 
SEXP lsCMatrix_chol(SEXP x, SEXP pivot)
{
    int piv = asLogical(pivot);
    SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lCholCMatrix")));
    int j, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
    int *Xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
	*Xp = INTEGER(GET_SLOT(x, Matrix_pSym)),
	*P, *Pinv = (int *) NULL, *Parent, *Lp;
    double *D = Calloc(n, double), *Tx, *Xx = Calloc(Xp[n], double);
 
    if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] != 'U')
	error(_("Must have uplo == 'U' in x argument to lsCMatrix_chol"));
    SET_SLOT(ans, Matrix_uploSym, mkString("L"));
    SET_SLOT(ans, Matrix_diagSym, mkString("U"));
    SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
    SET_SLOT(ans, Matrix_DimNamesSym, duplicate(GET_SLOT(x, Matrix_DimNamesSym)));
    P = INTEGER(ALLOC_SLOT(ans, Matrix_permSym, INTSXP, n));
    if (piv) {
	Pinv = Calloc(n, int);
	ssc_metis_order(n, Xp, Xi, P, Pinv);
    } else {
	int i;
	for (i = 0; i < n; i++) P[i] = i;
    }
    Lp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
    Parent = INTEGER(ALLOC_SLOT(ans, Matrix_ParentSym, INTSXP, n));
    R_ldl_symbolic(n, Xp, Xi, Lp, Parent,
		   (piv) ? P : (int *) NULL, (piv) ? Pinv : (int *) NULL);
    				/* Decompose the identity to get Li */
    for (j = 0; j < n; j++) {	/* Create an identity from Xp, Xi and Xx */
	int ii, ii2 = Xp[j + 1];
	for (ii = Xp[j]; ii < ii2; ii++)
	    Xx[ii] = (Xi[ii] == j) ? 1. : 0.;
    }
    Tx = Calloc(Lp[n], double);
    R_ldl_numeric(n, Xp, Xi, Xx, Lp, Parent,
		  INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, Lp[n])),
		  Tx, D, (piv) ? P : (int *) NULL,
		  (piv) ? Pinv : (int *) NULL);
    if (piv) Free(Pinv);
    Free(Xx); Free(Tx); Free(D);
    UNPROTECT(1);
    return ans;
}

The R Project SVN R-packages

(root)/trunk/Matrix/src/lsCMatrix.c @ 1663 – Rev 1535 → 1576

Rev 1535		Rev 1576
1	`#include "lsCMatrix.h"`	1	`#include "lsCMatrix.h"`
2		2
3	`/**`	3	`/**`
4	`* Check the validity of the slots of an lsCMatrix object`	4	`* Check the validity of the slots of an lsCMatrix object`
5	`*`	5	`*`
6	`* @param x Pointer to an lsCMatrix object`	6	`* @param x Pointer to an lsCMatrix object`
7	`*`	7	`*`
8	`* @return an SEXP that is either TRUE or a character string`	8	`* @return an SEXP that is either TRUE or a character string`
9	`* describing the way in which the object failed the validity check`	9	`* describing the way in which the object failed the validity check`
10	`*/`	10	`*/`
11	`SEXP lsCMatrix_validate(SEXP x)`	11	`SEXP lsCMatrix_validate(SEXP x)`
12	`{`	12	`{`
13	`SEXP val = check_scalar_string(GET_SLOT(x, Matrix_uploSym),`	13	`SEXP val = check_scalar_string(GET_SLOT(x, Matrix_uploSym),`
14	`"LU", "uplo");`	14	`"LU", "uplo");`
15	`int *Dim = INTEGER(GET_SLOT(x, Matrix_DimSym));`	15	`int *Dim = INTEGER(GET_SLOT(x, Matrix_DimSym));`
16		16
17	`if (isString(val)) return val;`	17	`if (isString(val)) return val;`
18	`if (Dim[0] != Dim[1])`	18	`if (Dim[0] != Dim[1])`
19	`return mkString(_("Symmetric matrix must be square"));`	19	`return mkString(_("Symmetric matrix must be square"));`
20	`csc_check_column_sorting(x);`	20	`csc_check_column_sorting(x);`
21	`return ScalarLogical(1);`	21	`return ScalarLogical(1);`
22	`}`	22	`}`
23		23
24	`/**`	24	`/**`
25	`* Transpose an lsCMatrix`	25	`* Transpose an lsCMatrix`
26	`*`	26	`*`
27	`* @param x Pointer to an lsCMatrix object`	27	`* @param x Pointer to an lsCMatrix object`
28	`*`	28	`*`
29	`* @return the transpose of x. It represents the same matrix but is`	29	`* @return the transpose of x. It represents the same matrix but is`
30	`* stored in the opposite triangle.`	30	`* stored in the opposite triangle.`
31	`*/`	31	`*/`
32	`SEXP lsCMatrix_trans(SEXP x)`	32	`SEXP lsCMatrix_trans(SEXP x)`
33	`{`	33	`{`
34	`SEXP Xi = GET_SLOT(x, Matrix_iSym),`	34	`SEXP Xi = GET_SLOT(x, Matrix_iSym),`
35	`xDim = GET_SLOT(x, Matrix_DimSym),`	35	`xDim = GET_SLOT(x, Matrix_DimSym),`
36	`ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lsCMatrix")));`	36	`ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lsCMatrix")));`
37	`int n = INTEGER(xDim)[0], nz = length(Xi);`	37	`int n = INTEGER(xDim)[0], nz = length(Xi);`
38	`int *xj = expand_cmprPt(n, INTEGER(GET_SLOT(x, Matrix_pSym)),`	38	`int *xj = expand_cmprPt(n, INTEGER(GET_SLOT(x, Matrix_pSym)),`
39	`Calloc(nz, int)),`	39	`Calloc(nz, int)),`
40	`*xi = Memcpy(Calloc(nz, int), Xi, nz);`	40	`*xi = Memcpy(Calloc(nz, int), Xi, nz);`
41	`int up = CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U';`	41	`int up = CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U';`
42		42
43	`SET_SLOT(ans, Matrix_DimSym, duplicate(xDim));`	43	`SET_SLOT(ans, Matrix_DimSym, duplicate(xDim));`
44	`SET_SLOT(ans, Matrix_DimNamesSym,`	44	`SET_SLOT(ans, Matrix_DimNamesSym,`
45	`duplicate(GET_SLOT(x, Matrix_DimNamesSym)));`	45	`duplicate(GET_SLOT(x, Matrix_DimNamesSym)));`
46	`SET_SLOT(ans, Matrix_uploSym, mkString(up ? "L" : "U"));`	46	`SET_SLOT(ans, Matrix_uploSym, mkString(up ? "L" : "U"));`
47	`make_upper_triangular(up ? xj : xi, up ? xi : xj, nz);`	47	`make_upper_triangular(up ? xj : xi, up ? xi : xj, nz);`
48	`triplet_to_col(n, n, nz, xi, xj, (double *) NULL,`	48	`triplet_to_col(n, n, nz, xi, xj, (double *) NULL,`
49	`INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1)),`	49	`INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1)),`
50	`INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)),`	50	`INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)),`
51	`(double *) NULL);`	51	`(double *) NULL);`
52	`Free(xj); Free(xi);`	52	`Free(xj); Free(xi);`
53	`UNPROTECT(1);`	53	`UNPROTECT(1);`
54	`return ans;`	54	`return ans;`
55	`}`	55	`}`
56		56
57	`SEXP lsCMatrix_chol(SEXP x, SEXP pivot)`	57	`SEXP lsCMatrix_chol(SEXP x, SEXP pivot)`
58	`{`	58	`{`
59	`int piv = asLogical(pivot);`	59	`int piv = asLogical(pivot);`
60	`SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lCholCMatrix")));`	60	`SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lCholCMatrix")));`
61	`int j, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];`	61	`int j, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];`
62	`int *Xi = INTEGER(GET_SLOT(x, Matrix_iSym)),`	62	`int *Xi = INTEGER(GET_SLOT(x, Matrix_iSym)),`
63	`*Xp = INTEGER(GET_SLOT(x, Matrix_pSym)),`	63	`*Xp = INTEGER(GET_SLOT(x, Matrix_pSym)),`
64	`P, Pinv, Parent, Lp;`	64	`P, Pinv = (int ) NULL, Parent, *Lp;`
65	`double D = Calloc(n, double), Tx, *Xx = Calloc(Xp[n], double);`	65	`double D = Calloc(n, double), Tx, *Xx = Calloc(Xp[n], double);`
66		66
67	`if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] != 'U')`	67	`if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] != 'U')`
68	`error(_("Must have uplo == 'U' in x argument to lsCMatrix_chol"));`	68	`error(_("Must have uplo == 'U' in x argument to lsCMatrix_chol"));`
69	`SET_SLOT(ans, Matrix_uploSym, mkString("L"));`	69	`SET_SLOT(ans, Matrix_uploSym, mkString("L"));`
70	`SET_SLOT(ans, Matrix_diagSym, mkString("U"));`	70	`SET_SLOT(ans, Matrix_diagSym, mkString("U"));`
71	`SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));`	71	`SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));`
72	`SET_SLOT(ans, Matrix_DimNamesSym, duplicate(GET_SLOT(x, Matrix_DimNamesSym)));`	72	`SET_SLOT(ans, Matrix_DimNamesSym, duplicate(GET_SLOT(x, Matrix_DimNamesSym)));`
73	`P = INTEGER(ALLOC_SLOT(ans, Matrix_permSym, INTSXP, n));`	73	`P = INTEGER(ALLOC_SLOT(ans, Matrix_permSym, INTSXP, n));`
74	`if (piv) {`	74	`if (piv) {`
75	`Pinv = Calloc(n, int);`	75	`Pinv = Calloc(n, int);`
76	`ssc_metis_order(n, Xp, Xi, P, Pinv);`	76	`ssc_metis_order(n, Xp, Xi, P, Pinv);`
77	`} else {`	77	`} else {`
78	`int i;`	78	`int i;`
79	`for (i = 0; i < n; i++) P[i] = i;`	79	`for (i = 0; i < n; i++) P[i] = i;`
80	`}`	80	`}`
81	`Lp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));`	81	`Lp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));`
82	`Parent = INTEGER(ALLOC_SLOT(ans, Matrix_ParentSym, INTSXP, n));`	82	`Parent = INTEGER(ALLOC_SLOT(ans, Matrix_ParentSym, INTSXP, n));`
83	`R_ldl_symbolic(n, Xp, Xi, Lp, Parent,`	83	`R_ldl_symbolic(n, Xp, Xi, Lp, Parent,`
84	`(piv) ? P : (int ) NULL, (piv) ? Pinv : (int ) NULL);`	84	`(piv) ? P : (int ) NULL, (piv) ? Pinv : (int ) NULL);`
85	`/* Decompose the identity to get Li */`	85	`/* Decompose the identity to get Li */`
86	`for (j = 0; j < n; j++) { /* Create an identity from Xp, Xi and Xx */`	86	`for (j = 0; j < n; j++) { /* Create an identity from Xp, Xi and Xx */`
87	`int ii, ii2 = Xp[j + 1];`	87	`int ii, ii2 = Xp[j + 1];`
88	`for (ii = Xp[j]; ii < ii2; ii++)`	88	`for (ii = Xp[j]; ii < ii2; ii++)`
89	`Xx[ii] = (Xi[ii] == j) ? 1. : 0.;`	89	`Xx[ii] = (Xi[ii] == j) ? 1. : 0.;`
90	`}`	90	`}`
91	`Tx = Calloc(Lp[n], double);`	91	`Tx = Calloc(Lp[n], double);`
92	`R_ldl_numeric(n, Xp, Xi, Xx, Lp, Parent,`	92	`R_ldl_numeric(n, Xp, Xi, Xx, Lp, Parent,`
93	`INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, Lp[n])),`	93	`INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, Lp[n])),`
94	`Tx, D, (piv) ? P : (int *) NULL,`	94	`Tx, D, (piv) ? P : (int *) NULL,`
95	`(piv) ? Pinv : (int *) NULL);`	95	`(piv) ? Pinv : (int *) NULL);`
96	`if (piv) Free(Pinv);`	96	`if (piv) Free(Pinv);`
97	`Free(Xx); Free(Tx); Free(D);`	97	`Free(Xx); Free(Tx); Free(D);`
98	`UNPROTECT(1);`	98	`UNPROTECT(1);`
99	`return ans;`	99	`return ans;`
100	`}`	100	`}`