The R Project SVN R-packages

Rev

Details | Last modification | View Log | RSS feed

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