SCM

SCM Repository

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

Annotation of /pkg/Matrix/src/dtCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 358 - (view) (download) (as text)
Original Path: pkg/src/tscMatrix.c

1 : bates 10 /* Sparse triangular matrices */
2 :     #include "tscMatrix.h"
3 :    
4 :     SEXP tsc_validate(SEXP x)
5 :     {
6 :     return ScalarLogical(1);
7 :     }
8 :    
9 :     SEXP tsc_transpose(SEXP x)
10 :     {
11 :     SEXP
12 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix"))),
13 :     islot = GET_SLOT(x, Matrix_iSym);
14 :     int nnz = length(islot),
15 : bates 338 *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
16 : bates 10
17 : bates 338
18 :     SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
19 :     adims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
20 : bates 10 adims[0] = xdims[1]; adims[1] = xdims[0];
21 :     if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'U')
22 :     SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));
23 :     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));
24 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
25 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
26 :     csc_components_transpose(xdims[0], xdims[1], nnz,
27 :     INTEGER(GET_SLOT(x, Matrix_pSym)),
28 :     INTEGER(islot),
29 :     REAL(GET_SLOT(x, Matrix_xSym)),
30 :     INTEGER(GET_SLOT(ans, Matrix_pSym)),
31 :     INTEGER(GET_SLOT(ans, Matrix_iSym)),
32 :     REAL(GET_SLOT(ans, Matrix_xSym)));
33 :     UNPROTECT(1);
34 :     return ans;
35 :     }
36 : bates 70
37 :     SEXP tsc_to_triplet(SEXP x)
38 :     {
39 :     SEXP ans;
40 :     if (toupper(CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0]) != 'U')
41 :     ans = csc_to_triplet(x);
42 :     else { /* unit triangular matrix */
43 :     SEXP islot = GET_SLOT(x, Matrix_iSym),
44 :     pslot = GET_SLOT(x, Matrix_pSym);
45 :     int *ai, *aj, j,
46 :     n = length(pslot) - 1,
47 :     nod = length(islot),
48 :     nout = n + nod,
49 :     *p = INTEGER(pslot);
50 :     double *ax;
51 :    
52 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tripletMatrix")));
53 :     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
54 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));
55 :     ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
56 :     Memcpy(ai, INTEGER(islot), nod);
57 :     SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nout));
58 :     aj = INTEGER(GET_SLOT(ans, Matrix_jSym));
59 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nout));
60 :     ax = REAL(GET_SLOT(ans, Matrix_xSym));
61 :     Memcpy(ax, REAL(GET_SLOT(x, Matrix_xSym)), nod);
62 :     for (j = 0; j < n; j++) {
63 :     int jj, npj = nod + j, p2 = p[j+1];
64 :     aj[npj] = ai[npj] = j;
65 :     ax[npj] = 1.;
66 :     for (jj = p[j]; jj < p2; jj++) aj[jj] = j;
67 :     }
68 :     UNPROTECT(1);
69 :     }
70 :     return ans;
71 :     }
72 : bates 338
73 : bates 358 /**
74 :     * Derive the column pointer vector for the inverse of L from the parent array
75 :     *
76 :     * @param n length of parent array
77 :     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
78 :     * @param pr parent vector describing the elimination tree
79 :     * @param ap array of length n+1 to be filled with the column pointers
80 :     *
81 :     * @return the number of non-zero entries (ap[n])
82 :     */
83 :     int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])
84 :     {
85 :     int *sz = Calloc(n, int), j;
86 :    
87 :     for (j = n - 1; j >= 0; j--) {
88 :     int parent = pr[j];
89 :     sz[j] = (parent < 0) ? countDiag : (1 + sz[parent]);
90 :     }
91 :     ap[0] = 0;
92 :     for (j = 0; j < n; j++)
93 :     ap[j+1] = ap[j] + sz[j];
94 :     Free(sz);
95 :     return ap[n];
96 :     }
97 :    
98 :     /**
99 :     * Derive the row index array for the inverse of L from the parent array
100 :     *
101 :     * @param n length of parent array
102 :     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
103 :     * @param pr parent vector describing the elimination tree
104 :     * @param ai row index vector of length ap[n]
105 :     */
106 :     void parent_inv_ai(int n, int countDiag, const int pr[], int ai[])
107 :     {
108 :     int j, k, pos = 0;
109 :     for (j = 0; j < n; j++) {
110 :     if (countDiag) ai[pos++] = j;
111 :     for (k = pr[j]; k >= 0; k = pr[k]) ai[pos++] = k;
112 :     }
113 :     }
114 :    
115 : bates 338 SEXP Parent_inverse(SEXP par, SEXP unitdiag)
116 :     {
117 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix")));
118 :     int *ap, *ai, *dims, *pr = INTEGER(par),
119 :     countDiag = 1 - asLogical(unitdiag),
120 : bates 358 j, k, n = length(par), nnz;
121 : bates 338 double *ax;
122 :    
123 :     if (!isInteger(par)) error("par argument must be an integer vector");
124 :     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));
125 :     ap = INTEGER(GET_SLOT(ans, Matrix_pSym));
126 : bates 358 nnz = parent_inv_ap(n, countDiag, pr, ap);
127 : bates 338 SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
128 :     ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
129 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
130 :     ax = REAL(GET_SLOT(ans, Matrix_xSym));
131 :     for (j = 0; j < nnz; j++) ax[j] = 1.;
132 :     SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
133 :     dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
134 :     dims[0] = dims[1] = n;
135 :     SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));
136 :     SET_SLOT(ans, Matrix_diagSym,
137 :     (countDiag ? ScalarString(mkChar("N")) :
138 :     ScalarString(mkChar("U"))));
139 : bates 358 parent_inv_ai(n, countDiag, pr, ai);
140 : bates 338 UNPROTECT(1);
141 :     return ans;
142 :     }

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