WebSVN - shark - Blame - Rev 3 - /shark/trunk/ports/fftw/rfftw/rfftwnd.c

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	/* $Id: rfftwnd.c,v 1.1.1.1 2002-03-29 14:12:59 pj Exp $ */
		21
		22	#include <ports/fftw-int.h>
		23	#include <ports/rfftw.h>
		24
		25	/******************** prototypes for rexec2 routines ********************/
		26
		27	extern void rfftw_real2c_aux(fftw_plan plan, int howmany,
		28	fftw_real *in, int istride, int idist,
		29	fftw_complex *out, int ostride, int odist,
		30	fftw_real *work);
		31	extern void rfftw_c2real_aux(fftw_plan plan, int howmany,
		32	fftw_complex *in, int istride, int idist,
		33	fftw_real *out, int ostride, int odist,
		34	fftw_real *work);
		35	extern void rfftw_real2c_overlap_aux(fftw_plan plan, int howmany,
		36	fftw_real *in, int istride, int idist,
		37	fftw_complex *out, int ostride, int odist,
		38	fftw_real *work);
		39	extern void rfftw_c2real_overlap_aux(fftw_plan plan, int howmany,
		40	fftw_complex *in, int istride, int idist,
		41	fftw_real *out, int ostride, int odist,
		42	fftw_real *work);
		43
		44	/******************** Initializing the RFFTWND Plan *********************/
		45
		46	/*
		47	* Create an fftwnd_plan specialized for specific arrays. (These
		48	* arrays are ignored, however, if they are NULL or if the flags
		49	* do not include FFTW_MEASURE.) The main advantage of being
		50	* provided arrays like this is that we can do runtime timing
		51	* measurements of our options, without worrying about allocating
		52	* excessive scratch space.
		53	*/
		54	fftwnd_plan rfftwnd_create_plan_specific(int rank, const int *n,
		55	fftw_direction dir, int flags,
		56	fftw_real *in, int istride,
		57	fftw_real *out, int ostride)
		58	{
		59	fftwnd_plan p;
		60	int i;
		61	int rflags = flags & ~FFTW_IN_PLACE;
		62	/* note that we always do rfftw transforms out-of-place in rexec2.c */
		63
		64	if (flags & FFTW_IN_PLACE) {
		65	out = NULL;
		66	ostride = istride;
		67	}
		68	istride = ostride = 1; /*
		69	* strides don't work yet, since it is not
		70	* clear whether they apply to real
		71	* or complex data
		72	*/
		73
		74	if (!(p = fftwnd_create_plan_aux(rank, n, dir, flags)))
		75	return 0;
		76
		77	for (i = 0; i < rank - 1; ++i)
		78	p->n_after[i] = (n[rank - 1]/2 + 1) * (p->n_after[i] / n[rank - 1]);
		79	if (rank > 0)
		80	p->n[rank - 1] = n[rank - 1] / 2 + 1;
		81
		82	p->plans = fftwnd_new_plan_array(rank);
		83	if (rank > 0 && !p->plans) {
		84	rfftwnd_destroy_plan(p);
		85	return 0;
		86	}
		87	if (rank > 0) {
		88	p->plans[rank - 1] = rfftw_create_plan(n[rank - 1], dir, rflags);
		89	if (!p->plans[rank - 1]) {
		90	rfftwnd_destroy_plan(p);
		91	return 0;
		92	}
		93	}
		94	if (rank > 1) {
		95	if (!(flags & FFTW_MEASURE) \|\| in == 0
		96	\|\| (!p->is_in_place && out == 0)) {
		97	if (!fftwnd_create_plans_generic(p->plans, rank - 1, n,
		98	dir, flags \| FFTW_IN_PLACE)) {
		99	rfftwnd_destroy_plan(p);
		100	return 0;
		101	}
		102	} else if (dir == FFTW_COMPLEX_TO_REAL \|\| (flags & FFTW_IN_PLACE)) {
		103	if (!fftwnd_create_plans_specific(p->plans, rank - 1, n,
		104	p->n_after,
		105	dir, flags \| FFTW_IN_PLACE,
		106	(fftw_complex *) in,
		107	istride,
		108	0, 0)) {
		109	rfftwnd_destroy_plan(p);
		110	return 0;
		111	}
		112	} else {
		113	if (!fftwnd_create_plans_specific(p->plans, rank - 1, n,
		114	p->n_after,
		115	dir, flags \| FFTW_IN_PLACE,
		116	(fftw_complex *) out,
		117	ostride,
		118	0, 0)) {
		119	rfftwnd_destroy_plan(p);
		120	return 0;
		121	}
		122	}
		123	}
		124	p->nbuffers = 0;
		125	p->nwork = fftwnd_work_size(rank, p->n, flags \| FFTW_IN_PLACE,
		126	p->nbuffers + 1);
		127	if (p->nwork && !(flags & FFTW_THREADSAFE)) {
		128	p->work = (fftw_complex *) fftw_malloc(p->nwork
		129	* sizeof(fftw_complex));
		130	if (!p->work) {
		131	rfftwnd_destroy_plan(p);
		132	return 0;
		133	}
		134	}
		135	return p;
		136	}
		137
		138	fftwnd_plan rfftw2d_create_plan_specific(int nx, int ny,
		139	fftw_direction dir, int flags,
		140	fftw_real *in, int istride,
		141	fftw_real *out, int ostride)
		142	{
		143	int n[2];
		144
		145	n[0] = nx;
		146	n[1] = ny;
		147
		148	return rfftwnd_create_plan_specific(2, n, dir, flags,
		149	in, istride, out, ostride);
		150	}
		151
		152	fftwnd_plan rfftw3d_create_plan_specific(int nx, int ny, int nz,
		153	fftw_direction dir, int flags,
		154	fftw_real *in, int istride,
		155	fftw_real *out, int ostride)
		156	{
		157	int n[3];
		158
		159	n[0] = nx;
		160	n[1] = ny;
		161	n[2] = nz;
		162
		163	return rfftwnd_create_plan_specific(3, n, dir, flags,
		164	in, istride, out, ostride);
		165	}
		166
		167	/* Create a generic fftwnd plan: */
		168
		169	fftwnd_plan rfftwnd_create_plan(int rank, const int *n,
		170	fftw_direction dir, int flags)
		171	{
		172	return rfftwnd_create_plan_specific(rank, n, dir, flags, 0, 1, 0, 1);
		173	}
		174
		175	fftwnd_plan rfftw2d_create_plan(int nx, int ny,
		176	fftw_direction dir, int flags)
		177	{
		178	return rfftw2d_create_plan_specific(nx, ny, dir, flags, 0, 1, 0, 1);
		179	}
		180
		181	fftwnd_plan rfftw3d_create_plan(int nx, int ny, int nz,
		182	fftw_direction dir, int flags)
		183	{
		184	return rfftw3d_create_plan_specific(nx, ny, nz, dir, flags, 0, 1, 0, 1);
		185	}
		186
		187	/********************** Freeing the RFFTWND Plan **********************/
		188
		189	void rfftwnd_destroy_plan(fftwnd_plan plan)
		190	{
		191	fftwnd_destroy_plan(plan);
		192	}
		193
		194	/********************** Printing the RFFTWND Plan **********************/
		195	/*
		196	void rfftwnd_fprint_plan(FILE *f, fftwnd_plan plan)
		197	{
		198	fftwnd_fprint_plan(f, plan);
		199	}
		200
		201	void rfftwnd_print_plan(fftwnd_plan plan)
		202	{
		203	rfftwnd_fprint_plan(stdout, plan);
		204	}
		205	*/
		206	/********* Computing the N-Dimensional FFT: Auxiliary Routines **********/
		207
		208	void rfftwnd_real2c_aux(fftwnd_plan p, int cur_dim,
		209	fftw_real *in, int istride,
		210	fftw_complex *out, int ostride,
		211	fftw_real *work)
		212	{
		213	int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
		214
		215	if (cur_dim == p->rank - 2) {
		216	/* just do the last dimension directly: */
		217	if (p->is_in_place)
		218	rfftw_real2c_aux(p->plans[p->rank - 1], n,
		219	in, istride, (n_after * istride) * 2,
		220	out, istride, n_after * istride,
		221	work);
		222	else
		223	rfftw_real2c_aux(p->plans[p->rank - 1], n,
		224	in, istride, p->plans[p->rank - 1]->n * istride,
		225	out, ostride, n_after * ostride,
		226	work);
		227	} else { /* we have at least two dimensions to go */
		228	int nr = p->plans[p->rank - 1]->n;
		229	int n_after_r = p->is_in_place ? n_after * 2
		230	: nr * (n_after / (nr/2 + 1));
		231	int i;
		232
		233	/*
		234	* process the subsequent dimensions recursively, in hyperslabs,
		235	* to get maximum locality:
		236	*/
		237	for (i = 0; i < n; ++i)
		238	rfftwnd_real2c_aux(p, cur_dim + 1,
		239	in + i * n_after_r * istride, istride,
		240	out + i * n_after * ostride, ostride, work);
		241	}
		242
		243	/* do the current dimension (in-place): */
		244	fftw(p->plans[cur_dim], n_after,
		245	out, n_after * ostride, ostride,
		246	(fftw_complex *) work, 1, 0);
		247	/* I hate this cast */
		248	}
		249
		250	void rfftwnd_c2real_aux(fftwnd_plan p, int cur_dim,
		251	fftw_complex *in, int istride,
		252	fftw_real *out, int ostride,
		253	fftw_real *work)
		254	{
		255	int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
		256
		257	/* do the current dimension (in-place): */
		258	fftw(p->plans[cur_dim], n_after,
		259	in, n_after * istride, istride,
		260	(fftw_complex *) work, 1, 0);
		261
		262	if (cur_dim == p->rank - 2) {
		263	/* just do the last dimension directly: */
		264	if (p->is_in_place)
		265	rfftw_c2real_aux(p->plans[p->rank - 1], n,
		266	in, istride, n_after * istride,
		267	out, istride, (n_after * istride) * 2,
		268	work);
		269	else
		270	rfftw_c2real_aux(p->plans[p->rank - 1], n,
		271	in, istride, n_after * istride,
		272	out, ostride, p->plans[p->rank - 1]->n * ostride,
		273	work);
		274	} else { /* we have at least two dimensions to go */
		275	int nr = p->plans[p->rank - 1]->n;
		276	int n_after_r = p->is_in_place ? n_after * 2 :
		277	nr * (n_after / (nr/2 + 1));
		278	int i;
		279
		280	/*
		281	* process the subsequent dimensions recursively, in hyperslabs,
		282	* to get maximum locality:
		283	*/
		284	for (i = 0; i < n; ++i)
		285	rfftwnd_c2real_aux(p, cur_dim + 1,
		286	in + i * n_after * istride, istride,
		287	out + i * n_after_r * ostride, ostride, work);
		288	}
		289	}
		290
		291	/*
		292	* alternate version of rfftwnd_aux -- this version pushes the howmany
		293	* loop down to the leaves of the computation, for greater locality
		294	* in cases where dist < stride. It is also required for correctness
		295	* if in==out, and we must call a special version of the executor.
		296	* Note that work must point to 'howmany' copies of its data
		297	* if in == out.
		298	*/
		299
		300	void rfftwnd_real2c_aux_howmany(fftwnd_plan p, int cur_dim,
		301	int howmany,
		302	fftw_real *in, int istride, int idist,
		303	fftw_complex *out, int ostride, int odist,
		304	fftw_complex *work)
		305	{
		306	int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
		307	int k;
		308
		309	if (cur_dim == p->rank - 2) {
		310	/* just do the last dimension directly: */
		311	if (p->is_in_place)
		312	for (k = 0; k < n; ++k)
		313	rfftw_real2c_overlap_aux(p->plans[p->rank - 1], howmany,
		314	in + (k * n_after * istride) * 2,
		315	istride, idist,
		316	out + (k * n_after * ostride),
		317	ostride, odist,
		318	(fftw_real *) work);
		319	else {
		320	int nlast = p->plans[p->rank - 1]->n;
		321	for (k = 0; k < n; ++k)
		322	rfftw_real2c_aux(p->plans[p->rank - 1], howmany,
		323	in + k * nlast * istride,
		324	istride, idist,
		325	out + k * n_after * ostride,
		326	ostride, odist,
		327	(fftw_real *) work);
		328	}
		329	} else { /* we have at least two dimensions to go */
		330	int nr = p->plans[p->rank - 1]->n;
		331	int n_after_r = p->is_in_place ? n_after * 2 :
		332	nr * (n_after / (nr/2 + 1));
		333	int i;
		334
		335	/*
		336	* process the subsequent dimensions recursively, in hyperslabs,
		337	* to get maximum locality:
		338	*/
		339	for (i = 0; i < n; ++i)
		340	rfftwnd_real2c_aux_howmany(p, cur_dim + 1, howmany,
		341	in + i * n_after_r * istride, istride, idist,
		342	out + i * n_after * ostride, ostride, odist,
		343	work);
		344	}
		345
		346	/* do the current dimension (in-place): */
		347	for (k = 0; k < n_after; ++k)
		348	fftw(p->plans[cur_dim], howmany,
		349	out + k * ostride, n_after * ostride, odist,
		350	work, 1, 0);
		351	}
		352
		353	void rfftwnd_c2real_aux_howmany(fftwnd_plan p, int cur_dim,
		354	int howmany,
		355	fftw_complex *in, int istride, int idist,
		356	fftw_real *out, int ostride, int odist,
		357	fftw_complex *work)
		358	{
		359	int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
		360	int k;
		361
		362	/* do the current dimension (in-place): */
		363	for (k = 0; k < n_after; ++k)
		364	fftw(p->plans[cur_dim], howmany,
		365	in + k * istride, n_after * istride, idist,
		366	work, 1, 0);
		367
		368	if (cur_dim == p->rank - 2) {
		369	/* just do the last dimension directly: */
		370	if (p->is_in_place)
		371	for (k = 0; k < n; ++k)
		372	rfftw_c2real_overlap_aux(p->plans[p->rank - 1], howmany,
		373	in + (k * n_after * istride),
		374	istride, idist,
		375	out + (k * n_after * ostride) * 2,
		376	ostride, odist,
		377	(fftw_real *) work);
		378	else {
		379	int nlast = p->plans[p->rank - 1]->n;
		380	for (k = 0; k < n; ++k)
		381	rfftw_c2real_aux(p->plans[p->rank - 1], howmany,
		382	in + k * n_after * istride,
		383	istride, idist,
		384	out + k * nlast * ostride,
		385	ostride, odist,
		386	(fftw_real *) work);
		387	}
		388	} else { /* we have at least two dimensions to go */
		389	int nr = p->plans[p->rank - 1]->n;
		390	int n_after_r = p->is_in_place ? n_after * 2
		391	: nr * (n_after / (nr/2 + 1));
		392	int i;
		393
		394	/*
		395	* process the subsequent dimensions recursively, in hyperslabs,
		396	* to get maximum locality:
		397	*/
		398	for (i = 0; i < n; ++i)
		399	rfftwnd_c2real_aux_howmany(p, cur_dim + 1, howmany,
		400	in + i * n_after * istride, istride, idist,
		401	out + i * n_after_r * ostride, ostride, odist,
		402	work);
		403	}
		404	}
		405
		406	/******** Computing the N-Dimensional FFT: User-Visible Routines ********/
		407
		408	void rfftwnd_real_to_complex(fftwnd_plan p, int howmany,
		409	fftw_real *in, int istride, int idist,
		410	fftw_complex *out, int ostride, int odist)
		411	{
		412	fftw_complex *work = p->work;
		413	int rank = p->rank;
		414	int free_work = 0;
		415
		416	if (p->dir != FFTW_REAL_TO_COMPLEX)
		417	fftw_die("rfftwnd_real_to_complex with complex-to-real plan");
		418
		419	#ifdef FFTW_DEBUG
		420	if (p->rank > 0 && (p->plans[0]->flags & FFTW_THREADSAFE)
		421	&& p->nwork && p->work)
		422	fftw_die("bug with FFTW_THREADSAFE flag");
		423	#endif
		424
		425	if (p->is_in_place) {
		426	ostride = istride;
		427	odist = (idist == 1) ? 1 : (idist / 2); /* ugh */
		428	out = (fftw_complex *) in;
		429	if (howmany > 1 && istride > idist && rank > 0) {
		430	int new_nwork;
		431
		432	new_nwork = p->n[rank - 1] * howmany;
		433	if (new_nwork > p->nwork) {
		434	work = (fftw_complex *)
		435	fftw_malloc(sizeof(fftw_complex) * new_nwork);
		436	if (!work)
		437	fftw_die("error allocating work array");
		438	free_work = 1;
		439	}
		440	}
		441	}
		442	if (p->nwork && !work) {
		443	work = (fftw_complex ) fftw_malloc(sizeof(fftw_complex) p->nwork);
		444	free_work = 1;
		445	}
		446	switch (rank) {
		447	case 0:
		448	break;
		449	case 1:
		450	if (p->is_in_place && howmany > 1 && istride > idist)
		451	rfftw_real2c_overlap_aux(p->plans[0], howmany,
		452	in, istride, idist,
		453	out, ostride, odist,
		454	(fftw_real *) work);
		455	else
		456	rfftw_real2c_aux(p->plans[0], howmany,
		457	in, istride, idist,
		458	out, ostride, odist,
		459	(fftw_real *) work);
		460	break;
		461	default: /* rank >= 2 */
		462	{
		463	if (howmany > 1 && ostride > odist)
		464	rfftwnd_real2c_aux_howmany(p, 0, howmany,
		465	in, istride, idist,
		466	out, ostride, odist,
		467	work);
		468	else {
		469	int i;
		470
		471	for (i = 0; i < howmany; ++i)
		472	rfftwnd_real2c_aux(p, 0,
		473	in + i * idist, istride,
		474	out + i * odist, ostride,
		475	(fftw_real *) work);
		476	}
		477	}
		478	}
		479
		480	if (free_work)
		481	fftw_free(work);
		482	}
		483
		484	void rfftwnd_complex_to_real(fftwnd_plan p, int howmany,
		485	fftw_complex *in, int istride, int idist,
		486	fftw_real *out, int ostride, int odist)
		487	{
		488	fftw_complex *work = p->work;
		489	int rank = p->rank;
		490	int free_work = 0;
		491
		492	if (p->dir != FFTW_COMPLEX_TO_REAL)
		493	fftw_die("rfftwnd_complex_to_real with real-to-complex plan");
		494
		495	#ifdef FFTW_DEBUG
		496	if (p->rank > 0 && (p->plans[0]->flags & FFTW_THREADSAFE)
		497	&& p->nwork && p->work)
		498	fftw_die("bug with FFTW_THREADSAFE flag");
		499	#endif
		500
		501	if (p->is_in_place) {
		502	ostride = istride;
		503	odist = idist;
		504	odist = (idist == 1) ? 1 : (idist * 2); /* ugh */
		505	out = (fftw_real *) in;
		506	if (howmany > 1 && istride > idist && rank > 0) {
		507	int new_nwork = p->n[rank - 1] * howmany;
		508	if (new_nwork > p->nwork) {
		509	work = (fftw_complex *)
		510	fftw_malloc(sizeof(fftw_complex) * new_nwork);
		511	if (!work)
		512	fftw_die("error allocating work array");
		513	free_work = 1;
		514	}
		515	}
		516	}
		517	if (p->nwork && !work) {
		518	work = (fftw_complex ) fftw_malloc(sizeof(fftw_complex) p->nwork);
		519	free_work = 1;
		520	}
		521	switch (rank) {
		522	case 0:
		523	break;
		524	case 1:
		525	if (p->is_in_place && howmany > 1 && istride > idist)
		526	rfftw_c2real_overlap_aux(p->plans[0], howmany,
		527	in, istride, idist,
		528	out, ostride, odist,
		529	(fftw_real *) work);
		530	else
		531	rfftw_c2real_aux(p->plans[0], howmany,
		532	in, istride, idist,
		533	out, ostride, odist,
		534	(fftw_real *) work);
		535	break;
		536	default: /* rank >= 2 */
		537	{
		538	if (howmany > 1 && ostride > odist)
		539	rfftwnd_c2real_aux_howmany(p, 0, howmany,
		540	in, istride, idist,
		541	out, ostride, odist,
		542	work);
		543	else {
		544	int i;
		545
		546	for (i = 0; i < howmany; ++i)
		547	rfftwnd_c2real_aux(p, 0,
		548	in + i * idist, istride,
		549	out + i * odist, ostride,
		550	(fftw_real *) work);
		551	}
		552	}
		553	}
		554
		555	if (free_work)
		556	fftw_free(work);
		557	}
		558
		559	void rfftwnd_one_real_to_complex(fftwnd_plan p,
		560	fftw_real in, fftw_complex out)
		561	{
		562	rfftwnd_real_to_complex(p, 1, in, 1, 1, out, 1, 1);
		563	}
		564
		565	void rfftwnd_one_complex_to_real(fftwnd_plan p,
		566	fftw_complex in, fftw_real out)
		567	{
		568	rfftwnd_complex_to_real(p, 1, in, 1, 1, out, 1, 1);
		569	}

Subversion Repositories shark

(root)/shark/trunk/ports/fftw/rfftw/rfftwnd.c - Rev 3