WebSVN - shark - Blame - Rev 2 - /shark/branches/pj/ports/fftw/tests/rfftw

Rev	Author	Line No.	Line
2	pj	1	/*
		2	* Copyright (c) 1997-1999 Massachusetts Institute of Technology
		3	*
		4	* This program is free software; you can redistribute it and/or modify
		5	* it under the terms of the GNU General Public License as published by
		6	* the Free Software Foundation; either version 2 of the License, or
		7	* (at your option) any later version.
		8	*
		9	* This program is distributed in the hope that it will be useful,
		10	* but WITHOUT ANY WARRANTY; without even the implied warranty of
		11	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
		12	* GNU General Public License for more details.
		13	*
		14	* You should have received a copy of the GNU General Public License
		15	* along with this program; if not, write to the Free Software
		16	* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
		17	*
		18	*/
		19
		20	/*
		21	* fftw_test.c : test program for complex-complex transforms
		22	*/
		23
		24	/* $Id: rfftw_te.c,v 1.1.1.1 2002-03-29 14:12:59 pj Exp $ */
		25	#include <fftw-int.h>
		26	#include <stdlib.h>
		27	#include <stdio.h>
		28	#include <string.h>
		29	#include <math.h>
		30	#include <time.h>
		31
		32	#include <rfftw.h>
		33
		34	#include "test_main.h"
		35
		36	char fftw_prefix[] = "rfftw";
		37
		38	/*************************************************
		39	* Speed tests
		40	*************************************************/
		41
		42	void zero_arr(int n, fftw_real * a)
		43	{
		44	int i;
		45	for (i = 0; i < n; ++i)
		46	a[i] = 0.0;
		47	}
		48
		49	void test_speed_aux(int n, fftw_direction dir, int flags, int specific)
		50	{
		51	fftw_real in, out;
		52	fftw_plan plan;
		53	double t;
		54	fftw_time begin, end;
		55
		56	in = (fftw_real ) fftw_malloc(n howmany_fields
		57	* sizeof(fftw_real));
		58	out = (fftw_real ) fftw_malloc(n howmany_fields
		59	* sizeof(fftw_real));
		60
		61	if (specific) {
		62	begin = fftw_get_time();
		63	plan = rfftw_create_plan_specific(n, dir,
		64	speed_flag \| flags \| wisdom_flag,
		65	in, howmany_fields,
		66	out, howmany_fields);
		67	end = fftw_get_time();
		68	} else {
		69	begin = fftw_get_time();
		70	plan = rfftw_create_plan(n, dir, speed_flag \| flags \| wisdom_flag);
		71	end = fftw_get_time();
		72	}
		73	CHECK(plan != NULL, "can't create plan");
		74
		75	t = fftw_time_to_sec(fftw_time_diff(end, begin));
		76	WHEN_VERBOSE(2, printf("time for planner: %f s\n", t));
		77
		78	WHEN_VERBOSE(2, rfftw_print_plan(plan));
		79
		80	FFTW_TIME_FFT(rfftw(plan, howmany_fields,
		81	in, howmany_fields, 1, out, howmany_fields, 1),
		82	in, n * howmany_fields, t);
		83
		84	rfftw_destroy_plan(plan);
		85
		86	WHEN_VERBOSE(1, printf("time for one fft: %s", smart_sprint_time(t)));
		87	WHEN_VERBOSE(1, printf(" (%s/point)\n", smart_sprint_time(t / n)));
		88	WHEN_VERBOSE(1, printf("\"mflops\" = 5/2 (n log2 n) / (t in microseconds)"
		89	" = %f\n", 0.5 * howmany_fields * mflops(t, n)));
		90
		91	fftw_free(in);
		92	fftw_free(out);
		93
		94	WHEN_VERBOSE(1, printf("\n"));
		95	}
		96
		97	void test_speed_nd_aux(struct size sz,
		98	fftw_direction dir, int flags, int specific)
		99	{
		100	fftw_real *in;
		101	fftwnd_plan plan;
		102	double t;
		103	fftw_time begin, end;
		104	int i, N;
		105
		106	/* only bench in-place multi-dim transforms */
		107	flags \|= FFTW_IN_PLACE;
		108
		109	N = 1;
		110	for (i = 0; i < sz.rank - 1; ++i)
		111	N *= sz.narray[i];
		112
		113	N *= (sz.narray[i] + 2);
		114
		115	in = (fftw_real ) fftw_malloc(N howmany_fields * sizeof(fftw_real));
		116
		117	if (specific) {
		118	begin = fftw_get_time();
		119	plan = rfftwnd_create_plan_specific(sz.rank, sz.narray, dir,
		120	speed_flag \| flags \| wisdom_flag,
		121	in, howmany_fields, 0, 1);
		122	} else {
		123	begin = fftw_get_time();
		124	plan = rfftwnd_create_plan(sz.rank, sz.narray,
		125	dir, speed_flag \| flags \| wisdom_flag);
		126	}
		127	end = fftw_get_time();
		128	CHECK(plan != NULL, "can't create plan");
		129
		130	t = fftw_time_to_sec(fftw_time_diff(end, begin));
		131	WHEN_VERBOSE(2, printf("time for planner: %f s\n", t));
		132
		133	WHEN_VERBOSE(2, printf("\n"));
		134	WHEN_VERBOSE(2, (rfftwnd_print_plan(plan)));
		135	WHEN_VERBOSE(2, printf("\n"));
		136
		137	if (dir == FFTW_REAL_TO_COMPLEX) {
		138	FFTW_TIME_FFT(rfftwnd_real_to_complex(plan, howmany_fields,
		139	in, howmany_fields, 1,
		140	0, 0, 0),
		141	in, N * howmany_fields, t);
		142	} else {
		143	FFTW_TIME_FFT(rfftwnd_complex_to_real(plan, howmany_fields,
		144	(fftw_complex *) in,
		145	howmany_fields, 1,
		146	0, 0, 0),
		147	in, N * howmany_fields, t);
		148	}
		149
		150	rfftwnd_destroy_plan(plan);
		151
		152	WHEN_VERBOSE(1, printf("time for one fft: %s", smart_sprint_time(t)));
		153	WHEN_VERBOSE(1, printf(" (%s/point)\n", smart_sprint_time(t / N)));
		154	WHEN_VERBOSE(1, printf("\"mflops\" = 5/2 (N log2 N) / (t in microseconds)"
		155	" = %f\n", 0.5 * howmany_fields * mflops(t, N)));
		156
		157	fftw_free(in);
		158
		159	WHEN_VERBOSE(1, printf("\n"));
		160	}
		161
		162	/*************************************************
		163	* correctness tests
		164	*************************************************/
		165
		166	void fill_random(fftw_real * a, int n, int stride)
		167	{
		168	int i;
		169
		170	/* generate random inputs */
		171	for (i = 0; i < n; ++i)
		172	a[i * stride] = DRAND();
		173	}
		174
		175	double compute_error(fftw_real * A, int astride,
		176	fftw_real * B, int bstride, int n)
		177	{
		178	/* compute the relative error */
		179	double error = 0.0;
		180	int i;
		181
		182	for (i = 0; i < n; ++i) {
		183	double a;
		184	double mag;
		185	a = fabs(A[i * astride] - B[i * bstride]);
		186	mag = 0.5 * (fabs(A[i * astride]) + fabs(B[i * bstride])) + TOLERANCE;
		187
		188	a /= mag;
		189	if (a > error)
		190	error = a;
		191
		192	#ifdef HAVE_ISNAN
		193	CHECK(!isnan(a), "NaN in answer");
		194	#endif
		195	}
		196	return error;
		197	}
		198
		199	void array_compare(fftw_real * A, fftw_real * B, int n)
		200	{
		201	CHECK(compute_error(A, 1, B, 1, n) < TOLERANCE,
		202	"failure in RFFTW verification");
		203	}
		204
		205	void test_out_of_place(int n, int istride, int ostride,
		206	int howmany, fftw_direction dir,
		207	fftw_plan validated_plan, int specific)
		208	{
		209	fftw_complex in2, out2;
		210	fftw_real in1, out1, *out3;
		211	fftw_plan plan;
		212	int i, j;
		213	int flags = measure_flag \| wisdom_flag;
		214
		215	if (coinflip())
		216	flags \|= FFTW_THREADSAFE;
		217
		218	in1 = (fftw_real ) fftw_malloc(istride n * sizeof(fftw_real) * howmany);
		219	in2 = (fftw_complex ) fftw_malloc(n sizeof(fftw_complex));
		220	out1 = (fftw_real ) fftw_malloc(ostride n * sizeof(fftw_real) * howmany);
		221	out2 = (fftw_complex ) fftw_malloc(n sizeof(fftw_complex));
		222	out3 = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		223
		224	if (!specific)
		225	plan = rfftw_create_plan(n, dir, flags);
		226	else
		227	plan = rfftw_create_plan_specific(n, dir, flags,
		228	in1, istride, out1, ostride);
		229	CHECK(plan != NULL, "can't create plan");
		230
		231	/* generate random inputs */
		232	fill_random(in1, n, istride);
		233	for (j = 1; j < howmany; ++j)
		234	for (i = 0; i < n; ++i)
		235	in1[(j * n + i) * istride] = in1[i * istride];
		236
		237	/* copy random inputs to complex array for comparison with fftw: */
		238	if (dir == FFTW_REAL_TO_COMPLEX)
		239	for (i = 0; i < n; ++i) {
		240	c_re(in2[i]) = in1[i * istride];
		241	c_im(in2[i]) = 0.0;
		242	} else {
		243	int n2 = (n + 1) / 2;
		244	c_re(in2[0]) = in1[0];
		245	c_im(in2[0]) = 0.0;
		246	for (i = 1; i < n2; ++i) {
		247	c_re(in2[i]) = in1[i * istride];
		248	c_im(in2[i]) = in1[(n - i) * istride];
		249	}
		250	if (n2 * 2 == n) {
		251	c_re(in2[n2]) = in1[n2 * istride];
		252	c_im(in2[n2]) = 0.0;
		253	++i;
		254	}
		255	for (; i < n; ++i) {
		256	c_re(in2[i]) = c_re(in2[n - i]);
		257	c_im(in2[i]) = -c_im(in2[n - i]);
		258	}
		259	}
		260
		261	/*
		262	* fill in other positions of the array, to make sure that
		263	* rfftw doesn't overwrite them
		264	*/
		265	for (j = 1; j < istride; ++j)
		266	for (i = 0; i < n * howmany; ++i)
		267	in1[i * istride + j] = i * istride + j;
		268
		269	for (j = 1; j < ostride; ++j)
		270	for (i = 0; i < n * howmany; ++i)
		271	out1[i * ostride + j] = -i * ostride + j;
		272
		273	WHEN_VERBOSE(2, rfftw_print_plan(plan));
		274
		275	/* fft-ize */
		276	if (howmany != 1 \|\| istride != 1 \|\| ostride != 1 \|\| coinflip())
		277	rfftw(plan, howmany, in1, istride, n * istride, out1, ostride,
		278	n * ostride);
		279	else
		280	rfftw_one(plan, in1, out1);
		281
		282	rfftw_destroy_plan(plan);
		283
		284	/* check for overwriting */
		285	for (j = 1; j < istride; ++j)
		286	for (i = 0; i < n * howmany; ++i)
		287	CHECK(in1[i * istride + j] == i * istride + j,
		288	"input has been overwritten");
		289	for (j = 1; j < ostride; ++j)
		290	for (i = 0; i < n * howmany; ++i)
		291	CHECK(out1[i * ostride + j] == -i * ostride + j,
		292	"output has been overwritten");
		293
		294	fftw(validated_plan, 1, in2, 1, n, out2, 1, n);
		295
		296	if (dir == FFTW_REAL_TO_COMPLEX) {
		297	int n2 = (n + 1) / 2;
		298	out3[0] = c_re(out2[0]);
		299	for (i = 1; i < n2; ++i) {
		300	out3[i] = c_re(out2[i]);
		301	out3[n - i] = c_im(out2[i]);
		302	}
		303	if (n2 * 2 == n)
		304	out3[n2] = c_re(out2[n2]);
		305	} else {
		306	for (i = 0; i < n; ++i)
		307	out3[i] = c_re(out2[i]);
		308	}
		309
		310	for (j = 0; j < howmany; ++j)
		311	CHECK(compute_error(out1 + j * n * ostride, ostride, out3, 1, n)
		312	< TOLERANCE,
		313	"test_out_of_place: wrong answer");
		314	WHEN_VERBOSE(2, printf("OK\n"));
		315
		316	fftw_free(in1);
		317	fftw_free(in2);
		318	fftw_free(out1);
		319	fftw_free(out2);
		320	fftw_free(out3);
		321	}
		322
		323	void test_in_place(int n, int istride,
		324	int howmany, fftw_direction dir,
		325	fftw_plan validated_plan, int specific)
		326	{
		327	fftw_complex in2, out2;
		328	fftw_real in1, out1, *out3;
		329	fftw_plan plan;
		330	int i, j;
		331	int ostride = istride;
		332	int flags = measure_flag \| wisdom_flag \| FFTW_IN_PLACE;
		333
		334	if (coinflip())
		335	flags \|= FFTW_THREADSAFE;
		336
		337	in1 = (fftw_real ) fftw_malloc(istride n * sizeof(fftw_real) * howmany);
		338	in2 = (fftw_complex ) fftw_malloc(n sizeof(fftw_complex));
		339	out1 = in1;
		340	out2 = (fftw_complex ) fftw_malloc(n sizeof(fftw_complex));
		341	out3 = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		342
		343	if (!specific)
		344	plan = rfftw_create_plan(n, dir, flags);
		345	else
		346	plan = rfftw_create_plan_specific(n, dir, flags,
		347	in1, istride, out1, ostride);
		348	CHECK(plan != NULL, "can't create plan");
		349
		350	/* generate random inputs */
		351	fill_random(in1, n, istride);
		352	for (j = 1; j < howmany; ++j)
		353	for (i = 0; i < n; ++i)
		354	in1[(j * n + i) * istride] = in1[i * istride];
		355
		356	/* copy random inputs to complex array for comparison with fftw: */
		357	if (dir == FFTW_REAL_TO_COMPLEX)
		358	for (i = 0; i < n; ++i) {
		359	c_re(in2[i]) = in1[i * istride];
		360	c_im(in2[i]) = 0.0;
		361	} else {
		362	int n2 = (n + 1) / 2;
		363	c_re(in2[0]) = in1[0];
		364	c_im(in2[0]) = 0.0;
		365	for (i = 1; i < n2; ++i) {
		366	c_re(in2[i]) = in1[i * istride];
		367	c_im(in2[i]) = in1[(n - i) * istride];
		368	}
		369	if (n2 * 2 == n) {
		370	c_re(in2[n2]) = in1[n2 * istride];
		371	c_im(in2[n2]) = 0.0;
		372	++i;
		373	}
		374	for (; i < n; ++i) {
		375	c_re(in2[i]) = c_re(in2[n - i]);
		376	c_im(in2[i]) = -c_im(in2[n - i]);
		377	}
		378	}
		379
		380	/*
		381	* fill in other positions of the array, to make sure that
		382	* rfftw doesn't overwrite them
		383	*/
		384	for (j = 1; j < istride; ++j)
		385	for (i = 0; i < n * howmany; ++i)
		386	in1[i * istride + j] = i * istride + j;
		387
		388	WHEN_VERBOSE(2, rfftw_print_plan(plan));
		389
		390	/* fft-ize */
		391	if (howmany != 1 \|\| istride != 1 \|\| coinflip())
		392	rfftw(plan, howmany, in1, istride, n * istride, 0, 0, 0);
		393	else
		394	rfftw_one(plan, in1, NULL);
		395
		396	rfftw_destroy_plan(plan);
		397
		398	/* check for overwriting */
		399	for (j = 1; j < ostride; ++j)
		400	for (i = 0; i < n * howmany; ++i)
		401	CHECK(out1[i * ostride + j] == i * ostride + j,
		402	"output has been overwritten");
		403
		404	fftw(validated_plan, 1, in2, 1, n, out2, 1, n);
		405
		406	if (dir == FFTW_REAL_TO_COMPLEX) {
		407	int n2 = (n + 1) / 2;
		408	out3[0] = c_re(out2[0]);
		409	for (i = 1; i < n2; ++i) {
		410	out3[i] = c_re(out2[i]);
		411	out3[n - i] = c_im(out2[i]);
		412	}
		413	if (n2 * 2 == n)
		414	out3[n2] = c_re(out2[n2]);
		415	} else {
		416	for (i = 0; i < n; ++i)
		417	out3[i] = c_re(out2[i]);
		418	}
		419
		420	for (j = 0; j < howmany; ++j)
		421	CHECK(compute_error(out1 + j * n * ostride, ostride, out3, 1, n)
		422	< TOLERANCE,
		423	"test_in_place: wrong answer");
		424	WHEN_VERBOSE(2, printf("OK\n"));
		425
		426	fftw_free(in1);
		427	fftw_free(in2);
		428	fftw_free(out2);
		429	fftw_free(out3);
		430	}
		431
		432	void test_out_of_place_both(int n, int istride, int ostride,
		433	int howmany,
		434	fftw_plan validated_plan_forward,
		435	fftw_plan validated_plan_backward)
		436	{
		437	int specific;
		438
		439	for (specific = 0; specific <= 1; ++specific) {
		440	WHEN_VERBOSE(2,
		441	printf("TEST CORRECTNESS (out of place, FFTW_FORWARD, %s)"
		442	" n = %d istride = %d ostride = %d howmany = %d\n",
		443	SPECIFICP(specific),
		444	n, istride, ostride, howmany));
		445	test_out_of_place(n, istride, ostride, howmany, FFTW_FORWARD,
		446	validated_plan_forward, specific);
		447
		448	WHEN_VERBOSE(2,
		449	printf("TEST CORRECTNESS (out of place, FFTW_BACKWARD, %s)"
		450	" n = %d istride = %d ostride = %d howmany = %d\n",
		451	SPECIFICP(specific),
		452	n, istride, ostride, howmany));
		453	test_out_of_place(n, istride, ostride, howmany, FFTW_BACKWARD,
		454	validated_plan_backward, specific);
		455	}
		456	}
		457
		458	void test_in_place_both(int n, int istride, int howmany,
		459	fftw_plan validated_plan_forward,
		460	fftw_plan validated_plan_backward)
		461	{
		462	int specific;
		463
		464	for (specific = 0; specific <= 1; ++specific) {
		465	WHEN_VERBOSE(2,
		466	printf("TEST CORRECTNESS (in place, FFTW_FORWARD, %s) "
		467	"n = %d istride = %d howmany = %d\n",
		468	SPECIFICP(specific),
		469	n, istride, howmany));
		470	test_in_place(n, istride, howmany, FFTW_FORWARD,
		471	validated_plan_forward, specific);
		472
		473	WHEN_VERBOSE(2,
		474	printf("TEST CORRECTNESS (in place, FFTW_BACKWARD, %s) "
		475	"n = %d istride = %d howmany = %d\n",
		476	SPECIFICP(specific),
		477	n, istride, howmany));
		478	test_in_place(n, istride, howmany, FFTW_BACKWARD,
		479	validated_plan_backward, specific);
		480	}
		481	}
		482
		483	void test_correctness(int n)
		484	{
		485	int istride, ostride, howmany;
		486	fftw_plan validated_plan_forward, validated_plan_backward;
		487
		488	WHEN_VERBOSE(1,
		489	printf("Testing correctness for n = %d...", n);
		490	fflush(stdout));
		491
		492	/* produce a good plan (validated by Ergun's test procedure) */
		493	validated_plan_forward =
		494	fftw_create_plan(n, FFTW_FORWARD, measure_flag \| wisdom_flag);
		495	validated_plan_backward =
		496	fftw_create_plan(n, FFTW_BACKWARD, measure_flag \| wisdom_flag);
		497	CHECK(validated_plan_forward != NULL, "can't create plan");
		498	CHECK(validated_plan_backward != NULL, "can't create plan");
		499
		500	for (istride = 1; istride <= MAX_STRIDE; ++istride)
		501	for (ostride = 1; ostride <= MAX_STRIDE; ++ostride)
		502	for (howmany = 1; howmany <= MAX_HOWMANY; ++howmany)
		503	test_out_of_place_both(n, istride, ostride, howmany,
		504	validated_plan_forward,
		505	validated_plan_backward);
		506
		507	for (istride = 1; istride <= MAX_STRIDE; ++istride)
		508	for (howmany = 1; howmany <= MAX_HOWMANY; ++howmany)
		509	test_in_place_both(n, istride, howmany,
		510	validated_plan_forward,
		511	validated_plan_backward);
		512
		513	fftw_destroy_plan(validated_plan_forward);
		514	fftw_destroy_plan(validated_plan_backward);
		515
		516	if (!(wisdom_flag & FFTW_USE_WISDOM) && chk_mem_leak)
		517	fftw_check_memory_leaks();
		518
		519	WHEN_VERBOSE(1, printf("OK\n"));
		520	}
		521
		522	/*************************************************
		523	* multi-dimensional correctness tests
		524	*************************************************/
		525
		526	void testnd_out_of_place(int rank, int *n, fftwnd_plan validated_plan)
		527	{
		528	int istride, ostride;
		529	int N, dim, i, j, k;
		530	int nc, nhc, nr;
		531	fftw_real in1, out3;
		532	fftw_complex in2, out1, *out2;
		533	fftwnd_plan p, ip;
		534	int flags = measure_flag \| wisdom_flag;
		535
		536	if (coinflip())
		537	flags \|= FFTW_THREADSAFE;
		538
		539	N = nc = nr = nhc = 1;
		540	for (dim = 0; dim < rank; ++dim)
		541	N *= n[dim];
		542	if (rank > 0) {
		543	nr = n[rank - 1];
		544	nc = N / nr;
		545	nhc = nr / 2 + 1;
		546	}
		547	in1 = (fftw_real ) fftw_malloc(N MAX_STRIDE * sizeof(fftw_real));
		548	out3 = (fftw_real ) fftw_malloc(N MAX_STRIDE * sizeof(fftw_real));
		549	out1 = (fftw_complex ) fftw_malloc(nhc nc * MAX_STRIDE
		550	* sizeof(fftw_complex));
		551	in2 = (fftw_complex ) fftw_malloc(N sizeof(fftw_complex));
		552	out2 = (fftw_complex ) fftw_malloc(N sizeof(fftw_complex));
		553
		554	p = rfftwnd_create_plan(rank, n, FFTW_REAL_TO_COMPLEX, flags);
		555	ip = rfftwnd_create_plan(rank, n, FFTW_COMPLEX_TO_REAL, flags);
		556	CHECK(p != NULL && ip != NULL, "can't create plan");
		557
		558	for (istride = 1; istride <= MAX_STRIDE; ++istride) {
		559	/* generate random inputs */
		560	for (i = 0; i < nc; ++i)
		561	for (j = 0; j < nr; ++j) {
		562	c_re(in2[i * nr + j]) = DRAND();
		563	c_im(in2[i * nr + j]) = 0.0;
		564	for (k = 0; k < istride; ++k)
		565	in1[(i * nr + j) * istride + k]
		566	= c_re(in2[i * nr + j]);
		567	}
		568	for (i = 0; i < N * istride; ++i)
		569	out3[i] = 0.0;
		570
		571	fftwnd(validated_plan, 1, in2, 1, 1, out2, 1, 1);
		572
		573	for (ostride = 1; ostride <= MAX_STRIDE; ++ostride) {
		574	int howmany = (istride < ostride) ? istride : ostride;
		575
		576	WHEN_VERBOSE(2, printf("\n testing stride %d/%d...",
		577	istride, ostride));
		578
		579	if (howmany != 1 \|\| istride != 1 \|\| ostride != 1 \|\| coinflip())
		580	rfftwnd_real_to_complex(p, howmany, in1, istride, 1,
		581	out1, ostride, 1);
		582	else
		583	rfftwnd_one_real_to_complex(p, in1, out1);
		584
		585	for (i = 0; i < nc; ++i)
		586	for (k = 0; k < howmany; ++k)
		587	CHECK(compute_error_complex(out1 + i * nhc * ostride + k,
		588	ostride,
		589	out2 + i * nr, 1,
		590	nhc) < TOLERANCE,
		591	"out-of-place (r2c): wrong answer");
		592
		593	if (howmany != 1 \|\| istride != 1 \|\| ostride != 1 \|\| coinflip())
		594	rfftwnd_complex_to_real(ip, howmany, out1, ostride, 1,
		595	out3, istride, 1);
		596	else
		597	rfftwnd_one_complex_to_real(ip, out1, out3);
		598
		599	for (i = 0; i < N * istride; ++i)
		600	out3[i] *= 1.0 / N;
		601
		602	if (istride == howmany)
		603	CHECK(compute_error(out3, 1, in1, 1, N * istride)
		604	< TOLERANCE, "out-of-place (c2r): wrong answer");
		605	for (i = 0; i < nc; ++i)
		606	for (k = 0; k < howmany; ++k)
		607	CHECK(compute_error(out3 + i * nr * istride + k,
		608	istride,
		609	(fftw_real ) (in2 + i nr), 2,
		610	nr) < TOLERANCE,
		611	"out-of-place (c2r): wrong answer (check 2)");
		612	}
		613	}
		614
		615	rfftwnd_destroy_plan(p);
		616	rfftwnd_destroy_plan(ip);
		617
		618	fftw_free(out3);
		619	fftw_free(out2);
		620	fftw_free(in2);
		621	fftw_free(out1);
		622	fftw_free(in1);
		623	}
		624
		625	void testnd_in_place(int rank, int *n, fftwnd_plan validated_plan,
		626	int alternate_api, int specific)
		627	{
		628	int istride, ostride, howmany;
		629	int N, dim, i, j, k;
		630	int nc, nhc, nr;
		631	fftw_real in1, out3;
		632	fftw_complex in2, out1, *out2;
		633	fftwnd_plan p, ip;
		634	int flags = measure_flag \| wisdom_flag \| FFTW_IN_PLACE;
		635
		636	if (coinflip())
		637	flags \|= FFTW_THREADSAFE;
		638
		639	N = nc = nr = nhc = 1;
		640	for (dim = 0; dim < rank; ++dim)
		641	N *= n[dim];
		642	if (rank > 0) {
		643	nr = n[rank - 1];
		644	nc = N / nr;
		645	nhc = nr / 2 + 1;
		646	}
		647	in1 = (fftw_real ) fftw_malloc(2 nhc * nc * MAX_STRIDE * sizeof(fftw_real));
		648	out3 = in1;
		649	out1 = (fftw_complex *) in1;
		650	in2 = (fftw_complex ) fftw_malloc(N sizeof(fftw_complex));
		651	out2 = (fftw_complex ) fftw_malloc(N sizeof(fftw_complex));
		652
		653	if (alternate_api && specific && (rank == 2 \|\| rank == 3)) {
		654	if (rank == 2) {
		655	p = rfftw2d_create_plan_specific(n[0], n[1],
		656	FFTW_REAL_TO_COMPLEX, flags,
		657	in1, MAX_STRIDE, 0, 0);
		658	ip = rfftw2d_create_plan_specific(n[0], n[1],
		659	FFTW_COMPLEX_TO_REAL, flags,
		660	in1, MAX_STRIDE, 0, 0);
		661	} else {
		662	p = rfftw3d_create_plan_specific(n[0], n[1], n[2],
		663	FFTW_REAL_TO_COMPLEX, flags,
		664	in1, MAX_STRIDE, 0, 0);
		665	ip = rfftw3d_create_plan_specific(n[0], n[1], n[2],
		666	FFTW_COMPLEX_TO_REAL, flags,
		667	in1, MAX_STRIDE, 0, 0);
		668	}
		669	} else if (specific) {
		670	p = rfftwnd_create_plan_specific(rank, n, FFTW_REAL_TO_COMPLEX,
		671	flags,
		672	in1, MAX_STRIDE, in1, MAX_STRIDE);
		673	ip = rfftwnd_create_plan_specific(rank, n, FFTW_COMPLEX_TO_REAL,
		674	flags,
		675	in1, MAX_STRIDE, in1, MAX_STRIDE);
		676	} else if (alternate_api && (rank == 2 \|\| rank == 3)) {
		677	if (rank == 2) {
		678	p = rfftw2d_create_plan(n[0], n[1], FFTW_REAL_TO_COMPLEX,
		679	flags);
		680	ip = rfftw2d_create_plan(n[0], n[1], FFTW_COMPLEX_TO_REAL,
		681	flags);
		682	} else {
		683	p = rfftw3d_create_plan(n[0], n[1], n[2], FFTW_REAL_TO_COMPLEX,
		684	flags);
		685	ip = rfftw3d_create_plan(n[0], n[1], n[2], FFTW_COMPLEX_TO_REAL,
		686	flags);
		687	}
		688	} else {
		689	p = rfftwnd_create_plan(rank, n, FFTW_REAL_TO_COMPLEX, flags);
		690	ip = rfftwnd_create_plan(rank, n, FFTW_COMPLEX_TO_REAL, flags);
		691	}
		692
		693	CHECK(p != NULL && ip != NULL, "can't create plan");
		694
		695	for (i = 0; i < nc * nhc * 2 * MAX_STRIDE; ++i)
		696	out3[i] = 0;
		697
		698	for (istride = 1; istride <= MAX_STRIDE; ++istride) {
		699	/* generate random inputs */
		700	for (i = 0; i < nc; ++i)
		701	for (j = 0; j < nr; ++j) {
		702	c_re(in2[i * nr + j]) = DRAND();
		703	c_im(in2[i * nr + j]) = 0.0;
		704	for (k = 0; k < istride; ++k)
		705	in1[(i * nhc * 2 + j) * istride + k]
		706	= c_re(in2[i * nr + j]);
		707	}
		708
		709	fftwnd(validated_plan, 1, in2, 1, 1, out2, 1, 1);
		710
		711	howmany = ostride = istride;
		712
		713	WHEN_VERBOSE(2, printf("\n testing in-place stride %d...",
		714	istride));
		715
		716	if (howmany != 1 \|\| istride != 1 \|\| ostride != 1 \|\| coinflip())
		717	rfftwnd_real_to_complex(p, howmany, in1, istride, 1,
		718	out1, ostride, 1);
		719	else
		720	rfftwnd_one_real_to_complex(p, in1, NULL);
		721
		722	for (i = 0; i < nc; ++i)
		723	for (k = 0; k < howmany; ++k)
		724	CHECK(compute_error_complex(out1 + i * nhc * ostride + k,
		725	ostride,
		726	out2 + i * nr, 1,
		727	nhc) < TOLERANCE,
		728	"in-place (r2c): wrong answer");
		729
		730	if (howmany != 1 \|\| istride != 1 \|\| ostride != 1 \|\| coinflip())
		731	rfftwnd_complex_to_real(ip, howmany, out1, ostride, 1,
		732	out3, istride, 1);
		733	else
		734	rfftwnd_one_complex_to_real(ip, out1, NULL);
		735
		736	for (i = 0; i < nc * nhc * 2 * istride; ++i)
		737	out3[i] *= 1.0 / N;
		738
		739	for (i = 0; i < nc; ++i)
		740	for (k = 0; k < howmany; ++k)
		741	CHECK(compute_error(out3 + i * nhc * 2 * istride + k,
		742	istride,
		743	(fftw_real ) (in2 + i nr), 2,
		744	nr) < TOLERANCE,
		745	"in-place (c2r): wrong answer (check 2)");
		746	}
		747
		748	rfftwnd_destroy_plan(p);
		749	rfftwnd_destroy_plan(ip);
		750
		751	fftw_free(out2);
		752	fftw_free(in2);
		753	fftw_free(in1);
		754	}
		755
		756	void testnd_correctness(struct size sz, fftw_direction dir,
		757	int alt_api, int specific, int force_buf)
		758	{
		759	fftwnd_plan validated_plan;
		760
		761	if (dir != FFTW_FORWARD)
		762	return;
		763	if (force_buf)
		764	return;
		765
		766	validated_plan = fftwnd_create_plan(sz.rank, sz.narray,
		767	dir, measure_flag \| wisdom_flag);
		768	CHECK(validated_plan != NULL, "can't create plan");
		769
		770	testnd_out_of_place(sz.rank, sz.narray, validated_plan);
		771	testnd_in_place(sz.rank, sz.narray,
		772	validated_plan, alt_api, specific);
		773
		774	fftwnd_destroy_plan(validated_plan);
		775	}
		776
		777	/*************************************************
		778	* planner tests
		779	*************************************************/
		780
		781	void test_planner(int rank)
		782	{
		783	/*
		784	* create and destroy many plans, at random. Check the
		785	* garbage-collecting allocator of twiddle factors
		786	*/
		787	int i, dim;
		788	int r, s;
		789	fftw_plan p[PLANNER_TEST_SIZE];
		790	fftwnd_plan pnd[PLANNER_TEST_SIZE];
		791	int *narr, maxdim;
		792
		793	chk_mem_leak = 0;
		794	verbose--;
		795
		796	please_wait();
		797	if (rank < 1)
		798	rank = 1;
		799
		800	narr = (int ) fftw_malloc(rank sizeof(int));
		801
		802	maxdim = (int) pow(8192.0, 1.0/rank);
		803
		804	for (i = 0; i < PLANNER_TEST_SIZE; ++i) {
		805	p[i] = (fftw_plan) 0;
		806	pnd[i] = (fftwnd_plan) 0;
		807	}
		808
		809	for (i = 0; i < PLANNER_TEST_SIZE * PLANNER_TEST_SIZE; ++i) {
		810	r = rand();
		811	if (r < 0)
		812	r = -r;
		813	r = r % PLANNER_TEST_SIZE;
		814
		815	for (dim = 0; dim < rank; ++dim) {
		816	do {
		817	s = rand();
		818	if (s < 0)
		819	s = -s;
		820	s = s % maxdim + 1;
		821	} while (s == 0);
		822	narr[dim] = s;
		823	}
		824
		825	if (rank == 1) {
		826	if (p[r])
		827	rfftw_destroy_plan(p[r]);
		828
		829	p[r] = rfftw_create_plan(narr[0], random_dir(), measure_flag \|
		830	wisdom_flag);
		831	if (paranoid && narr[0] < 200)
		832	test_correctness(narr[0]);
		833	}
		834	if (pnd[r])
		835	rfftwnd_destroy_plan(pnd[r]);
		836
		837	pnd[r] = rfftwnd_create_plan(rank, narr,
		838	random_dir(), measure_flag \|
		839	wisdom_flag);
		840
		841	if (i % (PLANNER_TEST_SIZE * PLANNER_TEST_SIZE / 20) == 0) {
		842	WHEN_VERBOSE(0, printf("test planner: so far so good\n"));
		843	WHEN_VERBOSE(0, printf("test planner: iteration %d out of %d\n",
		844	i, PLANNER_TEST_SIZE * PLANNER_TEST_SIZE));
		845	}
		846	}
		847
		848	for (i = 0; i < PLANNER_TEST_SIZE; ++i) {
		849	if (p[i])
		850	rfftw_destroy_plan(p[i]);
		851	if (pnd[i])
		852	rfftwnd_destroy_plan(pnd[i]);
		853	}
		854
		855	fftw_free(narr);
		856	verbose++;
		857	chk_mem_leak = 1;
		858	}
		859
		860
		861	/*************************************************
		862	* Ergun's test for real->complex transforms
		863	*************************************************/
		864	static void rfill_random(fftw_real *a, int n)
		865	{
		866	int i;
		867
		868	for (i = 0; i < n; ++i) {
		869	a[i] = DRAND();
		870	}
		871	}
		872
		873	static void rarray_copy(fftw_real out, fftw_real in, int n)
		874	{
		875	int i;
		876
		877	for (i = 0; i < n; ++i)
		878	out[i] = in[i];
		879	}
		880
		881	/* C = A + B */
		882	void rarray_add(fftw_real C, fftw_real A, fftw_real *B, int n)
		883	{
		884	int i;
		885
		886	for (i = 0; i < n; ++i) {
		887	C[i] = A[i] + B[i];
		888	}
		889	}
		890
		891	/* C = A - B */
		892	void rarray_sub(fftw_real C, fftw_real A, fftw_real *B, int n)
		893	{
		894	int i;
		895
		896	for (i = 0; i < n; ++i) {
		897	C[i] = A[i] - B[i];
		898	}
		899	}
		900
		901	/* B = rotate left A */
		902	void rarray_rol(fftw_real B, fftw_real A,
		903	int n, int n_before, int n_after)
		904	{
		905	int i, ib, ia;
		906
		907	for (ib = 0; ib < n_before; ++ib) {
		908	for (i = 0; i < n - 1; ++i)
		909	for (ia = 0; ia < n_after; ++ia)
		910	B[(ib * n + i) * n_after + ia] =
		911	A[(ib * n + i + 1) * n_after + ia];
		912
		913	for (ia = 0; ia < n_after; ++ia)
		914	B[(ib * n + n - 1) * n_after + ia] = A[ib * n * n_after + ia];
		915	}
		916	}
		917
		918	/* A = alpha * B (out of place) */
		919	void rarray_scale(fftw_real A, fftw_real B, fftw_real alpha, int n)
		920	{
		921	int i;
		922
		923	for (i = 0; i < n; ++i) {
		924	A[i] = B[i] * alpha;
		925	}
		926	}
		927
		928	void rarray_compare(fftw_real A, fftw_real B, int n)
		929	{
		930	double d = compute_error(A, 1, B, 1, n);
		931	if (d > TOLERANCE) {
		932	fflush(stdout);
		933	fprintf(stderr, "Found relative error %e\n", d);
		934	fftw_die("failure in Ergun's verification procedure\n");
		935	}
		936	}
		937
		938	/*
		939	* guaranteed out-of-place transform. Does the necessary
		940	* copying if the plan is in-place.
		941	*/
		942	static void rfftw_out_of_place(fftw_plan plan, int n,
		943	fftw_real in, fftw_real out)
		944	{
		945	if (plan->flags & FFTW_IN_PLACE) {
		946	rarray_copy(out, in, n);
		947	rfftw(plan, 1, out, 1, n, (fftw_real *)0, 1, n);
		948	} else {
		949	rfftw(plan, 1, in, 1, n, out, 1, n);
		950	}
		951	}
		952
		953	/*
		954	* This is a real (as opposed to complex) variation of the FFT tester
		955	* described in
		956	*
		957	* Funda Erg�n. Testing multivariate linear functions: Overcoming the
		958	* generator bottleneck. In Proceedings of the Twenty-Seventh Annual
		959	* ACM Symposium on the Theory of Computing, pages 407-416, Las Vegas,
		960	* Nevada, 29 May--1 June 1995.
		961	*/
		962	void test_ergun(int n, fftw_direction dir, fftw_plan plan)
		963	{
		964	fftw_real inA, inB, inC, outA, outB, outC;
		965	fftw_real inA1, inB1;
		966	fftw_real *tmp;
		967	int i;
		968	int rounds = 20;
		969	FFTW_TRIG_REAL twopin = FFTW_K2PI / (FFTW_TRIG_REAL) n;
		970
		971	inA = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		972	inB = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		973	inA1 = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		974	inB1 = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		975	inC = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		976	outA = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		977	outB = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		978	outC = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		979	tmp = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		980
		981	WHEN_VERBOSE(2,
		982	printf("Validating plan, n = %d, dir = %s\n", n,
		983	dir == FFTW_REAL_TO_COMPLEX ?
		984	"REAL_TO_COMPLEX" : "COMPLEX_TO_REAL"));
		985
		986	/* test 1: check linearity */
		987	for (i = 0; i < rounds; ++i) {
		988	fftw_real alpha, beta;
		989	alpha = DRAND();
		990	beta = DRAND();
		991	rfill_random(inA, n);
		992	rfill_random(inB, n);
		993	rarray_scale(inA1, inA, alpha, n);
		994	rarray_scale(inB1, inB, beta, n);
		995	rarray_add(inC, inA1, inB1, n);
		996	rfftw_out_of_place(plan, n, inA, outA);
		997	rfftw_out_of_place(plan, n, inB, outB);
		998	rarray_scale(outA, outA, alpha, n);
		999	rarray_scale(outB, outB, beta, n);
		1000	rarray_add(tmp, outA, outB, n);
		1001	rfftw_out_of_place(plan, n, inC, outC);
		1002	rarray_compare(outC, tmp, n);
		1003	}
		1004
		1005	/* test 2: check that the unit impulse is transformed properly */
		1006	for (i = 0; i < n; ++i) {
		1007	/* impulse */
		1008	inA[i] = 0.0;
		1009
		1010	/* transform of the impulse */
		1011	if (2 * i <= n)
		1012	outA[i] = 1.0;
		1013	else
		1014	outA[i] = 0.0;
		1015	}
		1016	inA[0] = 1.0;
		1017
		1018	if (dir == FFTW_REAL_TO_COMPLEX) {
		1019	for (i = 0; i < rounds; ++i) {
		1020	rfill_random(inB, n);
		1021	rarray_sub(inC, inA, inB, n);
		1022	rfftw_out_of_place(plan, n, inB, outB);
		1023	rfftw_out_of_place(plan, n, inC, outC);
		1024	rarray_add(tmp, outB, outC, n);
		1025	rarray_compare(tmp, outA, n);
		1026	}
		1027	} else {
		1028	for (i = 0; i < rounds; ++i) {
		1029	rfill_random(outB, n);
		1030	rarray_sub(outC, outA, outB, n);
		1031	rfftw_out_of_place(plan, n, outB, inB);
		1032	rfftw_out_of_place(plan, n, outC, inC);
		1033	rarray_add(tmp, inB, inC, n);
		1034	rarray_scale(tmp, tmp, 1.0 / ((double) n), n);
		1035	rarray_compare(tmp, inA, n);
		1036	}
		1037	}
		1038
		1039	/* test 3: check the time-shift property */
		1040	/* the paper performs more tests, but this code should be fine too */
		1041	if (dir == FFTW_REAL_TO_COMPLEX) {
		1042	for (i = 0; i < rounds; ++i) {
		1043	int j;
		1044
		1045	rfill_random(inA, n);
		1046	rarray_rol(inB, inA, n, 1, 1);
		1047	rfftw_out_of_place(plan, n, inA, outA);
		1048	rfftw_out_of_place(plan, n, inB, outB);
		1049	tmp[0] = outB[0];
		1050	for (j = 1; 2 * j < n; ++j) {
		1051	FFTW_TRIG_REAL s = dir * FFTW_TRIG_SIN(j * twopin);
		1052	FFTW_TRIG_REAL c = FFTW_TRIG_COS(j * twopin);
		1053	tmp[j] = outB[j] * c - outB[n - j] * s;
		1054	tmp[n - j] = outB[j] * s + outB[n - j] * c;
		1055	}
		1056	if (2 * j == n)
		1057	tmp[j] = -outB[j];
		1058
		1059	rarray_compare(tmp, outA, n);
		1060	}
		1061	} else {
		1062	for (i = 0; i < rounds; ++i) {
		1063	int j;
		1064
		1065	rfill_random(inA, n);
		1066	inB[0] = inA[0];
		1067	for (j = 1; 2 * j < n; ++j) {
		1068	FFTW_TRIG_REAL s = dir * FFTW_TRIG_SIN(j * twopin);
		1069	FFTW_TRIG_REAL c = FFTW_TRIG_COS(j * twopin);
		1070	inB[j] = inA[j] * c - inA[n - j] * s;
		1071	inB[n - j] = inA[j] * s + inA[n - j] * c;
		1072	}
		1073	if (2 * j == n)
		1074	inB[j] = -inA[j];
		1075
		1076	rfftw_out_of_place(plan, n, inA, outA);
		1077	rfftw_out_of_place(plan, n, inB, outB);
		1078	rarray_rol(tmp, outA, n, 1, 1);
		1079	rarray_compare(tmp, outB, n);
		1080	}
		1081	}
		1082
		1083	WHEN_VERBOSE(2, printf("Validation done\n"));
		1084
		1085	fftw_free(tmp);
		1086	fftw_free(outC);
		1087	fftw_free(outB);
		1088	fftw_free(outA);
		1089	fftw_free(inC);
		1090	fftw_free(inB1);
		1091	fftw_free(inA1);
		1092	fftw_free(inB);
		1093	fftw_free(inA);
		1094	}
		1095
		1096	static void rfftw_plan_hook_function(fftw_plan p)
		1097	{
		1098	WHEN_VERBOSE(3, printf("Validating tentative plan\n"));
		1099	WHEN_VERBOSE(3, fftw_print_plan(p));
		1100	test_ergun(p->n, p->dir, p);
		1101	}
		1102
		1103	/*************************************************
		1104	* test initialization
		1105	*************************************************/
		1106
		1107	void test_init(int argc, char *argv)
		1108	{
		1109	}
		1110
		1111	void test_finish(void)
		1112	{
		1113	}
		1114
		1115	extern void (*rfftw_plan_hook)(fftw_plan plan);
		1116	void enter_paranoid_mode(void)
		1117	{
		1118	rfftw_plan_hook = rfftw_plan_hook_function;
		1119	}
		1120
		1121	int get_option(int argc, char *argv, char argval, int argval_maxlen)
		1122	{
		1123	return default_get_option(argc,argv,argval,argval_maxlen);
		1124	}

Subversion Repositories shark

(root)/shark/branches/pj/ports/fftw/tests/rfftw_te.c - Rev 2