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 951 - (view) (download) (as text)
Original Path: pkg/src/dtCMatrix.c

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

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