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 896 - (view) (download) (as text)

1 : bates 737 #include "lsCMatrix.h"
2 :    
3 : maechler 896 /**
4 : bates 737 * Check the validity of the slots of an lsCMatrix object
5 : maechler 896 *
6 : bates 737 * @param x Pointer to an lsCMatrix object
7 : maechler 896 *
8 : bates 737 * @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 : maechler 896 SEXP val = symmetricMatrix_validate(x);
14 :     if(isString(val))
15 :     return(val);
16 :     else {
17 :     /* FIXME needed? ltC* inherits from lgC* which does this in validate*/
18 :     SEXP pslot = GET_SLOT(x, Matrix_pSym),
19 :     islot = GET_SLOT(x, Matrix_iSym);
20 :     int
21 :     ncol = length(pslot) - 1,
22 :     *xp = INTEGER(pslot),
23 :     *xi = INTEGER(islot);
24 : bates 737
25 : maechler 896 if (csc_unsorted_columns(ncol, xp, xi))
26 :     csc_sort_columns(ncol, xp, xi, (double *) NULL);
27 :    
28 :     return ScalarLogical(1);
29 :     }
30 : bates 737 }
31 :    
32 : maechler 896 /**
33 : bates 737 * Transpose an lsCMatrix
34 : maechler 896 *
35 : bates 737 * @param x Pointer to an lsCMatrix object
36 : maechler 896 *
37 : bates 737 * @return the transpose of x. It represents the same matrix but is
38 :     * stored in the opposite triangle.
39 :     */
40 :     SEXP lsCMatrix_trans(SEXP x)
41 :     {
42 :     SEXP Xi = GET_SLOT(x, Matrix_iSym),
43 :     xDim = GET_SLOT(x, Matrix_DimSym),
44 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lsCMatrix")));
45 :     int n = INTEGER(xDim)[0], nz = length(Xi);
46 :     int *xj = expand_cmprPt(n, INTEGER(GET_SLOT(x, Matrix_pSym)),
47 :     Calloc(nz, int)),
48 :     *xi = Memcpy(Calloc(nz, int), Xi, nz);
49 :     int up = CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U';
50 :    
51 :     SET_SLOT(ans, Matrix_DimSym, duplicate(xDim));
52 :     SET_SLOT(ans, Matrix_DimNamesSym,
53 :     duplicate(GET_SLOT(x, Matrix_DimNamesSym)));
54 :     SET_SLOT(ans, Matrix_uploSym, mkString(up ? "L" : "U"));
55 :     make_upper_triangular(up ? xj : xi, up ? xi : xj, nz);
56 :     triplet_to_col(n, n, nz, xi, xj, (double *) NULL,
57 :     INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1)),
58 :     INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)),
59 :     (double *) NULL);
60 :     Free(xj); Free(xi);
61 :     UNPROTECT(1);
62 :     return ans;
63 :     }
64 :    
65 : maechler 896 /**
66 : bates 883 * Create the symbolic Cholesky decomposition of an lsCMatrix object
67 : maechler 896 *
68 : bates 883 * @param x Pointer to an lsCMatrix object
69 :     * @param pivot Pointer to a scalar logical indicating if a
70 :     * fill-reducing permutation should be determined
71 : maechler 896 *
72 : bates 883 * @return an lCholCMatrix object
73 :     */
74 : bates 737 SEXP lsCMatrix_chol(SEXP x, SEXP pivot)
75 :     {
76 :     int piv = asLogical(pivot);
77 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lCholCMatrix")));
78 :     int j, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
79 :     int *Xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
80 :     *Xp = INTEGER(GET_SLOT(x, Matrix_pSym)),
81 : bates 765 *P, *Pinv = (int *) NULL, *Parent, *Lp;
82 : bates 737 double *D = Calloc(n, double), *Tx, *Xx = Calloc(Xp[n], double);
83 :    
84 :     if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] != 'U')
85 :     error(_("Must have uplo == 'U' in x argument to lsCMatrix_chol"));
86 :     SET_SLOT(ans, Matrix_uploSym, mkString("L"));
87 :     SET_SLOT(ans, Matrix_diagSym, mkString("U"));
88 :     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
89 :     SET_SLOT(ans, Matrix_DimNamesSym, duplicate(GET_SLOT(x, Matrix_DimNamesSym)));
90 :     P = INTEGER(ALLOC_SLOT(ans, Matrix_permSym, INTSXP, n));
91 :     if (piv) {
92 :     Pinv = Calloc(n, int);
93 :     ssc_metis_order(n, Xp, Xi, P, Pinv);
94 :     } else {
95 :     int i;
96 :     for (i = 0; i < n; i++) P[i] = i;
97 :     }
98 :     Lp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
99 :     Parent = INTEGER(ALLOC_SLOT(ans, Matrix_ParentSym, INTSXP, n));
100 :     R_ldl_symbolic(n, Xp, Xi, Lp, Parent,
101 :     (piv) ? P : (int *) NULL, (piv) ? Pinv : (int *) NULL);
102 :     /* Decompose the identity to get Li */
103 :     for (j = 0; j < n; j++) { /* Create an identity from Xp, Xi and Xx */
104 :     int ii, ii2 = Xp[j + 1];
105 :     for (ii = Xp[j]; ii < ii2; ii++)
106 :     Xx[ii] = (Xi[ii] == j) ? 1. : 0.;
107 :     }
108 :     Tx = Calloc(Lp[n], double);
109 :     R_ldl_numeric(n, Xp, Xi, Xx, Lp, Parent,
110 :     INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, Lp[n])),
111 :     Tx, D, (piv) ? P : (int *) NULL,
112 :     (piv) ? Pinv : (int *) NULL);
113 :     if (piv) Free(Pinv);
114 :     Free(Xx); Free(Tx); Free(D);
115 :     UNPROTECT(1);
116 :     return ans;
117 :     }

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