Matrix r5059
Loading...
Searching...
No Matches
band.c
Go to the documentation of this file.
1/* C implementation of methods for band */
2
3#include "Mdefines.h"
4#include "M5.h"
5#include "idz.h"
6
7SEXP dense_band(SEXP from, const char *class, int a, int b)
8{
9 int packed = class[2] == 'p';
10
11 /* Need tri[ul](<0-by-0>) and tri[ul](<1-by-1>) */
12 /* to be triangularMatrix */
13 int *pdim = DIM(from), m = pdim[0], n = pdim[1];
14 if ((m == 0 || n == 0 || (a <= 1 - m && b >= n - 1)) &&
15 (m != n || n > 1 || class[1] == 't'))
16 return from;
17
18 int ge, sy, tr;
19 tr = class[1] == 't' || (m == n && (a >= 0 || b <= 0));
20 sy = !tr && class[1] == 's' && a == -b;
21 ge = !tr && !sy;
22
23 char ul0 = '\0', ct0 = '\0', nu0 = '\0';
24 if (class[1] != 'g')
25 ul0 = UPLO(from);
26 if (class[1] == 's' && class[0] == 'z')
27 ct0 = TRANS(from);
28 if (class[1] == 't') {
29 /* Be fast if band contains entire triangle */
30 if ((ul0 == 'U') ? (a <= 0 && b >= n - 1) : (a <= 1 - m && b >= 0))
31 return from;
32 nu0 = DIAG(from);
33 }
34
35 char cl[] = "...Matrix";
36 cl[0] = class[0];
37 cl[1] = (ge) ? 'g' : ((sy) ? 's' : 't') ;
38 cl[2] = (ge) ? 'e' : ((packed) ? 'p' : ((sy) ? 'y' : 'r'));
39 SEXP to = PROTECT(newObject(cl));
40
41 SET_DIM(to, m, n);
42 SET_DIMNAMES(to, -(class[1] == 's' && !sy), DIMNAMES(from, 0));
43
44 char ul1 = (tr && class[1] != 't') ? ((a >= 0) ? 'U' : 'L') : ul0;
45 if (ul1 != '\0' && ul1 != 'U')
46 SET_UPLO(to);
47 if (ct0 != '\0' && ct0 != 'C' && sy)
48 SET_TRANS(to);
49 if (nu0 != '\0' && nu0 != 'N' && tr && a <= 0 && b >= 0)
50 SET_DIAG(to);
51
52 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
53 x1 = PROTECT(Rf_allocVector(TYPEOF(x0), (!packed || ge) ? (R_xlen_t) m * n : (R_xlen_t) PACKED_LENGTH((size_t) n)));
54
55 size_t
56 m_ = (size_t) m,
57 n_ = (size_t) n,
58 a_ = (size_t) ((int_fast64_t) m + a),
59 b_ = (size_t) ((int_fast64_t) m + b);
60
61#define TEMPLATE(c) \
62 do { \
63 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
64 if (ge && class[1] != 'g') { \
65 if (!packed) \
66 c##NAME(force2)(px1, px0, n_, ul0, ct0, nu0); \
67 else \
68 c##NAME( pack1)(px1, px0, n_, ul0, ct0, nu0); \
69 px0 = NULL; \
70 packed = 0; \
71 } else if (tr && class[1] == 's' && ul0 != ul1) { \
72 if (!packed) \
73 c##NAME(trans2)(px1, px0, m_, n_ , ct0); \
74 else \
75 c##NAME(trans1)(px1, px0, n_, ul0, ct0); \
76 px0 = NULL; \
77 ul0 = ul1; \
78 } \
79 if (!packed) \
80 c##NAME( band2)(px1, px0, m_, n_ , a_, b_); \
81 else \
82 c##NAME( band1)(px1, px0, n_, ul0, a_, b_); \
83 } while (0)
84
85 SWITCH4(class[0], TEMPLATE);
86
87#undef TEMPLATE
88
89 SET_SLOT(to, Matrix_xSym, x1);
90
91 UNPROTECT(3); /* x1, x0, to */
92 return to;
93}
94
95SEXP sparse_band(SEXP from, const char *class, int a, int b)
96{
97 /* Need tri[ul](<0-by-0>) and tri[ul](<1-by-1>) */
98 /* to be triangularMatrix */
99 int *pdim = DIM(from), m = pdim[0], n = pdim[1];
100 if ((m == 0 || n == 0 || (a <= 1 - m && b >= n - 1)) &&
101 (m != n || n > 1 || class[1] == 't'))
102 return from;
103
104 int ge, sy, tr;
105 tr = class[1] == 't' || (m == n && (a >= 0 || b <= 0));
106 sy = !tr && class[1] == 's' && a == -b;
107 ge = !tr && !sy;
108
109 char ul0 = '\0', ct0 = '\0', nu0 = '\0';
110 if (class[1] != 'g')
111 ul0 = UPLO(from);
112 if (class[1] == 's' && class[0] == 'z')
113 ct0 = TRANS(from);
114 if (class[1] == 't') {
115 /* Be fast if band contains entire triangle */
116 if ((ul0 == 'U') ? (a <= 0 && b >= n - 1) : (a <= 1 - m && b >= 0))
117 return from;
118 nu0 = DIAG(from);
119 }
120
121 char cl[] = "...Matrix";
122 cl[0] = class[0];
123 cl[1] = (ge) ? 'g' : ((sy) ? 's' : 't');
124 cl[2] = class[2];
125 SEXP to = PROTECT(newObject(cl));
126
127 SET_DIM(to, m, n);
128 SET_DIMNAMES(to, -(class[1] == 's' && !sy), DIMNAMES(from, 0));
129
130 char ul1 = (tr && class[1] != 't') ? ((a >= 0) ? 'U' : 'L') : ul0;
131 if (ul1 != '\0' && ul1 != 'U')
132 SET_UPLO(to);
133 if (ct0 != '\0' && ct0 != 'C' && sy)
134 SET_TRANS(to);
135 if (nu0 != '\0' && nu0 != 'N' && tr && a <= 0 && b >= 0)
136 SET_DIAG(to);
137
138 if (class[2] != 'T') {
139
140 if (class[2] == 'R') {
141 SWAP(m, n, int, );
142 SWAP(a, b, int, -);
143 }
144
145 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
146 p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
147 i0 = PROTECT(GET_SLOT(from, iSym)),
148 p1 = PROTECT(Rf_allocVector(INTSXP, XLENGTH(p0)));
149 int *pp0 = INTEGER(p0), *pi0 = INTEGER(i0),
150 *pp1 = INTEGER(p1), *iwork = NULL,
151 j, k, kend, nnz0 = pp0[n], nnz1 = 0, d;
152 pp0++; *(pp1++) = 0;
153 SET_SLOT(to, Matrix_pSym, p1);
154
155 if (class[1] != 's' || sy) {
156 for (j = 0, k = 0; j < n; ++j) {
157 kend = pp0[j];
158 while (k < kend) {
159 if ((d = j - pi0[k]) >= a && d <= b)
160 ++nnz1;
161 ++k;
162 }
163 pp1[j] = nnz1;
164 }
165 if (nnz0 == nnz1) {
166 SET_SLOT(to, iSym, i0);
167 if (class[0] != 'n') {
168 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
169 SET_SLOT(to, Matrix_xSym, x0);
170 UNPROTECT(1); /* x0 */
171 }
172 UNPROTECT(4); /* p1, i0, p0, to */
173 return to;
174 }
175 } else {
176 memset(pp1, 0, sizeof(int) * (size_t) n);
177 for (j = 0, k = 0; j < n; ++j) {
178 kend = pp0[j];
179 while (k < kend) {
180 if ((d = j - pi0[k]) >= a && d <= b)
181 ++pp1[j];
182 if (d != 0 && -d >= a && -d <= b)
183 ++pp1[pi0[k]];
184 ++k;
185 }
186 }
187 Matrix_Calloc(iwork, n, int);
188 for (j = 0; j < n; ++j) {
189 iwork[j] = nnz1;
190 nnz1 += pp1[j];
191 pp1[j] = nnz1;
192 }
193 }
194
195 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, nnz1));
196 int *pi1 = INTEGER(i1);
197 SET_SLOT(to, iSym, i1);
198
199#define TEMPLATE(c) \
200 do { \
201 c##IF_NPATTERN( \
202 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)), \
203 x1 = PROTECT(Rf_allocVector(c##TYPESXP, nnz1)); \
204 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
205 SET_SLOT(to, Matrix_xSym, x1); \
206 UNPROTECT(2); /* x1, x0 */ \
207 ); \
208 if (class[1] != 's' || sy) \
209 for (j = 0, k = 0; j < n; ++j) { \
210 kend = pp0[j]; \
211 while (k < kend) { \
212 if ((d = j - pi0[k]) >= a && d <= b) { \
213 *(pi1++) = pi0[k]; \
214 c##IF_NPATTERN( \
215 *(px1++) = px0[k]; \
216 ); \
217 } \
218 ++k; \
219 } \
220 } \
221 else { \
222 for (j = 0, k = 0; j < n; ++j) { \
223 kend = pp0[j]; \
224 while (k < kend) { \
225 if ((d = j - pi0[k]) >= a && d <= b) { \
226 c##IF_NPATTERN( \
227 if (ct0 == 'C' && d == 0) \
228 c##ASSIGN_PROJ_REAL(px1[iwork[j]], px0[k]); \
229 else \
230 c##ASSIGN_IDEN (px1[iwork[j]], px0[k]); \
231 ); \
232 pi1[iwork[j]++] = pi0[k]; \
233 } \
234 if (d != 0 && -d >= a && -d <= b) { \
235 c##IF_NPATTERN( \
236 if (ct0 == 'C') \
237 c##ASSIGN_CONJ(px1[iwork[pi0[k]]], px0[k]); \
238 else \
239 c##ASSIGN_IDEN(px1[iwork[pi0[k]]], px0[k]); \
240 ); \
241 pi1[iwork[pi0[k]]++] = j; \
242 } \
243 ++k; \
244 } \
245 } \
246 Matrix_Free(iwork, n); \
247 } \
248 } while (0)
249
250 SWITCH5(class[0], TEMPLATE);
251
252#undef TEMPLATE
253
254 UNPROTECT(4); /* i1, p1, i0, p0 */
255
256 } else {
257
258 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
259 j0 = PROTECT(GET_SLOT(from, Matrix_jSym));
260 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0), d;
261 R_xlen_t k, nnz0 = XLENGTH(i0), nnz1 = 0;
262
263 if (class[1] != 's' || sy) {
264 for (k = 0; k < nnz0; ++k)
265 if ((d = pj0[k] - pi0[k]) >= a && d <= b)
266 ++nnz1;
267 if (nnz0 == nnz1) {
268 SET_SLOT(to, Matrix_iSym, i0);
269 SET_SLOT(to, Matrix_jSym, j0);
270 if (class[0] != 'n') {
271 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
272 SET_SLOT(to, Matrix_xSym, x0);
273 UNPROTECT(1); /* x0 */
274 }
275 UNPROTECT(3); /* j0, i0, to */
276 return to;
277 }
278 } else
279 for (k = 0; k < nnz0; ++k) {
280 if ((d = pj0[k] - pi0[k]) >= a && d <= b)
281 ++nnz1;
282 if (d != 0 && -d >= a && -d <= b)
283 ++nnz1;
284 }
285
286 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
287 j1 = PROTECT(Rf_allocVector(INTSXP, nnz1));
288 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
289 SET_SLOT(to, Matrix_iSym, i1);
290 SET_SLOT(to, Matrix_jSym, j1);
291
292#define TEMPLATE(c) \
293 do { \
294 c##IF_NPATTERN( \
295 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)), \
296 x1 = PROTECT(Rf_allocVector(c##TYPESXP, nnz1)); \
297 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
298 SET_SLOT(to, Matrix_xSym, x1); \
299 UNPROTECT(2); /* x1, x0 */ \
300 ); \
301 if (class[1] != 's' || sy) \
302 for (k = 0; k < nnz0; ++k) { \
303 if ((d = pj0[k] - pi0[k]) >= a && d <= b) { \
304 *(pi1++) = pi0[k]; \
305 *(pj1++) = pj0[k]; \
306 c##IF_NPATTERN( \
307 *(px1++) = px0[k]; \
308 ); \
309 } \
310 } \
311 else \
312 for (k = 0; k < nnz0; ++k) { \
313 if ((d = pj0[k] - pi0[k]) >= a && d <= b) { \
314 *(pi1++) = pi0[k]; \
315 *(pj1++) = pj0[k]; \
316 c##IF_NPATTERN( \
317 if (ct0 == 'C' && d == 0) \
318 c##ASSIGN_PROJ_REAL(*px1, px0[k]); \
319 else \
320 c##ASSIGN_IDEN (*px1, px0[k]); \
321 px1++; \
322 ); \
323 } \
324 if (d != 0 && -d >= a && -d <= b) { \
325 *(pi1++) = pj0[k]; \
326 *(pj1++) = pi0[k]; \
327 c##IF_NPATTERN( \
328 if (ct0 == 'C') \
329 c##ASSIGN_CONJ(*px1, px0[k]); \
330 else \
331 c##ASSIGN_IDEN(*px1, px0[k]); \
332 px1++; \
333 ); \
334 } \
335 } \
336 } while (0)
337
338 SWITCH5(class[0], TEMPLATE);
339
340#undef TEMPLATE
341
342 UNPROTECT(4); /* j1, i1, j0, i0 */
343
344 }
345
346 UNPROTECT(1); /* to */
347 return to;
348}
349
350SEXP R_dense_band(SEXP s_from, SEXP s_a, SEXP s_b)
351{
352 if (TYPEOF(s_from) != OBJSXP) {
353 /* defined in ./coerce.c : */
354 SEXP matrix_as_dense(SEXP, const char *, char, char, char, int, int);
355 s_from = matrix_as_dense(s_from, ".ge", '\0', '\0', '\0', 1, 0);
356 }
357 PROTECT(s_from);
358 const char *class = Matrix_class(s_from, valid_dense, 6, __func__);
359
360 int *pdim = DIM(s_from), m = pdim[0], n = pdim[1], a, b;
361 if (s_a == R_NilValue)
362 a = -m;
363 else if ((a = Rf_asInteger(s_a)) == NA_INTEGER || a < -m || a > n)
364 Rf_error(_("'%s' (%d) must be an integer from %s (%d) to %s (%d)"),
365 "k1", a, "-Dim[1]", -m, "Dim[2]", n);
366 if (s_b == R_NilValue)
367 b = n;
368 else if ((b = Rf_asInteger(s_b)) == NA_INTEGER || b < -m || b > n)
369 Rf_error(_("'%s' (%d) must be an integer from %s (%d) to %s (%d)"),
370 "k2", b, "-Dim[1]", -m, "Dim[2]", n);
371 else if (b < a)
372 Rf_error(_("'%s' (%d) must be less than or equal to '%s' (%d)"),
373 "k1", a, "k2", b);
374
375 s_from = dense_band(s_from, class, a, b);
376 UNPROTECT(1);
377 return s_from;
378}
379
380SEXP R_sparse_band(SEXP s_from, SEXP s_a, SEXP s_b)
381{
382 const char *class = Matrix_class(s_from, valid_sparse, 6, __func__);
383
384 int *pdim = DIM(s_from), m = pdim[0], n = pdim[1], a, b;
385 if (s_a == R_NilValue)
386 a = -m;
387 else if ((a = Rf_asInteger(s_a)) == NA_INTEGER || a < -m || a > n)
388 Rf_error(_("'%s' (%d) must be an integer from %s (%d) to %s (%d)"),
389 "k1", a, "-Dim[1]", -m, "Dim[2]", n);
390 if (s_b == R_NilValue)
391 b = n;
392 else if ((b = Rf_asInteger(s_b)) == NA_INTEGER || b < -m || b > n)
393 Rf_error(_("'%s' (%d) must be an integer from %s (%d) to %s (%d)"),
394 "k2", b, "-Dim[1]", -m, "Dim[2]", n);
395 else if (b < a)
396 Rf_error(_("'%s' (%d) must be less than or equal to '%s' (%d)"),
397 "k1", a, "k2", b);
398
399 return sparse_band(s_from, class, a, b);
400}
SWITCH5
#define SWITCH5(c, template)
Definition M5.h:327
SWITCH4
#define SWITCH4(c, template)
Definition M5.h:315
SET_DIM
#define SET_DIM(x, m, n)
Definition Mdefines.h:87
SWAP
#define SWAP(a, b, t, op)
Definition Mdefines.h:138
valid_dense
const char * valid_dense[]
Definition objects.c:3
Matrix_Calloc
#define Matrix_Calloc(p, n, t)
Definition Mdefines.h:45
_
#define _(String)
Definition Mdefines.h:66
SET_UPLO
#define SET_UPLO(x)
Definition Mdefines.h:103
DIAG
#define DIAG(x)
Definition Mdefines.h:111
UPLO
#define UPLO(x)
Definition Mdefines.h:101
SET_DIMNAMES
#define SET_DIMNAMES(x, mode, value)
Definition Mdefines.h:98
Matrix_class
const char * Matrix_class(SEXP, const char **, int, const char *)
Definition objects.c:112
DIMNAMES
#define DIMNAMES(x, mode)
Definition Mdefines.h:96
SET_TRANS
#define SET_TRANS(x)
Definition Mdefines.h:108
TRANS
#define TRANS(x)
Definition Mdefines.h:106
SET_DIAG
#define SET_DIAG(x)
Definition Mdefines.h:113
DIM
#define DIM(x)
Definition Mdefines.h:85
GET_SLOT
#define GET_SLOT(x, name)
Definition Mdefines.h:72
PACKED_LENGTH
#define PACKED_LENGTH(n)
Definition Mdefines.h:132
newObject
SEXP newObject(const char *)
Definition objects.c:13
valid_sparse
const char * valid_sparse[]
Definition Mdefines.h:328
SET_SLOT
#define SET_SLOT(x, name, value)
Definition Mdefines.h:73
TYPEOF
#define TYPEOF(s)
Definition Mdefines.h:123
R_dense_band
SEXP R_dense_band(SEXP s_from, SEXP s_a, SEXP s_b)
Definition band.c:350
R_sparse_band
SEXP R_sparse_band(SEXP s_from, SEXP s_a, SEXP s_b)
Definition band.c:380
dense_band
SEXP dense_band(SEXP from, const char *class, int a, int b)
Definition band.c:7
TEMPLATE
#define TEMPLATE(c)
sparse_band
SEXP sparse_band(SEXP from, const char *class, int a, int b)
Definition band.c:95
cl
cholmod_common cl
Definition cholmod-etc.c:6
matrix_as_dense
SEXP matrix_as_dense(SEXP from, const char *zzz, char ul, char ct, char nu, int mg, int new)
Definition coerce.c:262
Matrix_xSym
SEXP Matrix_xSym
Definition init.c:635
Matrix_iSym
SEXP Matrix_iSym
Definition init.c:607
Matrix_jSym
SEXP Matrix_jSym
Definition init.c:610
Matrix_pSym
SEXP Matrix_pSym
Definition init.c:622
Generated by doxygen 1.12.0