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 3036 - (view) (download) (as text)

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 : bates 1251 SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
78 :     {
79 :     int cl = asLogical(classed);
80 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
81 : mmaechler 2227 CSP A = AS_CSP(a);
82 : bates 1251 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
83 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
84 :     getAttrib(b, R_DimSymbol));
85 :     int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
86 :     double *bx;
87 : maechler 1960 R_CheckStack();
88 : bates 1251
89 : mmaechler 3011 if (adims[0] != n || n != adims[1])
90 : bates 1251 error(_("Dimensions of system to be solved are inconsistent"));
91 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
92 : mmaechler 3011 // dimnames:
93 : mmaechler 3036 SEXP dn = PROTECT(allocVector(VECSXP, 2)), dn2;
94 : mmaechler 3011 SET_VECTOR_ELT(dn, 0, duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1)));
95 : mmaechler 3036 if(!cl) {
96 :     dn2 = getAttrib(b, R_DimNamesSymbol);
97 :     if(dn2 != R_NilValue) // either NULL or list(<dn1>, <dn2>)
98 :     dn2 = VECTOR_ELT(dn2, 1);
99 :     }
100 : mmaechler 3011 SET_VECTOR_ELT(dn, 1, duplicate(cl // b can be "Matrix" or not:
101 :     ? VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)
102 : mmaechler 3036 : dn2));
103 : mmaechler 3011 SET_SLOT(ans, Matrix_DimNamesSym, dn);
104 :     UNPROTECT(1);
105 :     if(n >= 1 && nrhs >=1) {
106 :     bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
107 :     REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
108 :     for (j = 0; j < nrhs; j++)
109 :     lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
110 :     }
111 : mmaechler 2770 RETURN(ans);
112 : bates 1251 }
113 : bates 1262
114 : dmbates 2210 SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)
115 : bates 1262 {
116 : dmbates 2217 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
117 : mmaechler 2227 CSP A = AS_CSP(a), B = AS_CSP(b);
118 : mmaechler 2770 R_CheckStack();
119 :     if (A->m != A->n || B->n < 1 || A->n < 1 || A->n != B->m)
120 :     error(_("Dimensions of system to be solved are inconsistent"));
121 : mmaechler 3011 // *before* Calloc()ing below [memory leak]! -- FIXME: 0-extent should work
122 : mmaechler 2770
123 : dmbates 2210 int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
124 : dmbates 2217 xnz = 10 * B->p[B->n]; /* initial estimate of nnz in x */
125 : mmaechler 2889 int k, lo = uplo_P(a)[0] == 'L', pos = 0;
126 :     int *ti = Calloc(xnz, int), *xi = Calloc(2*A->n, int); /* for cs_reach */
127 :     double *tx = Calloc(xnz, double), *wrk = Calloc( A->n, double);
128 : mmaechler 2227
129 : dmbates 2210 slot_dup(ans, b, Matrix_DimSym);
130 : mmaechler 3011
131 : dmbates 2210 xp[0] = 0;
132 : dmbates 2217 for (k = 0; k < B->n; k++) {
133 :     int top = cs_spsolve (A, B, k, xi, wrk, (int *)NULL, lo);
134 : mmaechler 2889 int nz = A->n - top;
135 : dmbates 2217
136 :     xp[k + 1] = nz + xp[k];
137 :     if (xp[k + 1] > xnz) {
138 :     while (xp[k + 1] > xnz) xnz *= 2;
139 : dmbates 2210 ti = Realloc(ti, xnz, int);
140 :     tx = Realloc(tx, xnz, double);
141 : bates 1262 }
142 : dmbates 2217 if (lo) /* increasing row order */
143 : mmaechler 2889 for(int p = top; p < A->n; p++, pos++) {
144 : dmbates 2210 ti[pos] = xi[p];
145 :     tx[pos] = wrk[xi[p]];
146 :     }
147 : dmbates 2217 else /* decreasing order, reverse copy */
148 : mmaechler 2889 for(int p = A->n - 1; p >= top; p--, pos++) {
149 : dmbates 2210 ti[pos] = xi[p];
150 :     tx[pos] = wrk[xi[p]];
151 :     }
152 : bates 1262 }
153 : dmbates 2220 xnz = xp[B->n];
154 : dmbates 2217 Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, xnz)), ti, xnz);
155 :     Memcpy( REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, xnz)), tx, xnz);
156 : mmaechler 2227
157 : mmaechler 2889 Free(ti); Free(tx);
158 : mmaechler 2770 Free(wrk); Free(xi);
159 :    
160 : mmaechler 3011 // dimnames:
161 :     SEXP dn = PROTECT(allocVector(VECSXP, 2));
162 :     SET_VECTOR_ELT(dn, 0, duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1)));
163 :     SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));
164 :     SET_SLOT(ans, Matrix_DimNamesSym, dn);
165 :     UNPROTECT(1);
166 :    
167 : mmaechler 2770 RETURN(ans);
168 : bates 1262 }
169 : mmaechler 2770 #undef RETURN
170 :    

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