1 /* Copyright (C) 2003-2007 CAMP
2 * Please see the accompanying LICENSE file for further information. */
3
4 #include "../extensions.h"
5 #include "bmgs.h"
6
7 #ifdef K
8
9 void
IP1D(const T * a,const int n,const int m,T * restrict b,int skip[2])10 IP1D(const T* a, const int n, const int m, T* restrict b, int skip[2])
11 {
12 a += K / 2 - 1;
13
14 for (int j = 0; j < m; j++) {
15 const T* aa = a + j * (K - 1 - skip[1] + n);
16 T* restrict bb = b + j;
17
18 for (int i = 0; i < n; i++) {
19 if (i == 0 && skip[0])
20 bb -= m;
21 else
22 bb[0] = aa[0];
23
24 if (i == n - 1 && skip[1])
25 bb -= m;
26 else {
27 if (K == 2)
28 bb[m] = 0.5 * (aa[0] + aa[1]);
29 else if (K == 4)
30 bb[m] = ( 0.5625 * (aa[ 0] + aa[1]) +
31 -0.0625 * (aa[-1] + aa[2]));
32 else if (K == 6)
33 bb[m] = ( 0.58593750 * (aa[ 0] + aa[1]) +
34 -0.09765625 * (aa[-1] + aa[2]) +
35 0.01171875 * (aa[-2] + aa[3]));
36 else
37 bb[m] = ( 0.59814453125 * (aa[ 0] + aa[1]) +
38 -0.11962890625 * (aa[-1] + aa[2]) +
39 0.02392578125 * (aa[-2] + aa[3]) +
40 -0.00244140625 * (aa[-3] + aa[4]));
41 }
42 aa++;
43 bb += 2 * m;
44 }
45 }
46 }
47
48 #else
49 # define K 2
50 # define IP1D Z(bmgs_interpolate1D2)
51 # include "interpolate.c"
52 # undef IP1D
53 # undef K
54 # define K 4
55 # define IP1D Z(bmgs_interpolate1D4)
56 # include "interpolate.c"
57 # undef IP1D
58 # undef K
59 # define K 6
60 # define IP1D Z(bmgs_interpolate1D6)
61 # include "interpolate.c"
62 # undef IP1D
63 # undef K
64 # define K 8
65 # define IP1D Z(bmgs_interpolate1D8)
66 # include "interpolate.c"
67 # undef IP1D
68 # undef K
69
70 void
Z(bmgs_interpolate)71 Z(bmgs_interpolate)(int k, int skip[3][2],
72 const T* a, const int size[3], T* restrict b, T* restrict w)
73 {
74 void (*ip)(const T*, int, int, T*, int[2]);
75 int e;
76
77 if (k == 2)
78 ip = Z(bmgs_interpolate1D2);
79 else if (k == 4)
80 ip = Z(bmgs_interpolate1D4);
81 else if (k == 6)
82 ip = Z(bmgs_interpolate1D6);
83 else
84 ip = Z(bmgs_interpolate1D8);
85
86 e = k - 1;
87
88 ip(a, size[2] - e + skip[2][1],
89 size[0] *
90 size[1],
91 b, skip[2]);
92 ip(b, size[1] - e + skip[1][1],
93 size[0] *
94 ((size[2] - e) * 2 - skip[2][0] + skip[2][1]),
95 w, skip[1]);
96 ip(w, size[0] - e + skip[0][1],
97 ((size[1] - e) * 2 - skip[1][0] + skip[1][1]) *
98 ((size[2] - e) * 2 - skip[2][0] + skip[2][1]),
99 b, skip[0]);
100 }
101
102 #endif
103