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 2210 - (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 1251 #include "cs_utils.h"
4 : bates 10
5 : maechler 1968 #define RETURN(_CH_) UNPROTECT(1); return (_CH_);
6 :    
7 :     /* This is used for *BOTH* triangular and symmetric Csparse: */
8 : maechler 1710 SEXP tCMatrix_validate(SEXP x)
9 :     {
10 :     SEXP val = xCMatrix_validate(x);/* checks x slot */
11 :     if(isString(val))
12 :     return(val);
13 :     else {
14 :     SEXP
15 :     islot = GET_SLOT(x, Matrix_iSym),
16 :     pslot = GET_SLOT(x, Matrix_pSym);
17 :     int uploT = (*uplo_P(x) == 'U'),
18 :     k, nnz = length(islot),
19 :     *xi = INTEGER(islot),
20 :     *xj = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
21 :    
22 :     expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xj);
23 :    
24 :     /* Maybe FIXME: ">" should be ">=" for diag = 'U' (uplo = 'U') */
25 :     if(uploT) {
26 :     for (k = 0; k < nnz; k++)
27 :     if(xi[k] > xj[k]) {
28 : dmbates 2208 RETURN(mkString(_("uplo='U' must not have sparse entries below the diagonal")));
29 : maechler 1710 }
30 :     }
31 :     else {
32 :     for (k = 0; k < nnz; k++)
33 :     if(xi[k] < xj[k]) {
34 : dmbates 2208 RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
35 : maechler 1710 }
36 :     }
37 :    
38 :     RETURN(ScalarLogical(1));
39 :     }
40 :     }
41 : maechler 1968
42 :     /* This is used for *BOTH* triangular and symmetric Rsparse: */
43 :     SEXP tRMatrix_validate(SEXP x)
44 :     {
45 :     SEXP val = xRMatrix_validate(x);/* checks x slot */
46 :     if(isString(val))
47 :     return(val);
48 :     else {
49 :     SEXP
50 :     jslot = GET_SLOT(x, Matrix_jSym),
51 :     pslot = GET_SLOT(x, Matrix_pSym);
52 :     int uploT = (*uplo_P(x) == 'U'),
53 :     k, nnz = length(jslot),
54 :     *xj = INTEGER(jslot),
55 :     *xi = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
56 :    
57 :     expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xi);
58 :    
59 :     /* Maybe FIXME: ">" should be ">=" for diag = 'U' (uplo = 'U') */
60 :     if(uploT) {
61 :     for (k = 0; k < nnz; k++)
62 :     if(xi[k] > xj[k]) {
63 : dmbates 2208 RETURN(mkString(_("uplo='U' must not have sparse entries below the diagonal")));
64 : maechler 1968 }
65 :     }
66 :     else {
67 :     for (k = 0; k < nnz; k++)
68 :     if(xi[k] < xj[k]) {
69 : dmbates 2208 RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
70 : maechler 1968 }
71 :     }
72 :    
73 :     RETURN(ScalarLogical(1));
74 :     }
75 :     }
76 :    
77 : maechler 1710 #undef RETURN
78 :    
79 : bates 1248 SEXP dtCMatrix_solve(SEXP a)
80 :     {
81 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
82 : dmbates 2210 CSP A = AS_CSP(Csparse_diagU2N(a));
83 : bates 1251 int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
84 : maechler 1960 bnz = 10 * A->n, /* initial estimate of nnz in b */
85 :     lo = uplo_P(a)[0] == 'L', top;
86 :     /* These arrays must use Calloc because of possible Realloc */
87 : maechler 1798 int *ti = Calloc(bnz, int), p, j, nz, pos = 0;
88 : maechler 1960 double *tx = Calloc(bnz, double);
89 : maechler 1798 cs *u = cs_spalloc(A->n, 1,1,1,0); /* Sparse unit vector */
90 : maechler 1960 double *wrk = Alloca(A->n, double);
91 :     int *xi = Alloca(2*A->n, int); /* for cs_reach */
92 :     R_CheckStack();
93 : bates 1251
94 : maechler 1968 slot_dup(ans, a, Matrix_DimSym);
95 : maechler 1736 SET_DimNames(ans, a);
96 : maechler 1968 slot_dup(ans, a, Matrix_uploSym);
97 :     slot_dup(ans, a, Matrix_diagSym);
98 : maechler 1798 /* initialize the "constant part" of the sparse unit vector */
99 :     u->x[0] = 1.;
100 :     u->p[0] = 0; u->p[1] = 1;
101 : bates 1251 bp[0] = 0;
102 :     for (j = 0; j < A->n; j++) {
103 : maechler 1798 u->i[0] = j; /* u := j'th unit vector */
104 :     /* (wrk[top:n],xi[top:n]) := A^{-1} u : */
105 :     top = cs_spsolve (A, u, 0, xi, wrk, 0, lo);
106 :     nz = A->n - top;
107 : bates 1251 bp[j + 1] = nz + bp[j];
108 :     if (bp[j + 1] > bnz) {
109 :     while (bp[j + 1] > bnz) bnz *= 2;
110 :     ti = Realloc(ti, bnz, int);
111 :     tx = Realloc(tx, bnz, double);
112 :     }
113 : maechler 1798 if (lo)
114 :     for(p = top; p < A->n; p++, pos++) {
115 :     ti[pos] = xi[p];
116 :     tx[pos] = wrk[xi[p]];
117 :     }
118 :     else /* upper triagonal */
119 :     for(p = A->n - 1; p >= top; p--, pos++) {
120 :     ti[pos] = xi[p];
121 :     tx[pos] = wrk[xi[p]];
122 :     }
123 : bates 1251 }
124 :     nz = bp[A->n];
125 : maechler 1798 Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
126 :     Memcpy( REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
127 : bates 1251
128 : maechler 1960 Free(ti); Free(tx); cs_spfree(u);
129 : bates 1251 UNPROTECT(1);
130 :     return ans;
131 :     }
132 :    
133 :     SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
134 :     {
135 :     int cl = asLogical(classed);
136 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
137 : dmbates 2210 CSP A = AS_CSP(Csparse_diagU2N(a));
138 : bates 1251 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
139 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
140 :     getAttrib(b, R_DimSymbol));
141 :     int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
142 :     double *bx;
143 : maechler 1960 R_CheckStack();
144 : bates 1251
145 :     if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
146 :     error(_("Dimensions of system to be solved are inconsistent"));
147 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
148 : maechler 1736 /* FIXME: copy dimnames or Dimnames as well */
149 : bates 1251 bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
150 :     REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
151 :     for (j = 0; j < nrhs; j++)
152 :     lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
153 :     UNPROTECT(1);
154 :     return ans;
155 :     }
156 : bates 1262
157 : dmbates 2210 SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)
158 : bates 1262 {
159 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
160 : dmbates 2210 CSP A = AS_CSP(Csparse_diagU2N(a)), B = AS_CSP(Csparse_diagU2N(b));
161 :     int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
162 :     xnz = 10 * B->p[B->n], /* initial estimate of nnz in x */
163 :     lo = uplo_P(a)[0] == 'L', top;
164 :     /* These arrays must use Calloc because of possible Realloc */
165 :     int *ti = Calloc(xnz, int), p, j, nz, pos = 0;
166 :     double *tx = Calloc(xnz, double);
167 :     double *wrk = Alloca(A->n, double);
168 :     int *xi = Alloca(2*A->n, int); /* for cs_reach */
169 : maechler 1960 R_CheckStack();
170 : dmbates 2210
171 :     if (A->m != A->n || B->n < 1 || A->n < 1 || A->n != B->m)
172 :     error(_("Dimensions of system to be solved are inconsistent"));
173 :     slot_dup(ans, b, Matrix_DimSym);
174 :     SET_DimNames(ans, b);
175 :     xp[0] = 0;
176 :     for (j = 0; j < B->n; j++) {
177 :     /* (wrk[top:n],xi[top:n]) := A^{-1} B[,j] */
178 :     top = cs_spsolve (A, B, j, xi, wrk, 0, lo);
179 :     nz = A->n - top;
180 :     xp[j + 1] = nz + xp[j];
181 :     if (xp[j + 1] > xnz) {
182 :     while (xp[j + 1] > xnz) xnz *= 2;
183 :     ti = Realloc(ti, xnz, int);
184 :     tx = Realloc(tx, xnz, double);
185 : bates 1262 }
186 : dmbates 2210 if (lo)
187 :     for(p = top; p < A->n; p++, pos++) {
188 :     ti[pos] = xi[p];
189 :     tx[pos] = wrk[xi[p]];
190 :     }
191 :     else /* upper triagonal */
192 :     for(p = A->n - 1; p >= top; p--, pos++) {
193 :     ti[pos] = xi[p];
194 :     tx[pos] = wrk[xi[p]];
195 :     }
196 : bates 1262 }
197 : dmbates 2210 nz = xp[A->n];
198 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
199 :     Memcpy( REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
200 :    
201 :     Free(ti); Free(tx);
202 : bates 1262 UNPROTECT(1);
203 :     return ans;
204 :     }

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