Matrix r4655
Loading...
Searching...
No Matches
idz.c
Go to the documentation of this file.
1#include "Mdefines.h"
2#include "idz.h"
3
4#define IDZ \
5TEMPLATE(i, int, 0 , 1 ) \
6TEMPLATE(d, double, 0.0, 1.0) \
7TEMPLATE(z, Rcomplex, Matrix_zzero, Matrix_zone)
8
9#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
10static void _PREFIX_ ## \
11swap(int n, _CTYPE_ *x, int incx, _CTYPE_ *y, int incy) \
12{ \
13 _CTYPE_ tmp; \
14 while (n--) { \
15 tmp = *x; \
16 *x = *y; \
17 *y = tmp; \
18 x += incx; \
19 y += incy; \
20 } \
21 return; \
22}
23IDZ
24#undef TEMPLATE
25
26#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
27static void _PREFIX_ ## \
28syswapr(char uplo, int n, _CTYPE_ *x, int k0, int k1) \
29{ \
30 _CTYPE_ *x0 = x + (R_xlen_t) k0 * n, *x1 = x + (R_xlen_t) k1 * n, \
31 tmp; \
32 if (uplo == 'U') { \
33 _PREFIX_ ## swap(k0, x0, 1, x1, 1); \
34 tmp = x0[k0]; \
35 x0[k0] = x1[k1]; \
36 x1[k1] = tmp; \
37 _PREFIX_ ## swap(k1 - k0 - 1, x0 + k0 + n, n, x1 + k0 + 1, 1); \
38 _PREFIX_ ## swap(n - k1 - 1, x1 + k0 + n, n, x1 + k1 + n, n); \
39 } else { \
40 _PREFIX_ ## swap(k0, x + k0, n, x + k1, n); \
41 tmp = x0[k0]; \
42 x0[k0] = x1[k1]; \
43 x1[k1] = tmp; \
44 _PREFIX_ ## swap(k1 - k0 - 1, x0 + k0 + 1, 1, x0 + k1 + n, n); \
45 _PREFIX_ ## swap(n - k1 - 1, x0 + k1 + 1, 1, x1 + k1 + 1, 1); \
46 } \
47 return; \
48}
49IDZ
50#undef TEMPLATE
51
52#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
53void _PREFIX_ ## \
54rowperm2(_CTYPE_ *x, int m, int n, int *p, int off, int invert) \
55{ \
56 int i, k0, k1; \
57 for (i = 0; i < m; ++i) \
58 p[i] = -(p[i] - off + 1); \
59 if (!invert) { \
60 for (i = 0; i < m; ++i) { \
61 if (p[i] > 0) \
62 continue; \
63 k0 = i; \
64 p[k0] = -p[k0]; \
65 k1 = p[k0] - 1; \
66 while (p[k1] < 0) { \
67 _PREFIX_ ## swap(n, x + k0, m, x + k1, m); \
68 k0 = k1; \
69 p[k0] = -p[k0]; \
70 k1 = p[k0] - 1; \
71 } \
72 } \
73 } else { \
74 for (i = 0; i < m; ++i) { \
75 if (p[i] > 0) \
76 continue; \
77 k0 = i; \
78 p[k0] = -p[k0]; \
79 k1 = p[k0] - 1; \
80 while (k1 != k0) { \
81 _PREFIX_ ## swap(n, x + k0, m, x + k1, m); \
82 p[k1] = -p[k1]; \
83 k1 = p[k1] - 1; \
84 } \
85 } \
86 } \
87 for (i = 0; i < m; ++i) \
88 p[i] = p[i] + off - 1; \
89 return; \
90}
91IDZ
92#undef TEMPLATE
93
94#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
95void _PREFIX_ ## \
96symperm2(_CTYPE_ *x, int n, char uplo, int *p, int off, int invert) \
97{ \
98 int i, k0, k1; \
99 for (i = 0; i < n; ++i) \
100 p[i] = -(p[i] - off + 1); \
101 if (!invert) { \
102 for (i = 0; i < n; ++i) { \
103 if (p[i] > 0) \
104 continue; \
105 k0 = i; \
106 p[k0] = -p[k0]; \
107 k1 = p[k0] - 1; \
108 while (p[k1] < 0) { \
109 if (k0 < k1) \
110 _PREFIX_ ## syswapr(uplo, n, x, k0, k1); \
111 else \
112 _PREFIX_ ## syswapr(uplo, n, x, k1, k0); \
113 k0 = k1; \
114 p[k0] = -p[k0]; \
115 k1 = p[k0] - 1; \
116 } \
117 } \
118 } else { \
119 for (i = 0; i < n; ++i) { \
120 if (p[i] > 0) \
121 continue; \
122 k0 = i; \
123 p[k0] = -p[k0]; \
124 k1 = p[k0] - 1; \
125 while (k1 != k0) { \
126 if (k0 < k1) \
127 _PREFIX_ ## syswapr(uplo, n, x, k0, k1); \
128 else \
129 _PREFIX_ ## syswapr(uplo, n, x, k1, k0); \
130 p[k1] = -p[k1]; \
131 k1 = p[k1] - 1; \
132 } \
133 } \
134 } \
135 for (i = 0; i < n; ++i) \
136 p[i] = p[i] + off - 1; \
137 return; \
138}
139IDZ
140#undef TEMPLATE
141
142#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
143void _PREFIX_ ## \
144pack2(_CTYPE_ *dest, const _CTYPE_ *src, int n, char uplo, char diag) \
145{ \
146 int i, j; \
147 R_xlen_t dpos = 0, spos = 0; \
148 if (uplo == 'U') { \
149 for (j = 0; j < n; spos += n-(++j)) \
150 for (i = 0; i <= j; ++i) \
151 dest[dpos++] = src[spos++]; \
152 if (diag != 'N') { \
153 dpos = 0; \
154 for (j = 0; j < n; dpos += (++j)+1) \
155 dest[dpos] = _ONE_; \
156 } \
157 } else { \
158 for (j = 0; j < n; spos += (++j)) \
159 for (i = j; i < n; ++i) \
160 dest[dpos++] = src[spos++]; \
161 if (diag != 'N') { \
162 dpos = 0; \
163 for (j = 0; j < n; dpos += n-(j++)) \
164 dest[dpos] = _ONE_; \
165 } \
166 } \
167 return; \
168}
169IDZ
170#undef TEMPLATE
171
172#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
173void _PREFIX_ ## \
174unpack1(_CTYPE_ *dest, const _CTYPE_ *src, int n, char uplo, char diag) \
175{ \
176 int i, j; \
177 R_xlen_t dpos = 0, spos = 0; \
178 if (uplo == 'U') { \
179 for (j = 0; j < n; dpos += n-(++j)) \
180 for (i = 0; i <= j; ++i) \
181 dest[dpos++] = src[spos++]; \
182 } else { \
183 for (j = 0; j < n; dpos += (++j)) \
184 for (i = j; i < n; ++i) \
185 dest[dpos++] = src[spos++]; \
186 } \
187 if (diag != 'N') { \
188 dpos = 0; \
189 R_xlen_t n1a = (R_xlen_t) n + 1; \
190 for (j = 0; j < n; ++j, dpos += n1a) \
191 dest[dpos] = _ONE_; \
192 } \
193 return; \
194}
195IDZ
196#undef TEMPLATE
197
198#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
199void _PREFIX_ ## \
200transpose2(_CTYPE_ *dest, const _CTYPE_ *src, int m, int n) \
201{ \
202 R_xlen_t mn1s = (R_xlen_t) m * n - 1; \
203 int i, j; \
204 for (j = 0; j < m; ++j, src -= mn1s) \
205 for (i = 0; i < n; ++i, src += m) \
206 *(dest++) = *src; \
207 return; \
208}
209IDZ
210#undef TEMPLATE
211
212#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
213void _PREFIX_ ## \
214transpose1(_CTYPE_ *dest, const _CTYPE_ *src, int n, char uplo) \
215{ \
216 int i, j; \
217 if (uplo == 'U') { \
218 for (j = 0; j < n; ++j) \
219 for (i = j; i < n; ++i) \
220 *(dest++) = *(src + PACKED_AR21_UP(j, i)); \
221 } else { \
222 R_xlen_t n2 = (R_xlen_t) n * 2; \
223 for (j = 0; j < n; ++j) \
224 for (i = 0; i <= j; ++i) \
225 *(dest++) = *(src + PACKED_AR21_LO(j, i, n2)); \
226 } \
227 return; \
228}
229IDZ
230#undef TEMPLATE
231
232#define ASSIGN_JJ_i(_X_)
233#define ASSIGN_JJ_d(_X_)
234#define ASSIGN_JJ_z(_X_) \
235 _X_.i = 0.0
236#define ASSIGN_JI_i(_X_, _Y_) \
237 _X_ = _Y_
238#define ASSIGN_JI_d(_X_, _Y_) \
239 _X_ = _Y_
240#define ASSIGN_JI_z(_X_, _Y_) \
241 do { \
242 _X_.r = _Y_.r; \
243 _X_.i = -_Y_.i; \
244 } while (0)
245
246#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
247void _PREFIX_ ## \
248syforce2(_CTYPE_ *x, int n, char uplo) \
249{ \
250 _CTYPE_ *y = x; \
251 int i, j; \
252 if (uplo == 'U') { \
253 for (j = 0; j < n; ++j) { \
254 ASSIGN_JJ_ ## _PREFIX_((*x)); \
255 x += 1; \
256 y += n; \
257 for (i = j + 1; i < n; ++i) { \
258 ASSIGN_JI_ ## _PREFIX_((*x), (*y)); \
259 x += 1; \
260 y += n; \
261 } \
262 x = y = x + j + 1; \
263 } \
264 } else { \
265 for (j = 0; j < n; ++j) { \
266 ASSIGN_JJ_ ## _PREFIX_((*y)); \
267 x += 1; \
268 y += n; \
269 for (i = j + 1; i < n; ++i) { \
270 ASSIGN_JI_ ## _PREFIX_((*y), (*x)); \
271 x += 1; \
272 y += n; \
273 } \
274 x = y = x + j + 1; \
275 } \
276 } \
277 return; \
278}
279IDZ
280#undef TEMPLATE
281
282#undef ASSIGN_JJ_i
283#undef ASSIGN_JJ_d
284#undef ASSIGN_JJ_z
285#undef ASSIGN_JI_i
286#undef ASSIGN_JI_d
287#undef ASSIGN_JI_z
288
289#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
290void _PREFIX_ ## \
291trforce2(_CTYPE_ *x, int m, int n, char uplo, char diag) \
292{ \
293 _CTYPE_ *y = x; \
294 int i, j, r = (m < n) ? m : n; \
295 if (uplo == 'U') { \
296 for (j = 0; j < r; ++j) { \
297 for (i = j + 1; i < m; ++i) \
298 *(++x) = _ZERO_; \
299 x += j + 2; \
300 } \
301 } else { \
302 for (j = 0; j < r; ++j) { \
303 for (i = 0; i < j; ++i) \
304 *(x++) = _ZERO_; \
305 x += m - j; \
306 } \
307 for (j = r; j < n; ++j) \
308 for (i = 0; i < m; ++i) \
309 *(x++) = _ZERO_; \
310 } \
311 if (diag != 'N') { \
312 R_xlen_t m1a = (R_xlen_t) m + 1; \
313 for (j = 0; j < r; ++j) { \
314 *y = _ONE_; \
315 y += m1a; \
316 } \
317 } \
318 return; \
319}
320IDZ
321#undef TEMPLATE
322
323#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
324void _PREFIX_ ## \
325band2(_CTYPE_ *x, int m, int n, int a, int b, char diag) \
326{ \
327 if (m == 0 || n == 0) \
328 return; \
329 if (a > b || a >= n || b <= -m) { \
330 Matrix_memset(x, 0, (R_xlen_t) m * n, sizeof(_CTYPE_)); \
331 return; \
332 } \
333 if (a <= -m) a = 1-m; \
334 if (b >= n) b = n-1; \
335 \
336 int i, j, i0, i1, \
337 j0 = (a < 0) ? 0 : a, \
338 j1 = (b < n-m) ? m+b : n; \
339 \
340 if (j0 > 0) { \
341 R_xlen_t dx = (R_xlen_t) m * j0; \
342 Matrix_memset(x, 0, dx, sizeof(_CTYPE_)); \
343 x += dx; \
344 } \
345 for (j = j0; j < j1; ++j, x += m) { \
346 i0 = j - b; \
347 i1 = j - a + 1; \
348 for (i = 0; i < i0; ++i) \
349 *(x + i) = _ZERO_; \
350 for (i = i1; i < m; ++i) \
351 *(x + i) = _ZERO_; \
352 } \
353 if (j1 < n) \
354 Matrix_memset(x, 0, (R_xlen_t) m * (n - j1), sizeof(_CTYPE_)); \
355 if (diag != 'N' && a <= 0 && b >= 0) { \
356 x -= m * (R_xlen_t) j; \
357 R_xlen_t m1a = (R_xlen_t) m + 1; \
358 for (j = 0; j < n; ++j, x += m1a) \
359 *x = _ONE_; \
360 } \
361 return; \
362}
363IDZ
364#undef TEMPLATE
365
366#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
367void _PREFIX_ ## \
368band1(_CTYPE_ *x, int n, int a, int b, char uplo, char diag) \
369{ \
370 if (n == 0) \
371 return; \
372 if (a > b || a >= n || b <= -n) { \
373 Matrix_memset(x, 0, PACKED_LENGTH(n), sizeof(_CTYPE_)); \
374 return; \
375 } \
376 if (uplo == 'U') { \
377 if (a < 0) a = 0; \
378 if (b >= n) b = n-1; \
379 } else { \
380 if (b > 0) b = 0; \
381 if (a <= -n) a = 1-n; \
382 } \
383 \
384 int i, j, i0, i1, \
385 j0 = (a < 0) ? 0 : a, \
386 j1 = (b < 0) ? n+b : n; \
387 \
388 if (uplo == 'U') { \
389 if (j0 > 0) { \
390 R_xlen_t dx; \
391 Matrix_memset(x, 0, dx = PACKED_LENGTH(j0), \
392 sizeof(_CTYPE_)); \
393 x += dx; \
394 } \
395 for (j = j0; j < j1; x += (++j)) { \
396 i0 = j - b; \
397 i1 = j - a + 1; \
398 for (i = 0; i < i0; ++i) \
399 *(x + i) = _ZERO_; \
400 for (i = i1; i <= j; ++i) \
401 *(x + i) = _ZERO_; \
402 } \
403 if (j1 < n) \
404 Matrix_memset(x, 0, PACKED_LENGTH(n) - PACKED_LENGTH(j1), \
405 sizeof(_CTYPE_)); \
406 if (diag != 'N' && a == 0) { \
407 x -= PACKED_LENGTH(j); \
408 for (j = 0; j < n; x += (++j)+1) \
409 *x = _ONE_; \
410 } \
411 } else { \
412 if (j0 > 0) { \
413 R_xlen_t dx = PACKED_LENGTH(n) - PACKED_LENGTH(j0); \
414 Matrix_memset(x, 0, dx, sizeof(_CTYPE_)); \
415 x += dx; \
416 } \
417 for (j = j0; j < j1; x += n-(j++)) { \
418 i0 = j - b; \
419 i1 = j - a + 1; \
420 for (i = j; i < i0; ++i) \
421 *(x + i - j) = _ZERO_; \
422 for (i = i1; i < n; ++i) \
423 *(x + i - j) = _ZERO_; \
424 } \
425 if (j1 < n) \
426 Matrix_memset(x, 0, PACKED_LENGTH(n - j1), \
427 sizeof(_CTYPE_)); \
428 if (diag != 'N' && b == 0) { \
429 x -= PACKED_LENGTH(n) - PACKED_LENGTH(j); \
430 for (j = 0; j < n; x += n-(j++)) \
431 *x = _ONE_; \
432 } \
433 } \
434 return; \
435}
436IDZ
437#undef TEMPLATE
438
439#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
440void _PREFIX_ ## \
441dcpy2(_CTYPE_ *dest, const _CTYPE_ *src, int n, R_xlen_t length, \
442 char uplo, char diag) \
443{ \
444 int j; \
445 R_xlen_t n1a = (R_xlen_t) n + 1; \
446 if (diag != 'N') { \
447 for (j = 0; j < n; ++j, dest += n1a) \
448 *dest = _ONE_; \
449 } else if (length == n) { \
450 /* copying from diagonalMatrix */ \
451 for (j = 0; j < n; ++j, dest += n1a, ++src) \
452 *dest = *src; \
453 } else if (length == (n * n1a) / 2) { \
454 /* copying from packedMatrix */ \
455 if (uplo == 'U') { \
456 for (j = 0; j < n; dest += n1a, src += (++j)+1) \
457 *dest = *src; \
458 } else { \
459 for (j = 0; j < n; dest += n1a, src += n-(j++)) \
460 *dest = *src; \
461 } \
462 } else if (length == (R_xlen_t) n * n) { \
463 /* copying from square unpackedMatrix */ \
464 for (j = 0; j < n; ++j, dest += n1a, src += n1a) \
465 *dest = *src; \
466 } else { \
467 error(_("incompatible '%s' and '%s' in '%s'"), \
468 "n", "length", __func__); \
469 } \
470 return; \
471}
472IDZ
473#undef TEMPLATE
474
475#define TEMPLATE(_PREFIX_, _CTYPE_, _ZERO_, _ONE_) \
476void _PREFIX_ ## \
477dcpy1(_CTYPE_ *dest, const _CTYPE_ *src, int n, R_xlen_t length, \
478 char uplo_dest, char uplo_src, char diag) \
479{ \
480 int j; \
481 if (diag != 'N') { \
482 if (uplo_dest == 'U') { \
483 for (j = 0; j < n; dest += (++j)+1) \
484 *dest = _ONE_; \
485 } else { \
486 for (j = 0; j < n; dest += n-(j++)) \
487 *dest = _ONE_; \
488 } \
489 } else if (length == n) { \
490 /* copying from diagonalMatrix */ \
491 if (uplo_dest == 'U') { \
492 for (j = 0; j < n; dest += (++j)+1, ++src) \
493 *dest = *src; \
494 } else { \
495 for (j = 0; j < n; dest += n-(j++), ++src) \
496 *dest = *src; \
497 } \
498 } else if (length == PACKED_LENGTH(n)) { \
499 /* copying from packedMatrix */ \
500 if (uplo_dest == 'U') { \
501 if (uplo_src == 'U') { \
502 for (j = 0; j < n; src += (++j)+1, dest += j+1) \
503 *dest = *src; \
504 } else { \
505 for (j = 0; j < n; src += n-j, dest += (++j)+1) \
506 *dest = *src; \
507 } \
508 } else { \
509 if (uplo_src == 'U') { \
510 for (j = 0; j < n; dest += n-(j++), src += j+1) \
511 *dest = *src; \
512 } else { \
513 for (j = 0; j < n; dest += n-j, src += n-(j++)) \
514 *dest = *src; \
515 } \
516 } \
517 } else if (length == (R_xlen_t) n * n) { \
518 /* copying from square unpackedMatrix */ \
519 R_xlen_t n1a = (R_xlen_t) n + 1; \
520 if (uplo_dest == 'U') { \
521 for (j = 0; j < n; dest += (++j)+1, src += n1a) \
522 *dest = *src; \
523 } else { \
524 for (j = 0; j < n; dest += n-(j++), src += n1a) \
525 *dest = *src; \
526 } \
527 } else { \
528 error(_("incompatible '%s' and '%s' in '%s'"), \
529 "n", "length", __func__); \
530 } \
531 return; \
532}
533IDZ
534#undef TEMPLATE
535
536#undef IDZ
IDZ
#define IDZ
Definition idz.c:4
Generated by doxygen 1.9.8