SCM

SCM Repository

[matrix] Annotation of /pkg/src/lsCMatrix.c
ViewVC logotype

Annotation of /pkg/src/lsCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 737 #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)),
64 :     *P, *Pinv, *Parent, *Lp;
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 :     }

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge