Skip to content

Commit 9d6e65a

Browse files
DiamonDinoiaDiamonDinoia
committed
drafted-test
1 parent 5dde122 commit 9d6e65a

File tree

2 files changed

+176
-1
lines changed

2 files changed

+176
-1
lines changed

‎test/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
# Each of these source test files is instantiated in single and double precision
2-
set(TESTS basicpassfail dumbinputs finufft1d_test finufft1dmany_test finufft2d_test finufft2dmany_test finufft3d_test finufft3dmany_test)
2+
set(TESTS basicpassfail dumbinputs finufft1d_test finufft1dmany_test finufft2d_test finufft2dmany_test finufft3d_test finufft3dmany_test finufft3dkernel_test)
33

44
foreach(TEST ${TESTS})
55
add_executable(${TEST} ${TEST}.cpp)

‎test/finufft3dkernel_test.cpp‎

Lines changed: 175 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,175 @@
1+
#include <finufft/test_defs.h>
2+
// this enforces recompilation, responding to SINGLE...
3+
#include "directft/dirft3d.cpp"
4+
using namespace std;
5+
using namespace finufft::utils;
6+
7+
const char *help[] = {"Tester for FINUFFT in 3d, all 3 types, either precision.",
8+
"",
9+
"Usage: finufft3d_test Nmodes1 Nmodes2 Nmodes3 Nsrc [tol [debug "
10+
"[spread_sort [upsampfac [errfail]]]]]",
11+
"\teg:\tfinufft3d_test 100 200 50 1e6 1e-12 0 2 0.0 1e-11",
12+
"\tnotes:\tif errfail present, exit code 1 if any error > errfail",
13+
NULL};
14+
// Barnett 2/2/17 onwards.
15+
16+
int main(int argc, char *argv[]) {
17+
BIGINT M, N1, N2, N3; // M = # srcs, N1,N2,N3 = # modes
18+
double w, tol = 1e-6; // default
19+
double err, errmax = 0;
20+
finufft_opts opts0, opts1;
21+
FINUFFT_DEFAULT_OPTS(&opts0);
22+
FINUFFT_DEFAULT_OPTS(&opts1);
23+
opts0.spread_kerevalmeth = 0;
24+
opts1.spread_kerevalmeth = 1;
25+
// opts.fftw = FFTW_MEASURE; // change from usual FFTW_ESTIMATE
26+
// opts.spread_max_sp_size = 3e4; // override test
27+
// opts.spread_nthr_atomic = 15; // "
28+
int isign = +1; // choose which exponential sign to test
29+
if (argc < 5 || argc > 10) {
30+
for (int i = 0; help[i]; ++i) fprintf(stderr, "%s\n", help[i]);
31+
return 2;
32+
}
33+
sscanf(argv[1], "%lf", &w);
34+
N1 = (BIGINT)w;
35+
sscanf(argv[2], "%lf", &w);
36+
N2 = (BIGINT)w;
37+
sscanf(argv[3], "%lf", &w);
38+
N3 = (BIGINT)w;
39+
sscanf(argv[4], "%lf", &w);
40+
M = (BIGINT)w;
41+
if (argc > 5) sscanf(argv[5], "%lf", &tol);
42+
if (argc > 6) sscanf(argv[6], "%d", &opts0.debug); // can be 0,1 or 2
43+
opts0.spread_debug = (opts0.debug > 1) ? 1 : 0; // see output from spreader
44+
if (argc > 7) sscanf(argv[7], "%d", &opts0.spread_sort);
45+
if (argc > 8) {
46+
sscanf(argv[8], "%lf", &w);
47+
opts0.upsampfac = (FLT)w;
48+
}
49+
50+
opts1 = opts0;
51+
opts0.spread_kerevalmeth = 0;
52+
opts1.spread_kerevalmeth = 1;
53+
54+
cout << scientific << setprecision(15);
55+
const BIGINT N = N1 * N2 * N3;
56+
57+
std::vector<FLT> x(M); // NU pts x coords
58+
std::vector<FLT> y(M); // NU pts y coords
59+
std::vector<FLT> z(M); // NU pts z coords
60+
std::vector<CPX> c0(M), c1(N); // strengths
61+
std::vector<CPX> F0(N); // mode ampls kereval 0
62+
std::vector<CPX> F1(N); // mode ampls kereval 1
63+
#pragma omp parallel
64+
{
65+
unsigned int se = MY_OMP_GET_THREAD_NUM(); // needed for parallel random #s
66+
#pragma omp for schedule(static, TEST_RANDCHUNK)
67+
for (BIGINT j = 0; j < M; ++j) {
68+
x[j] = M_PI * randm11r(&se);
69+
y[j] = M_PI * randm11r(&se);
70+
z[j] = M_PI * randm11r(&se);
71+
c0[j] = crandm11r(&se);
72+
}
73+
}
74+
c1 = c0; // copy strengths
75+
printf("test 3d type 1:\n"); // -------------- type 1
76+
printf("kerevalmeth 0:\n");
77+
CNTime timer{};
78+
timer.start();
79+
int ier = FINUFFT3D1(M, x.data(), y.data(), z.data(), c0.data(), isign, tol, N1, N2, N3,
80+
F0.data(), &opts0);
81+
double ti = timer.elapsedsec();
82+
if (ier > 1) {
83+
printf("error (ier=%d)!\n", ier);
84+
return ier;
85+
} else
86+
printf(" %lld NU pts to (%lld,%lld,%lld) modes in %.3g s \t%.3g NU pts/s\n",
87+
(long long)M, (long long)N1, (long long)N2, (long long)N3, ti, M / ti);
88+
printf("kerevalmeth 1:\n");
89+
timer.restart();
90+
ier = FINUFFT3D1(M, x.data(), y.data(), z.data(), c0.data(), isign, tol, N1, N2, N3,
91+
F1.data(), &opts1);
92+
ti = timer.elapsedsec();
93+
if (ier > 1) {
94+
printf("error (ier=%d)!\n", ier);
95+
return ier;
96+
} else
97+
printf(" %lld NU pts to (%lld,%lld,%lld) modes in %.3g s \t%.3g NU pts/s\n",
98+
(long long)M, (long long)N1, (long long)N2, (long long)N3, ti, M / ti);
99+
100+
err = relerrtwonorm(N, F0.data(), F1.data());
101+
errmax = max(err, errmax);
102+
printf("\ttype 1 rel l2-err in F is %.3g\n", err);
103+
// copy F0 to F1 so that we can test type 2
104+
F1 = F0;
105+
printf("kerevalmeth 0:\n");
106+
timer.restart();
107+
ier = FINUFFT3D2(M, x.data(), y.data(), z.data(), c0.data(), isign, tol, N1, N2, N3,
108+
F0.data(), &opts0);
109+
ti = timer.elapsedsec();
110+
if (ier > 1) {
111+
printf("error (ier=%d)!\n", ier);
112+
return ier;
113+
} else
114+
printf(" (%lld,%lld,%lld) modes to %lld NU pts in %.3g s \t%.3g NU pts/s\n",
115+
(long long)N1, (long long)N2, (long long)N3, (long long)M, ti, M / ti);
116+
printf("kerevalmeth 1:\n");
117+
timer.restart();
118+
ier = FINUFFT3D2(M, x.data(), y.data(), z.data(), c1.data(), isign, tol, N1, N2, N3,
119+
F0.data(), &opts1);
120+
ti = timer.elapsedsec();
121+
if (ier > 1) {
122+
printf("error (ier=%d)!\n", ier);
123+
return ier;
124+
} else
125+
printf(" (%lld,%lld,%lld) modes to %lld NU pts in %.3g s \t%.3g NU pts/s\n",
126+
(long long)N1, (long long)N2, (long long)N3, (long long)M, ti, M / ti);
127+
err = relerrtwonorm(M, c0.data(), c1.data());
128+
errmax = std::max(err, errmax);
129+
printf("\ttype 2 rel l2-err in c is %.3g\n", err);
130+
131+
printf("test 3d type 3:\n"); // -------------- type 3
132+
#pragma omp parallel
133+
{
134+
unsigned int se = MY_OMP_GET_THREAD_NUM();
135+
#pragma omp for schedule(static, TEST_RANDCHUNK)
136+
for (BIGINT j = 0; j < M; ++j) {
137+
x[j] = 2.0 + M_PI * randm11r(&se); // new x_j srcs, offset from origin
138+
y[j] = -3.0 + M_PI * randm11r(&se); // " y_j
139+
z[j] = 1.0 + M_PI * randm11r(&se); // " z_j
140+
}
141+
}
142+
std::vector<FLT> s(N); // targ freqs (1-cmpt)
143+
std::vector<FLT> t(N); // targ freqs (2-cmpt)
144+
std::vector<FLT> u(N); // targ freqs (3-cmpt)
145+
FLT S1 = (FLT)N1 / 2; // choose freq range sim to type 1
146+
FLT S2 = (FLT)N2 / 2;
147+
FLT S3 = (FLT)N3 / 2;
148+
149+
timer.restart();
150+
printf("kerevalmeth 0:\n");
151+
ier = FINUFFT3D3(M, x.data(), y.data(), z.data(), c0.data(), isign, tol, N, s.data(),
152+
t.data(), u.data(), F0.data(), &opts0);
153+
ti = timer.elapsedsec();
154+
if (ier > 1) {
155+
printf("error (ier=%d)!\n", ier);
156+
return ier;
157+
} else
158+
printf("\t%lld NU to %lld NU in %.3g s \t%.3g tot NU pts/s\n", (long long)M,
159+
(long long)N, ti, (M + N) / ti);
160+
timer.restart();
161+
printf("kerevalmeth 1:\n");
162+
ier = FINUFFT3D3(M, x.data(), y.data(), z.data(), c0.data(), isign, tol, N, s.data(),
163+
t.data(), u.data(), F1.data(), &opts1);
164+
ti = timer.elapsedsec();
165+
if (ier > 1) {
166+
printf("error (ier=%d)!\n", ier);
167+
return ier;
168+
} else
169+
printf("\t%lld NU to %lld NU in %.3g s \t%.3g tot NU pts/s\n", (long long)M,
170+
(long long)N, ti, (M + N) / ti);
171+
err = relerrtwonorm(N, F0.data(), F1.data());
172+
errmax = max(err, errmax);
173+
printf("\ttype 3 rel l2-err in F is %.3g\n", err);
174+
return (errmax > tol);
175+
}

0 commit comments

Comments
 (0)