Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.8/libbench2/verify-r2r.c @ 84:08ae793730bd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add null config files
author Chris Cannam
date Mon, 02 Mar 2020 14:03:47 +0000
parents d0c2a83c1364
children
rev   line source
Chris@82 1 /*
Chris@82 2 * Copyright (c) 2003, 2007-14 Matteo Frigo
Chris@82 3 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
Chris@82 4 *
Chris@82 5 * This program is free software; you can redistribute it and/or modify
Chris@82 6 * it under the terms of the GNU General Public License as published by
Chris@82 7 * the Free Software Foundation; either version 2 of the License, or
Chris@82 8 * (at your option) any later version.
Chris@82 9 *
Chris@82 10 * This program is distributed in the hope that it will be useful,
Chris@82 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@82 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
Chris@82 13 * GNU General Public License for more details.
Chris@82 14 *
Chris@82 15 * You should have received a copy of the GNU General Public License
Chris@82 16 * along with this program; if not, write to the Free Software
Chris@82 17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Chris@82 18 *
Chris@82 19 */
Chris@82 20
Chris@82 21 /* Lots of ugly duplication from verify-lib.c, plus lots of ugliness in
Chris@82 22 general for all of the r2r variants...oh well, for now */
Chris@82 23
Chris@82 24 #include "verify.h"
Chris@82 25 #include <math.h>
Chris@82 26 #include <stdlib.h>
Chris@82 27 #include <stdio.h>
Chris@82 28
Chris@82 29 typedef struct {
Chris@82 30 bench_problem *p;
Chris@82 31 bench_tensor *probsz;
Chris@82 32 bench_tensor *totalsz;
Chris@82 33 bench_tensor *pckdsz;
Chris@82 34 bench_tensor *pckdvecsz;
Chris@82 35 } info;
Chris@82 36
Chris@82 37 /*
Chris@82 38 * Utility functions:
Chris@82 39 */
Chris@82 40
Chris@82 41 static double dabs(double x) { return (x < 0.0) ? -x : x; }
Chris@82 42 static double dmin(double x, double y) { return (x < y) ? x : y; }
Chris@82 43
Chris@82 44 static double raerror(R *a, R *b, int n)
Chris@82 45 {
Chris@82 46 if (n > 0) {
Chris@82 47 /* compute the relative Linf error */
Chris@82 48 double e = 0.0, mag = 0.0;
Chris@82 49 int i;
Chris@82 50
Chris@82 51 for (i = 0; i < n; ++i) {
Chris@82 52 e = dmax(e, dabs(a[i] - b[i]));
Chris@82 53 mag = dmax(mag, dmin(dabs(a[i]), dabs(b[i])));
Chris@82 54 }
Chris@82 55 if (dabs(mag) < 1e-14 && dabs(e) < 1e-14)
Chris@82 56 e = 0.0;
Chris@82 57 else
Chris@82 58 e /= mag;
Chris@82 59
Chris@82 60 #ifdef HAVE_ISNAN
Chris@82 61 BENCH_ASSERT(!isnan(e));
Chris@82 62 #endif
Chris@82 63 return e;
Chris@82 64 } else
Chris@82 65 return 0.0;
Chris@82 66 }
Chris@82 67
Chris@82 68 #define by2pi(m, n) ((K2PI * (m)) / (n))
Chris@82 69
Chris@82 70 /*
Chris@82 71 * Improve accuracy by reducing x to range [0..1/8]
Chris@82 72 * before multiplication by 2 * PI.
Chris@82 73 */
Chris@82 74
Chris@82 75 static trigreal bench_sincos(trigreal m, trigreal n, int sinp)
Chris@82 76 {
Chris@82 77 /* waiting for C to get tail recursion... */
Chris@82 78 trigreal half_n = n * 0.5;
Chris@82 79 trigreal quarter_n = half_n * 0.5;
Chris@82 80 trigreal eighth_n = quarter_n * 0.5;
Chris@82 81 trigreal sgn = 1.0;
Chris@82 82
Chris@82 83 if (sinp) goto sin;
Chris@82 84 cos:
Chris@82 85 if (m < 0) { m = -m; /* goto cos; */ }
Chris@82 86 if (m > half_n) { m = n - m; goto cos; }
Chris@82 87 if (m > eighth_n) { m = quarter_n - m; goto sin; }
Chris@82 88 return sgn * COS(by2pi(m, n));
Chris@82 89
Chris@82 90 msin:
Chris@82 91 sgn = -sgn;
Chris@82 92 sin:
Chris@82 93 if (m < 0) { m = -m; goto msin; }
Chris@82 94 if (m > half_n) { m = n - m; goto msin; }
Chris@82 95 if (m > eighth_n) { m = quarter_n - m; goto cos; }
Chris@82 96 return sgn * SIN(by2pi(m, n));
Chris@82 97 }
Chris@82 98
Chris@82 99 static trigreal cos2pi(int m, int n)
Chris@82 100 {
Chris@82 101 return bench_sincos((trigreal)m, (trigreal)n, 0);
Chris@82 102 }
Chris@82 103
Chris@82 104 static trigreal sin2pi(int m, int n)
Chris@82 105 {
Chris@82 106 return bench_sincos((trigreal)m, (trigreal)n, 1);
Chris@82 107 }
Chris@82 108
Chris@82 109 static trigreal cos00(int i, int j, int n)
Chris@82 110 {
Chris@82 111 return cos2pi(i * j, n);
Chris@82 112 }
Chris@82 113
Chris@82 114 static trigreal cos01(int i, int j, int n)
Chris@82 115 {
Chris@82 116 return cos00(i, 2*j + 1, 2*n);
Chris@82 117 }
Chris@82 118
Chris@82 119 static trigreal cos10(int i, int j, int n)
Chris@82 120 {
Chris@82 121 return cos00(2*i + 1, j, 2*n);
Chris@82 122 }
Chris@82 123
Chris@82 124 static trigreal cos11(int i, int j, int n)
Chris@82 125 {
Chris@82 126 return cos00(2*i + 1, 2*j + 1, 4*n);
Chris@82 127 }
Chris@82 128
Chris@82 129 static trigreal sin00(int i, int j, int n)
Chris@82 130 {
Chris@82 131 return sin2pi(i * j, n);
Chris@82 132 }
Chris@82 133
Chris@82 134 static trigreal sin01(int i, int j, int n)
Chris@82 135 {
Chris@82 136 return sin00(i, 2*j + 1, 2*n);
Chris@82 137 }
Chris@82 138
Chris@82 139 static trigreal sin10(int i, int j, int n)
Chris@82 140 {
Chris@82 141 return sin00(2*i + 1, j, 2*n);
Chris@82 142 }
Chris@82 143
Chris@82 144 static trigreal sin11(int i, int j, int n)
Chris@82 145 {
Chris@82 146 return sin00(2*i + 1, 2*j + 1, 4*n);
Chris@82 147 }
Chris@82 148
Chris@82 149 static trigreal realhalf(int i, int j, int n)
Chris@82 150 {
Chris@82 151 UNUSED(i);
Chris@82 152 if (j <= n - j)
Chris@82 153 return 1.0;
Chris@82 154 else
Chris@82 155 return 0.0;
Chris@82 156 }
Chris@82 157
Chris@82 158 static trigreal coshalf(int i, int j, int n)
Chris@82 159 {
Chris@82 160 if (j <= n - j)
Chris@82 161 return cos00(i, j, n);
Chris@82 162 else
Chris@82 163 return cos00(i, n - j, n);
Chris@82 164 }
Chris@82 165
Chris@82 166 static trigreal unity(int i, int j, int n)
Chris@82 167 {
Chris@82 168 UNUSED(i);
Chris@82 169 UNUSED(j);
Chris@82 170 UNUSED(n);
Chris@82 171 return 1.0;
Chris@82 172 }
Chris@82 173
Chris@82 174 typedef trigreal (*trigfun)(int, int, int);
Chris@82 175
Chris@82 176 static void rarand(R *a, int n)
Chris@82 177 {
Chris@82 178 int i;
Chris@82 179
Chris@82 180 /* generate random inputs */
Chris@82 181 for (i = 0; i < n; ++i) {
Chris@82 182 a[i] = mydrand();
Chris@82 183 }
Chris@82 184 }
Chris@82 185
Chris@82 186 /* C = A + B */
Chris@82 187 static void raadd(R *c, R *a, R *b, int n)
Chris@82 188 {
Chris@82 189 int i;
Chris@82 190
Chris@82 191 for (i = 0; i < n; ++i) {
Chris@82 192 c[i] = a[i] + b[i];
Chris@82 193 }
Chris@82 194 }
Chris@82 195
Chris@82 196 /* C = A - B */
Chris@82 197 static void rasub(R *c, R *a, R *b, int n)
Chris@82 198 {
Chris@82 199 int i;
Chris@82 200
Chris@82 201 for (i = 0; i < n; ++i) {
Chris@82 202 c[i] = a[i] - b[i];
Chris@82 203 }
Chris@82 204 }
Chris@82 205
Chris@82 206 /* B = rotate left A + rotate right A */
Chris@82 207 static void rarolr(R *b, R *a, int n, int nb, int na,
Chris@82 208 r2r_kind_t k)
Chris@82 209 {
Chris@82 210 int isL0 = 0, isL1 = 0, isR0 = 0, isR1 = 0;
Chris@82 211 int i, ib, ia;
Chris@82 212
Chris@82 213 for (ib = 0; ib < nb; ++ib) {
Chris@82 214 for (i = 0; i < n - 1; ++i)
Chris@82 215 for (ia = 0; ia < na; ++ia)
Chris@82 216 b[(ib * n + i) * na + ia] =
Chris@82 217 a[(ib * n + i + 1) * na + ia];
Chris@82 218
Chris@82 219 /* ugly switch to do boundary conditions for various r2r types */
Chris@82 220 switch (k) {
Chris@82 221 /* periodic boundaries */
Chris@82 222 case R2R_DHT:
Chris@82 223 case R2R_R2HC:
Chris@82 224 for (ia = 0; ia < na; ++ia) {
Chris@82 225 b[(ib * n + n - 1) * na + ia] =
Chris@82 226 a[(ib * n + 0) * na + ia];
Chris@82 227 b[(ib * n + 0) * na + ia] +=
Chris@82 228 a[(ib * n + n - 1) * na + ia];
Chris@82 229 }
Chris@82 230 break;
Chris@82 231
Chris@82 232 case R2R_HC2R: /* ugh (hermitian halfcomplex boundaries) */
Chris@82 233 if (n > 2) {
Chris@82 234 if (n % 2 == 0)
Chris@82 235 for (ia = 0; ia < na; ++ia) {
Chris@82 236 b[(ib * n + n - 1) * na + ia] = 0.0;
Chris@82 237 b[(ib * n + 0) * na + ia] +=
Chris@82 238 a[(ib * n + 1) * na + ia];
Chris@82 239 b[(ib * n + n/2) * na + ia] +=
Chris@82 240 + a[(ib * n + n/2 - 1) * na + ia]
Chris@82 241 - a[(ib * n + n/2 + 1) * na + ia];
Chris@82 242 b[(ib * n + n/2 + 1) * na + ia] +=
Chris@82 243 - a[(ib * n + n/2) * na + ia];
Chris@82 244 }
Chris@82 245 else
Chris@82 246 for (ia = 0; ia < na; ++ia) {
Chris@82 247 b[(ib * n + n - 1) * na + ia] = 0.0;
Chris@82 248 b[(ib * n + 0) * na + ia] +=
Chris@82 249 a[(ib * n + 1) * na + ia];
Chris@82 250 b[(ib * n + n/2) * na + ia] +=
Chris@82 251 + a[(ib * n + n/2) * na + ia]
Chris@82 252 - a[(ib * n + n/2 + 1) * na + ia];
Chris@82 253 b[(ib * n + n/2 + 1) * na + ia] +=
Chris@82 254 - a[(ib * n + n/2 + 1) * na + ia]
Chris@82 255 - a[(ib * n + n/2) * na + ia];
Chris@82 256 }
Chris@82 257 } else /* n <= 2 */ {
Chris@82 258 for (ia = 0; ia < na; ++ia) {
Chris@82 259 b[(ib * n + n - 1) * na + ia] =
Chris@82 260 a[(ib * n + 0) * na + ia];
Chris@82 261 b[(ib * n + 0) * na + ia] +=
Chris@82 262 a[(ib * n + n - 1) * na + ia];
Chris@82 263 }
Chris@82 264 }
Chris@82 265 break;
Chris@82 266
Chris@82 267 /* various even/odd boundary conditions */
Chris@82 268 case R2R_REDFT00:
Chris@82 269 isL1 = isR1 = 1;
Chris@82 270 goto mirrors;
Chris@82 271 case R2R_REDFT01:
Chris@82 272 isL1 = 1;
Chris@82 273 goto mirrors;
Chris@82 274 case R2R_REDFT10:
Chris@82 275 isL0 = isR0 = 1;
Chris@82 276 goto mirrors;
Chris@82 277 case R2R_REDFT11:
Chris@82 278 isL0 = 1;
Chris@82 279 isR0 = -1;
Chris@82 280 goto mirrors;
Chris@82 281 case R2R_RODFT00:
Chris@82 282 goto mirrors;
Chris@82 283 case R2R_RODFT01:
Chris@82 284 isR1 = 1;
Chris@82 285 goto mirrors;
Chris@82 286 case R2R_RODFT10:
Chris@82 287 isL0 = isR0 = -1;
Chris@82 288 goto mirrors;
Chris@82 289 case R2R_RODFT11:
Chris@82 290 isL0 = -1;
Chris@82 291 isR0 = 1;
Chris@82 292 goto mirrors;
Chris@82 293
Chris@82 294 mirrors:
Chris@82 295
Chris@82 296 for (ia = 0; ia < na; ++ia)
Chris@82 297 b[(ib * n + n - 1) * na + ia] =
Chris@82 298 isR0 * a[(ib * n + n - 1) * na + ia]
Chris@82 299 + (n > 1 ? isR1 * a[(ib * n + n - 2) * na + ia]
Chris@82 300 : 0);
Chris@82 301
Chris@82 302 for (ia = 0; ia < na; ++ia)
Chris@82 303 b[(ib * n) * na + ia] +=
Chris@82 304 isL0 * a[(ib * n) * na + ia]
Chris@82 305 + (n > 1 ? isL1 * a[(ib * n + 1) * na + ia] : 0);
Chris@82 306
Chris@82 307 }
Chris@82 308
Chris@82 309 for (i = 1; i < n; ++i)
Chris@82 310 for (ia = 0; ia < na; ++ia)
Chris@82 311 b[(ib * n + i) * na + ia] +=
Chris@82 312 a[(ib * n + i - 1) * na + ia];
Chris@82 313 }
Chris@82 314 }
Chris@82 315
Chris@82 316 static void raphase_shift(R *b, R *a, int n, int nb, int na,
Chris@82 317 int n0, int k0, trigfun t)
Chris@82 318 {
Chris@82 319 int j, jb, ja;
Chris@82 320
Chris@82 321 for (jb = 0; jb < nb; ++jb)
Chris@82 322 for (j = 0; j < n; ++j) {
Chris@82 323 trigreal c = 2.0 * t(1, j + k0, n0);
Chris@82 324
Chris@82 325 for (ja = 0; ja < na; ++ja) {
Chris@82 326 int k = (jb * n + j) * na + ja;
Chris@82 327 b[k] = a[k] * c;
Chris@82 328 }
Chris@82 329 }
Chris@82 330 }
Chris@82 331
Chris@82 332 /* A = alpha * A (real, in place) */
Chris@82 333 static void rascale(R *a, R alpha, int n)
Chris@82 334 {
Chris@82 335 int i;
Chris@82 336
Chris@82 337 for (i = 0; i < n; ++i) {
Chris@82 338 a[i] *= alpha;
Chris@82 339 }
Chris@82 340 }
Chris@82 341
Chris@82 342 /*
Chris@82 343 * compute rdft:
Chris@82 344 */
Chris@82 345
Chris@82 346 /* copy real A into real B, using output stride of A and input stride of B */
Chris@82 347 typedef struct {
Chris@82 348 dotens2_closure k;
Chris@82 349 R *ra;
Chris@82 350 R *rb;
Chris@82 351 } cpyr_closure;
Chris@82 352
Chris@82 353 static void cpyr0(dotens2_closure *k_,
Chris@82 354 int indxa, int ondxa, int indxb, int ondxb)
Chris@82 355 {
Chris@82 356 cpyr_closure *k = (cpyr_closure *)k_;
Chris@82 357 k->rb[indxb] = k->ra[ondxa];
Chris@82 358 UNUSED(indxa); UNUSED(ondxb);
Chris@82 359 }
Chris@82 360
Chris@82 361 static void cpyr(R *ra, bench_tensor *sza, R *rb, bench_tensor *szb)
Chris@82 362 {
Chris@82 363 cpyr_closure k;
Chris@82 364 k.k.apply = cpyr0;
Chris@82 365 k.ra = ra; k.rb = rb;
Chris@82 366 bench_dotens2(sza, szb, &k.k);
Chris@82 367 }
Chris@82 368
Chris@82 369 static void dofft(info *nfo, R *in, R *out)
Chris@82 370 {
Chris@82 371 cpyr(in, nfo->pckdsz, (R *) nfo->p->in, nfo->totalsz);
Chris@82 372 after_problem_rcopy_from(nfo->p, (bench_real *)nfo->p->in);
Chris@82 373 doit(1, nfo->p);
Chris@82 374 after_problem_rcopy_to(nfo->p, (bench_real *)nfo->p->out);
Chris@82 375 cpyr((R *) nfo->p->out, nfo->totalsz, out, nfo->pckdsz);
Chris@82 376 }
Chris@82 377
Chris@82 378 static double racmp(R *a, R *b, int n, const char *test, double tol)
Chris@82 379 {
Chris@82 380 double d = raerror(a, b, n);
Chris@82 381 if (d > tol) {
Chris@82 382 ovtpvt_err("Found relative error %e (%s)\n", d, test);
Chris@82 383 {
Chris@82 384 int i, N;
Chris@82 385 N = n > 300 && verbose <= 2 ? 300 : n;
Chris@82 386 for (i = 0; i < N; ++i)
Chris@82 387 ovtpvt_err("%8d %16.12f %16.12f\n", i,
Chris@82 388 (double) a[i],
Chris@82 389 (double) b[i]);
Chris@82 390 }
Chris@82 391 bench_exit(EXIT_FAILURE);
Chris@82 392 }
Chris@82 393 return d;
Chris@82 394 }
Chris@82 395
Chris@82 396 /***********************************************************************/
Chris@82 397
Chris@82 398 typedef struct {
Chris@82 399 int n; /* physical size */
Chris@82 400 int n0; /* "logical" transform size */
Chris@82 401 int i0, k0; /* shifts of input/output */
Chris@82 402 trigfun ti, ts; /* impulse/shift trig functions */
Chris@82 403 } dim_stuff;
Chris@82 404
Chris@82 405 static void impulse_response(int rnk, dim_stuff *d, R impulse_amp,
Chris@82 406 R *A, int N)
Chris@82 407 {
Chris@82 408 if (rnk == 0)
Chris@82 409 A[0] = impulse_amp;
Chris@82 410 else {
Chris@82 411 int i;
Chris@82 412 N /= d->n;
Chris@82 413 for (i = 0; i < d->n; ++i) {
Chris@82 414 impulse_response(rnk - 1, d + 1,
Chris@82 415 impulse_amp * d->ti(d->i0, d->k0 + i, d->n0),
Chris@82 416 A + i * N, N);
Chris@82 417 }
Chris@82 418 }
Chris@82 419 }
Chris@82 420
Chris@82 421 /***************************************************************************/
Chris@82 422
Chris@82 423 /*
Chris@82 424 * Implementation of the FFT tester described in
Chris@82 425 *
Chris@82 426 * Funda Ergün. Testing multivariate linear functions: Overcoming the
Chris@82 427 * generator bottleneck. In Proceedings of the Twenty-Seventh Annual
Chris@82 428 * ACM Symposium on the Theory of Computing, pages 407-416, Las Vegas,
Chris@82 429 * Nevada, 29 May--1 June 1995.
Chris@82 430 *
Chris@82 431 * Also: F. Ergun, S. R. Kumar, and D. Sivakumar, "Self-testing without
Chris@82 432 * the generator bottleneck," SIAM J. on Computing 29 (5), 1630-51 (2000).
Chris@82 433 */
Chris@82 434
Chris@82 435 static double rlinear(int n, info *nfo, R *inA, R *inB, R *inC, R *outA,
Chris@82 436 R *outB, R *outC, R *tmp, int rounds, double tol)
Chris@82 437 {
Chris@82 438 double e = 0.0;
Chris@82 439 int j;
Chris@82 440
Chris@82 441 for (j = 0; j < rounds; ++j) {
Chris@82 442 R alpha, beta;
Chris@82 443 alpha = mydrand();
Chris@82 444 beta = mydrand();
Chris@82 445 rarand(inA, n);
Chris@82 446 rarand(inB, n);
Chris@82 447 dofft(nfo, inA, outA);
Chris@82 448 dofft(nfo, inB, outB);
Chris@82 449
Chris@82 450 rascale(outA, alpha, n);
Chris@82 451 rascale(outB, beta, n);
Chris@82 452 raadd(tmp, outA, outB, n);
Chris@82 453 rascale(inA, alpha, n);
Chris@82 454 rascale(inB, beta, n);
Chris@82 455 raadd(inC, inA, inB, n);
Chris@82 456 dofft(nfo, inC, outC);
Chris@82 457
Chris@82 458 e = dmax(e, racmp(outC, tmp, n, "linear", tol));
Chris@82 459 }
Chris@82 460 return e;
Chris@82 461 }
Chris@82 462
Chris@82 463 static double rimpulse(dim_stuff *d, R impulse_amp,
Chris@82 464 int n, int vecn, info *nfo,
Chris@82 465 R *inA, R *inB, R *inC,
Chris@82 466 R *outA, R *outB, R *outC,
Chris@82 467 R *tmp, int rounds, double tol)
Chris@82 468 {
Chris@82 469 double e = 0.0;
Chris@82 470 int N = n * vecn;
Chris@82 471 int i;
Chris@82 472 int j;
Chris@82 473
Chris@82 474 /* test 2: check that the unit impulse is transformed properly */
Chris@82 475
Chris@82 476 for (i = 0; i < N; ++i) {
Chris@82 477 /* pls */
Chris@82 478 inA[i] = 0.0;
Chris@82 479 }
Chris@82 480 for (i = 0; i < vecn; ++i) {
Chris@82 481 inA[i * n] = (i+1) / (double)(vecn+1);
Chris@82 482
Chris@82 483 /* transform of the pls */
Chris@82 484 impulse_response(nfo->probsz->rnk, d, impulse_amp * inA[i * n],
Chris@82 485 outA + i * n, n);
Chris@82 486 }
Chris@82 487
Chris@82 488 dofft(nfo, inA, tmp);
Chris@82 489 e = dmax(e, racmp(tmp, outA, N, "impulse 1", tol));
Chris@82 490
Chris@82 491 for (j = 0; j < rounds; ++j) {
Chris@82 492 rarand(inB, N);
Chris@82 493 rasub(inC, inA, inB, N);
Chris@82 494 dofft(nfo, inB, outB);
Chris@82 495 dofft(nfo, inC, outC);
Chris@82 496 raadd(tmp, outB, outC, N);
Chris@82 497 e = dmax(e, racmp(tmp, outA, N, "impulse", tol));
Chris@82 498 }
Chris@82 499 return e;
Chris@82 500 }
Chris@82 501
Chris@82 502 static double t_shift(int n, int vecn, info *nfo,
Chris@82 503 R *inA, R *inB, R *outA, R *outB, R *tmp,
Chris@82 504 int rounds, double tol,
Chris@82 505 dim_stuff *d)
Chris@82 506 {
Chris@82 507 double e = 0.0;
Chris@82 508 int nb, na, dim, N = n * vecn;
Chris@82 509 int i, j;
Chris@82 510 bench_tensor *sz = nfo->probsz;
Chris@82 511
Chris@82 512 /* test 3: check the time-shift property */
Chris@82 513 /* the paper performs more tests, but this code should be fine too */
Chris@82 514
Chris@82 515 nb = 1;
Chris@82 516 na = n;
Chris@82 517
Chris@82 518 /* check shifts across all SZ dimensions */
Chris@82 519 for (dim = 0; dim < sz->rnk; ++dim) {
Chris@82 520 int ncur = sz->dims[dim].n;
Chris@82 521
Chris@82 522 na /= ncur;
Chris@82 523
Chris@82 524 for (j = 0; j < rounds; ++j) {
Chris@82 525 rarand(inA, N);
Chris@82 526
Chris@82 527 for (i = 0; i < vecn; ++i) {
Chris@82 528 rarolr(inB + i * n, inA + i*n, ncur, nb,na,
Chris@82 529 nfo->p->k[dim]);
Chris@82 530 }
Chris@82 531 dofft(nfo, inA, outA);
Chris@82 532 dofft(nfo, inB, outB);
Chris@82 533 for (i = 0; i < vecn; ++i)
Chris@82 534 raphase_shift(tmp + i * n, outA + i * n, ncur,
Chris@82 535 nb, na, d[dim].n0, d[dim].k0, d[dim].ts);
Chris@82 536 e = dmax(e, racmp(tmp, outB, N, "time shift", tol));
Chris@82 537 }
Chris@82 538
Chris@82 539 nb *= ncur;
Chris@82 540 }
Chris@82 541 return e;
Chris@82 542 }
Chris@82 543
Chris@82 544 /***********************************************************************/
Chris@82 545
Chris@82 546 void verify_r2r(bench_problem *p, int rounds, double tol, errors *e)
Chris@82 547 {
Chris@82 548 R *inA, *inB, *inC, *outA, *outB, *outC, *tmp;
Chris@82 549 info nfo;
Chris@82 550 int n, vecn, N;
Chris@82 551 double impulse_amp = 1.0;
Chris@82 552 dim_stuff *d;
Chris@82 553 int i;
Chris@82 554
Chris@82 555 if (rounds == 0)
Chris@82 556 rounds = 20; /* default value */
Chris@82 557
Chris@82 558 n = tensor_sz(p->sz);
Chris@82 559 vecn = tensor_sz(p->vecsz);
Chris@82 560 N = n * vecn;
Chris@82 561
Chris@82 562 d = (dim_stuff *) bench_malloc(sizeof(dim_stuff) * p->sz->rnk);
Chris@82 563 for (i = 0; i < p->sz->rnk; ++i) {
Chris@82 564 int n0, i0, k0;
Chris@82 565 trigfun ti, ts;
Chris@82 566
Chris@82 567 d[i].n = n0 = p->sz->dims[i].n;
Chris@82 568 if (p->k[i] > R2R_DHT)
Chris@82 569 n0 = 2 * (n0 + (p->k[i] == R2R_REDFT00 ? -1 :
Chris@82 570 (p->k[i] == R2R_RODFT00 ? 1 : 0)));
Chris@82 571
Chris@82 572 switch (p->k[i]) {
Chris@82 573 case R2R_R2HC:
Chris@82 574 i0 = k0 = 0;
Chris@82 575 ti = realhalf;
Chris@82 576 ts = coshalf;
Chris@82 577 break;
Chris@82 578 case R2R_DHT:
Chris@82 579 i0 = k0 = 0;
Chris@82 580 ti = unity;
Chris@82 581 ts = cos00;
Chris@82 582 break;
Chris@82 583 case R2R_HC2R:
Chris@82 584 i0 = k0 = 0;
Chris@82 585 ti = unity;
Chris@82 586 ts = cos00;
Chris@82 587 break;
Chris@82 588 case R2R_REDFT00:
Chris@82 589 i0 = k0 = 0;
Chris@82 590 ti = ts = cos00;
Chris@82 591 break;
Chris@82 592 case R2R_REDFT01:
Chris@82 593 i0 = k0 = 0;
Chris@82 594 ti = ts = cos01;
Chris@82 595 break;
Chris@82 596 case R2R_REDFT10:
Chris@82 597 i0 = k0 = 0;
Chris@82 598 ti = cos10; impulse_amp *= 2.0;
Chris@82 599 ts = cos00;
Chris@82 600 break;
Chris@82 601 case R2R_REDFT11:
Chris@82 602 i0 = k0 = 0;
Chris@82 603 ti = cos11; impulse_amp *= 2.0;
Chris@82 604 ts = cos01;
Chris@82 605 break;
Chris@82 606 case R2R_RODFT00:
Chris@82 607 i0 = k0 = 1;
Chris@82 608 ti = sin00; impulse_amp *= 2.0;
Chris@82 609 ts = cos00;
Chris@82 610 break;
Chris@82 611 case R2R_RODFT01:
Chris@82 612 i0 = 1; k0 = 0;
Chris@82 613 ti = sin01; impulse_amp *= n == 1 ? 1.0 : 2.0;
Chris@82 614 ts = cos01;
Chris@82 615 break;
Chris@82 616 case R2R_RODFT10:
Chris@82 617 i0 = 0; k0 = 1;
Chris@82 618 ti = sin10; impulse_amp *= 2.0;
Chris@82 619 ts = cos00;
Chris@82 620 break;
Chris@82 621 case R2R_RODFT11:
Chris@82 622 i0 = k0 = 0;
Chris@82 623 ti = sin11; impulse_amp *= 2.0;
Chris@82 624 ts = cos01;
Chris@82 625 break;
Chris@82 626 default:
Chris@82 627 BENCH_ASSERT(0);
Chris@82 628 return;
Chris@82 629 }
Chris@82 630
Chris@82 631 d[i].n0 = n0;
Chris@82 632 d[i].i0 = i0;
Chris@82 633 d[i].k0 = k0;
Chris@82 634 d[i].ti = ti;
Chris@82 635 d[i].ts = ts;
Chris@82 636 }
Chris@82 637
Chris@82 638
Chris@82 639 inA = (R *) bench_malloc(N * sizeof(R));
Chris@82 640 inB = (R *) bench_malloc(N * sizeof(R));
Chris@82 641 inC = (R *) bench_malloc(N * sizeof(R));
Chris@82 642 outA = (R *) bench_malloc(N * sizeof(R));
Chris@82 643 outB = (R *) bench_malloc(N * sizeof(R));
Chris@82 644 outC = (R *) bench_malloc(N * sizeof(R));
Chris@82 645 tmp = (R *) bench_malloc(N * sizeof(R));
Chris@82 646
Chris@82 647 nfo.p = p;
Chris@82 648 nfo.probsz = p->sz;
Chris@82 649 nfo.totalsz = tensor_append(p->vecsz, nfo.probsz);
Chris@82 650 nfo.pckdsz = verify_pack(nfo.totalsz, 1);
Chris@82 651 nfo.pckdvecsz = verify_pack(p->vecsz, tensor_sz(nfo.probsz));
Chris@82 652
Chris@82 653 e->i = rimpulse(d, impulse_amp, n, vecn, &nfo,
Chris@82 654 inA, inB, inC, outA, outB, outC, tmp, rounds, tol);
Chris@82 655 e->l = rlinear(N, &nfo, inA, inB, inC, outA, outB, outC, tmp, rounds,tol);
Chris@82 656 e->s = t_shift(n, vecn, &nfo, inA, inB, outA, outB, tmp,
Chris@82 657 rounds, tol, d);
Chris@82 658
Chris@82 659 /* grr, verify-lib.c:preserves_input() only works for complex */
Chris@82 660 if (!p->in_place && !p->destroy_input) {
Chris@82 661 bench_tensor *totalsz_swap, *pckdsz_swap;
Chris@82 662 totalsz_swap = tensor_copy_swapio(nfo.totalsz);
Chris@82 663 pckdsz_swap = tensor_copy_swapio(nfo.pckdsz);
Chris@82 664
Chris@82 665 for (i = 0; i < rounds; ++i) {
Chris@82 666 rarand(inA, N);
Chris@82 667 dofft(&nfo, inA, outB);
Chris@82 668 cpyr((R *) nfo.p->in, totalsz_swap, inB, pckdsz_swap);
Chris@82 669 racmp(inB, inA, N, "preserves_input", 0.0);
Chris@82 670 }
Chris@82 671
Chris@82 672 tensor_destroy(totalsz_swap);
Chris@82 673 tensor_destroy(pckdsz_swap);
Chris@82 674 }
Chris@82 675
Chris@82 676 tensor_destroy(nfo.totalsz);
Chris@82 677 tensor_destroy(nfo.pckdsz);
Chris@82 678 tensor_destroy(nfo.pckdvecsz);
Chris@82 679 bench_free(tmp);
Chris@82 680 bench_free(outC);
Chris@82 681 bench_free(outB);
Chris@82 682 bench_free(outA);
Chris@82 683 bench_free(inC);
Chris@82 684 bench_free(inB);
Chris@82 685 bench_free(inA);
Chris@82 686 bench_free(d);
Chris@82 687 }
Chris@82 688
Chris@82 689
Chris@82 690 typedef struct {
Chris@82 691 dofft_closure k;
Chris@82 692 bench_problem *p;
Chris@82 693 int n0;
Chris@82 694 } dofft_r2r_closure;
Chris@82 695
Chris@82 696 static void cpyr1(int n, R *in, int is, R *out, int os, R scale)
Chris@82 697 {
Chris@82 698 int i;
Chris@82 699 for (i = 0; i < n; ++i)
Chris@82 700 out[i * os] = in[i * is] * scale;
Chris@82 701 }
Chris@82 702
Chris@82 703 static void mke00(C *a, int n, int c)
Chris@82 704 {
Chris@82 705 int i;
Chris@82 706 for (i = 1; i + i < n; ++i)
Chris@82 707 a[n - i][c] = a[i][c];
Chris@82 708 }
Chris@82 709
Chris@82 710 static void mkre00(C *a, int n)
Chris@82 711 {
Chris@82 712 mkreal(a, n);
Chris@82 713 mke00(a, n, 0);
Chris@82 714 }
Chris@82 715
Chris@82 716 static void mkimag(C *a, int n)
Chris@82 717 {
Chris@82 718 int i;
Chris@82 719 for (i = 0; i < n; ++i)
Chris@82 720 c_re(a[i]) = 0.0;
Chris@82 721 }
Chris@82 722
Chris@82 723 static void mko00(C *a, int n, int c)
Chris@82 724 {
Chris@82 725 int i;
Chris@82 726 a[0][c] = 0.0;
Chris@82 727 for (i = 1; i + i < n; ++i)
Chris@82 728 a[n - i][c] = -a[i][c];
Chris@82 729 if (i + i == n)
Chris@82 730 a[i][c] = 0.0;
Chris@82 731 }
Chris@82 732
Chris@82 733 static void mkro00(C *a, int n)
Chris@82 734 {
Chris@82 735 mkreal(a, n);
Chris@82 736 mko00(a, n, 0);
Chris@82 737 }
Chris@82 738
Chris@82 739 static void mkio00(C *a, int n)
Chris@82 740 {
Chris@82 741 mkimag(a, n);
Chris@82 742 mko00(a, n, 1);
Chris@82 743 }
Chris@82 744
Chris@82 745 static void mkre01(C *a, int n) /* n should be be multiple of 4 */
Chris@82 746 {
Chris@82 747 R a0;
Chris@82 748 a0 = c_re(a[0]);
Chris@82 749 mko00(a, n/2, 0);
Chris@82 750 c_re(a[n/2]) = -(c_re(a[0]) = a0);
Chris@82 751 mkre00(a, n);
Chris@82 752 }
Chris@82 753
Chris@82 754 static void mkro01(C *a, int n) /* n should be be multiple of 4 */
Chris@82 755 {
Chris@82 756 c_re(a[0]) = c_im(a[0]) = 0.0;
Chris@82 757 mkre00(a, n/2);
Chris@82 758 mkro00(a, n);
Chris@82 759 }
Chris@82 760
Chris@82 761 static void mkoddonly(C *a, int n)
Chris@82 762 {
Chris@82 763 int i;
Chris@82 764 for (i = 0; i < n; i += 2)
Chris@82 765 c_re(a[i]) = c_im(a[i]) = 0.0;
Chris@82 766 }
Chris@82 767
Chris@82 768 static void mkre10(C *a, int n)
Chris@82 769 {
Chris@82 770 mkoddonly(a, n);
Chris@82 771 mkre00(a, n);
Chris@82 772 }
Chris@82 773
Chris@82 774 static void mkio10(C *a, int n)
Chris@82 775 {
Chris@82 776 mkoddonly(a, n);
Chris@82 777 mkio00(a, n);
Chris@82 778 }
Chris@82 779
Chris@82 780 static void mkre11(C *a, int n)
Chris@82 781 {
Chris@82 782 mkoddonly(a, n);
Chris@82 783 mko00(a, n/2, 0);
Chris@82 784 mkre00(a, n);
Chris@82 785 }
Chris@82 786
Chris@82 787 static void mkro11(C *a, int n)
Chris@82 788 {
Chris@82 789 mkoddonly(a, n);
Chris@82 790 mkre00(a, n/2);
Chris@82 791 mkro00(a, n);
Chris@82 792 }
Chris@82 793
Chris@82 794 static void mkio11(C *a, int n)
Chris@82 795 {
Chris@82 796 mkoddonly(a, n);
Chris@82 797 mke00(a, n/2, 1);
Chris@82 798 mkio00(a, n);
Chris@82 799 }
Chris@82 800
Chris@82 801 static void r2r_apply(dofft_closure *k_, bench_complex *in, bench_complex *out)
Chris@82 802 {
Chris@82 803 dofft_r2r_closure *k = (dofft_r2r_closure *)k_;
Chris@82 804 bench_problem *p = k->p;
Chris@82 805 bench_real *ri, *ro;
Chris@82 806 int n, is, os;
Chris@82 807
Chris@82 808 n = p->sz->dims[0].n;
Chris@82 809 is = p->sz->dims[0].is;
Chris@82 810 os = p->sz->dims[0].os;
Chris@82 811
Chris@82 812 ri = (bench_real *) p->in;
Chris@82 813 ro = (bench_real *) p->out;
Chris@82 814
Chris@82 815 switch (p->k[0]) {
Chris@82 816 case R2R_R2HC:
Chris@82 817 cpyr1(n, &c_re(in[0]), 2, ri, is, 1.0);
Chris@82 818 break;
Chris@82 819 case R2R_HC2R:
Chris@82 820 cpyr1(n/2 + 1, &c_re(in[0]), 2, ri, is, 1.0);
Chris@82 821 cpyr1((n+1)/2 - 1, &c_im(in[n-1]), -2, ri + is*(n-1), -is, 1.0);
Chris@82 822 break;
Chris@82 823 case R2R_REDFT00:
Chris@82 824 cpyr1(n, &c_re(in[0]), 2, ri, is, 1.0);
Chris@82 825 break;
Chris@82 826 case R2R_RODFT00:
Chris@82 827 cpyr1(n, &c_re(in[1]), 2, ri, is, 1.0);
Chris@82 828 break;
Chris@82 829 case R2R_REDFT01:
Chris@82 830 cpyr1(n, &c_re(in[0]), 2, ri, is, 1.0);
Chris@82 831 break;
Chris@82 832 case R2R_REDFT10:
Chris@82 833 cpyr1(n, &c_re(in[1]), 4, ri, is, 1.0);
Chris@82 834 break;
Chris@82 835 case R2R_RODFT01:
Chris@82 836 cpyr1(n, &c_re(in[1]), 2, ri, is, 1.0);
Chris@82 837 break;
Chris@82 838 case R2R_RODFT10:
Chris@82 839 cpyr1(n, &c_im(in[1]), 4, ri, is, 1.0);
Chris@82 840 break;
Chris@82 841 case R2R_REDFT11:
Chris@82 842 cpyr1(n, &c_re(in[1]), 4, ri, is, 1.0);
Chris@82 843 break;
Chris@82 844 case R2R_RODFT11:
Chris@82 845 cpyr1(n, &c_re(in[1]), 4, ri, is, 1.0);
Chris@82 846 break;
Chris@82 847 default:
Chris@82 848 BENCH_ASSERT(0); /* not yet implemented */
Chris@82 849 }
Chris@82 850
Chris@82 851 after_problem_rcopy_from(p, ri);
Chris@82 852 doit(1, p);
Chris@82 853 after_problem_rcopy_to(p, ro);
Chris@82 854
Chris@82 855 switch (p->k[0]) {
Chris@82 856 case R2R_R2HC:
Chris@82 857 if (k->k.recopy_input)
Chris@82 858 cpyr1(n, ri, is, &c_re(in[0]), 2, 1.0);
Chris@82 859 cpyr1(n/2 + 1, ro, os, &c_re(out[0]), 2, 1.0);
Chris@82 860 cpyr1((n+1)/2 - 1, ro + os*(n-1), -os, &c_im(out[1]), 2, 1.0);
Chris@82 861 c_im(out[0]) = 0.0;
Chris@82 862 if (n % 2 == 0)
Chris@82 863 c_im(out[n/2]) = 0.0;
Chris@82 864 mkhermitian1(out, n);
Chris@82 865 break;
Chris@82 866 case R2R_HC2R:
Chris@82 867 if (k->k.recopy_input) {
Chris@82 868 cpyr1(n/2 + 1, ri, is, &c_re(in[0]), 2, 1.0);
Chris@82 869 cpyr1((n+1)/2 - 1, ri + is*(n-1), -is, &c_im(in[1]), 2,1.0);
Chris@82 870 }
Chris@82 871 cpyr1(n, ro, os, &c_re(out[0]), 2, 1.0);
Chris@82 872 mkreal(out, n);
Chris@82 873 break;
Chris@82 874 case R2R_REDFT00:
Chris@82 875 if (k->k.recopy_input)
Chris@82 876 cpyr1(n, ri, is, &c_re(in[0]), 2, 1.0);
Chris@82 877 cpyr1(n, ro, os, &c_re(out[0]), 2, 1.0);
Chris@82 878 mkre00(out, k->n0);
Chris@82 879 break;
Chris@82 880 case R2R_RODFT00:
Chris@82 881 if (k->k.recopy_input)
Chris@82 882 cpyr1(n, ri, is, &c_im(in[1]), 2, -1.0);
Chris@82 883 cpyr1(n, ro, os, &c_im(out[1]), 2, -1.0);
Chris@82 884 mkio00(out, k->n0);
Chris@82 885 break;
Chris@82 886 case R2R_REDFT01:
Chris@82 887 if (k->k.recopy_input)
Chris@82 888 cpyr1(n, ri, is, &c_re(in[0]), 2, 1.0);
Chris@82 889 cpyr1(n, ro, os, &c_re(out[1]), 4, 2.0);
Chris@82 890 mkre10(out, k->n0);
Chris@82 891 break;
Chris@82 892 case R2R_REDFT10:
Chris@82 893 if (k->k.recopy_input)
Chris@82 894 cpyr1(n, ri, is, &c_re(in[1]), 4, 2.0);
Chris@82 895 cpyr1(n, ro, os, &c_re(out[0]), 2, 1.0);
Chris@82 896 mkre01(out, k->n0);
Chris@82 897 break;
Chris@82 898 case R2R_RODFT01:
Chris@82 899 if (k->k.recopy_input)
Chris@82 900 cpyr1(n, ri, is, &c_re(in[1]), 2, 1.0);
Chris@82 901 cpyr1(n, ro, os, &c_im(out[1]), 4, -2.0);
Chris@82 902 mkio10(out, k->n0);
Chris@82 903 break;
Chris@82 904 case R2R_RODFT10:
Chris@82 905 if (k->k.recopy_input)
Chris@82 906 cpyr1(n, ri, is, &c_im(in[1]), 4, -2.0);
Chris@82 907 cpyr1(n, ro, os, &c_re(out[1]), 2, 1.0);
Chris@82 908 mkro01(out, k->n0);
Chris@82 909 break;
Chris@82 910 case R2R_REDFT11:
Chris@82 911 if (k->k.recopy_input)
Chris@82 912 cpyr1(n, ri, is, &c_re(in[1]), 4, 2.0);
Chris@82 913 cpyr1(n, ro, os, &c_re(out[1]), 4, 2.0);
Chris@82 914 mkre11(out, k->n0);
Chris@82 915 break;
Chris@82 916 case R2R_RODFT11:
Chris@82 917 if (k->k.recopy_input)
Chris@82 918 cpyr1(n, ri, is, &c_im(in[1]), 4, -2.0);
Chris@82 919 cpyr1(n, ro, os, &c_im(out[1]), 4, -2.0);
Chris@82 920 mkio11(out, k->n0);
Chris@82 921 break;
Chris@82 922 default:
Chris@82 923 BENCH_ASSERT(0); /* not yet implemented */
Chris@82 924 }
Chris@82 925 }
Chris@82 926
Chris@82 927 void accuracy_r2r(bench_problem *p, int rounds, int impulse_rounds,
Chris@82 928 double t[6])
Chris@82 929 {
Chris@82 930 dofft_r2r_closure k;
Chris@82 931 int n, n0 = 1;
Chris@82 932 C *a, *b;
Chris@82 933 aconstrain constrain = 0;
Chris@82 934
Chris@82 935 BENCH_ASSERT(p->kind == PROBLEM_R2R);
Chris@82 936 BENCH_ASSERT(p->sz->rnk == 1);
Chris@82 937 BENCH_ASSERT(p->vecsz->rnk == 0);
Chris@82 938
Chris@82 939 k.k.apply = r2r_apply;
Chris@82 940 k.k.recopy_input = 0;
Chris@82 941 k.p = p;
Chris@82 942 n = tensor_sz(p->sz);
Chris@82 943
Chris@82 944 switch (p->k[0]) {
Chris@82 945 case R2R_R2HC: constrain = mkreal; n0 = n; break;
Chris@82 946 case R2R_HC2R: constrain = mkhermitian1; n0 = n; break;
Chris@82 947 case R2R_REDFT00: constrain = mkre00; n0 = 2*(n-1); break;
Chris@82 948 case R2R_RODFT00: constrain = mkro00; n0 = 2*(n+1); break;
Chris@82 949 case R2R_REDFT01: constrain = mkre01; n0 = 4*n; break;
Chris@82 950 case R2R_REDFT10: constrain = mkre10; n0 = 4*n; break;
Chris@82 951 case R2R_RODFT01: constrain = mkro01; n0 = 4*n; break;
Chris@82 952 case R2R_RODFT10: constrain = mkio10; n0 = 4*n; break;
Chris@82 953 case R2R_REDFT11: constrain = mkre11; n0 = 8*n; break;
Chris@82 954 case R2R_RODFT11: constrain = mkro11; n0 = 8*n; break;
Chris@82 955 default: BENCH_ASSERT(0); /* not yet implemented */
Chris@82 956 }
Chris@82 957 k.n0 = n0;
Chris@82 958
Chris@82 959 a = (C *) bench_malloc(n0 * sizeof(C));
Chris@82 960 b = (C *) bench_malloc(n0 * sizeof(C));
Chris@82 961 accuracy_test(&k.k, constrain, -1, n0, a, b, rounds, impulse_rounds, t);
Chris@82 962 bench_free(b);
Chris@82 963 bench_free(a);
Chris@82 964 }