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

Parent Directory | Revision Log

---

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

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 :			}

---

root@r-forge.r-project.org	ViewVC Help
Powered by ViewVC 1.0.0