Mercurial > hg > sv-dependency-builds

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