1 : |
bates |
10 |
/* Sparse triangular matrices */
|
2 : |
bates |
478 |
#include "dtCMatrix.h"
|
3 : |
bates |
10 |
|
4 : |
|
|
SEXP tsc_validate(SEXP x)
|
5 : |
|
|
{
|
6 : |
|
|
return ScalarLogical(1);
|
7 : |
|
|
}
|
8 : |
|
|
|
9 : |
|
|
SEXP tsc_transpose(SEXP x)
|
10 : |
|
|
{
|
11 : |
|
|
SEXP
|
12 : |
bates |
478 |
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 : |
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 : |
maechler |
534 |
if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U')
|
22 : |
|
|
SET_SLOT(ans, Matrix_uploSym, mkString("L"));
|
23 : |
bates |
10 |
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 : |
bates |
478 |
SEXP tsc_to_dgTMatrix(SEXP x)
|
38 : |
bates |
70 |
{
|
39 : |
|
|
SEXP ans;
|
40 : |
maechler |
534 |
if (CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0] != 'U')
|
41 : |
bates |
478 |
ans = csc_to_dgTMatrix(x);
|
42 : |
bates |
70 |
else { /* unit triangular matrix */
|
43 : |
maechler |
534 |
SEXP islot = GET_SLOT(x, Matrix_iSym),
|
44 : |
bates |
70 |
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 : |
maechler |
534 |
|
52 : |
bates |
478 |
ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix")));
|
53 : |
bates |
70 |
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 : |
maechler |
534 |
/**
|
74 : |
bates |
358 |
* Derive the column pointer vector for the inverse of L from the parent array
|
75 : |
maechler |
534 |
*
|
76 : |
bates |
358 |
* @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 : |
maechler |
534 |
*
|
81 : |
bates |
358 |
* @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 : |
maechler |
534 |
/**
|
99 : |
bates |
358 |
* Derive the row index array for the inverse of L from the parent array
|
100 : |
maechler |
534 |
*
|
101 : |
bates |
358 |
* @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 : |
maechler |
534 |
|
115 : |
bates |
338 |
SEXP Parent_inverse(SEXP par, SEXP unitdiag)
|
116 : |
|
|
{
|
117 : |
bates |
478 |
SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
|
118 : |
bates |
338 |
int *ap, *ai, *dims, *pr = INTEGER(par),
|
119 : |
|
|
countDiag = 1 - asLogical(unitdiag),
|
120 : |
bates |
367 |
j, n = length(par), nnz;
|
121 : |
bates |
338 |
double *ax;
|
122 : |
maechler |
534 |
|
123 : |
bates |
338 |
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 : |
maechler |
534 |
SET_SLOT(ans, Matrix_uploSym, mkString("L"));
|
136 : |
|
|
SET_SLOT(ans, Matrix_diagSym, (countDiag ? mkString("N") : mkString("U")));
|
137 : |
bates |
358 |
parent_inv_ai(n, countDiag, pr, ai);
|
138 : |
bates |
338 |
UNPROTECT(1);
|
139 : |
|
|
return ans;
|
140 : |
|
|
}
|