WebSVN - shark - Blame - Rev 2 - /shark/branches/pj/ports/fftw/rfftw/rexec.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	/*
		21	* rexec.c -- execute the fft
		22	*/
		23
		24	/* $Id: rexec.c,v 1.1.1.1 2002-03-29 14:12:59 pj Exp $ */
		25	//#include <ports/stdio.h>
		26	#include <stdlib.h>
		27
		28	#include <ports/fftw-int.h>
		29	#include <ports/rfftw.h>
		30
		31	void rfftw_strided_copy(int n, fftw_real *in, int ostride,
		32	fftw_real *out)
		33	{
		34	int i;
		35	fftw_real r0, r1, r2, r3;
		36
		37	i = 0;
		38	for (; i < (n & 3); ++i) {
		39	out[i * ostride] = in[i];
		40	}
		41	for (; i < n; i += 4) {
		42	r0 = in[i];
		43	r1 = in[i + 1];
		44	r2 = in[i + 2];
		45	r3 = in[i + 3];
		46	out[i * ostride] = r0;
		47	out[(i + 1) * ostride] = r1;
		48	out[(i + 2) * ostride] = r2;
		49	out[(i + 3) * ostride] = r3;
		50	}
		51	}
		52
		53	void rfftw_executor_simple(int n, fftw_real *in,
		54	fftw_real *out,
		55	fftw_plan_node *p,
		56	int istride,
		57	int ostride)
		58	{
		59	switch (p->type) {
		60	case FFTW_REAL2HC:
		61	HACK_ALIGN_STACK_ODD();
		62	(p->nodeu.real2hc.codelet) (in, out, out + n * ostride,
		63	istride, ostride, -ostride);
		64	break;
		65
		66	case FFTW_HC2REAL:
		67	HACK_ALIGN_STACK_ODD();
		68	(p->nodeu.hc2real.codelet) (in, in + n * istride, out,
		69	istride, -istride, ostride);
		70	break;
		71
		72	case FFTW_HC2HC:
		73	{
		74	int r = p->nodeu.hc2hc.size;
		75	int m = n / r;
		76	int i;
		77	/*
		78	* please do resist the temptation of initializing
		79	* these variables here. Doing so forces the
		80	* compiler to keep a live variable across the
		81	* recursive call.
		82	*/
		83	fftw_hc2hc_codelet *codelet;
		84	fftw_complex *W;
		85
		86	switch (p->nodeu.hc2hc.dir) {
		87	case FFTW_REAL_TO_COMPLEX:
		88	for (i = 0; i < r; ++i)
		89	rfftw_executor_simple(m, in + i * istride,
		90	out + i * (m*ostride),
		91	p->nodeu.hc2hc.recurse,
		92	istride * r, ostride);
		93
		94	W = p->nodeu.hc2hc.tw->twarray;
		95	codelet = p->nodeu.hc2hc.codelet;
		96	HACK_ALIGN_STACK_EVEN();
		97	codelet(out, W, m * ostride, m, ostride);
		98	break;
		99	case FFTW_COMPLEX_TO_REAL:
		100	W = p->nodeu.hc2hc.tw->twarray;
		101	codelet = p->nodeu.hc2hc.codelet;
		102	HACK_ALIGN_STACK_EVEN();
		103	codelet(in, W, m * istride, m, istride);
		104
		105	for (i = 0; i < r; ++i)
		106	rfftw_executor_simple(m, in + i * (m*istride),
		107	out + i * ostride,
		108	p->nodeu.hc2hc.recurse,
		109	istride, ostride * r);
		110	break;
		111	default:
		112	goto bug;
		113	}
		114
		115	break;
		116	}
		117
		118	case FFTW_RGENERIC:
		119	{
		120	int r = p->nodeu.rgeneric.size;
		121	int m = n / r;
		122	int i;
		123	fftw_rgeneric_codelet *codelet = p->nodeu.rgeneric.codelet;
		124	fftw_complex *W = p->nodeu.rgeneric.tw->twarray;
		125
		126	switch (p->nodeu.rgeneric.dir) {
		127	case FFTW_REAL_TO_COMPLEX:
		128	for (i = 0; i < r; ++i)
		129	rfftw_executor_simple(m, in + i * istride,
		130	out + i * (m * ostride),
		131	p->nodeu.rgeneric.recurse,
		132	istride * r, ostride);
		133
		134	codelet(out, W, m, r, n, ostride);
		135	break;
		136	case FFTW_COMPLEX_TO_REAL:
		137	codelet(in, W, m, r, n, istride);
		138
		139	for (i = 0; i < r; ++i)
		140	rfftw_executor_simple(m, in + i * m * istride,
		141	out + i * ostride,
		142	p->nodeu.rgeneric.recurse,
		143	istride, ostride * r);
		144	break;
		145	default:
		146	goto bug;
		147	}
		148
		149	break;
		150	}
		151
		152	default:
		153	bug:
		154	fftw_die("BUG in rexecutor: invalid plan\n");
		155	break;
		156	}
		157	}
		158
		159	static void rexecutor_simple_inplace(int n, fftw_real *in,
		160	fftw_real *out,
		161	fftw_plan_node *p,
		162	int istride)
		163	{
		164	switch (p->type) {
		165	case FFTW_REAL2HC:
		166	HACK_ALIGN_STACK_ODD();
		167	(p->nodeu.real2hc.codelet) (in, in, in + n * istride,
		168	istride, istride, -istride);
		169	break;
		170
		171	case FFTW_HC2REAL:
		172	HACK_ALIGN_STACK_ODD();
		173	(p->nodeu.hc2real.codelet) (in, in + n * istride, in,
		174	istride, -istride, istride);
		175	break;
		176
		177	default:
		178	{
		179	fftw_real *tmp;
		180
		181	if (out)
		182	tmp = out;
		183	else
		184	tmp = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		185
		186	rfftw_executor_simple(n, in, tmp, p, istride, 1);
		187	rfftw_strided_copy(n, tmp, istride, in);
		188
		189	if (!out)
		190	fftw_free(tmp);
		191	}
		192	}
		193	}
		194
		195	static void rexecutor_many(int n, fftw_real *in,
		196	fftw_real *out,
		197	fftw_plan_node *p,
		198	int istride,
		199	int ostride,
		200	int howmany, int idist, int odist)
		201	{
		202	switch (p->type) {
		203	case FFTW_REAL2HC:
		204	{
		205	fftw_real2hc_codelet *codelet = p->nodeu.real2hc.codelet;
		206	int s;
		207
		208	HACK_ALIGN_STACK_ODD();
		209	for (s = 0; s < howmany; ++s)
		210	codelet(in + s * idist, out + s * odist,
		211	out + n * ostride + s * odist,
		212	istride, ostride, -ostride);
		213	break;
		214	}
		215
		216	case FFTW_HC2REAL:
		217	{
		218	fftw_hc2real_codelet *codelet = p->nodeu.hc2real.codelet;
		219	int s;
		220
		221	HACK_ALIGN_STACK_ODD();
		222	for (s = 0; s < howmany; ++s)
		223	codelet(in + s * idist, in + n * istride + s * idist,
		224	out + s * odist,
		225	istride, -istride, ostride);
		226	break;
		227	}
		228
		229	default:
		230	{
		231	int s;
		232	for (s = 0; s < howmany; ++s) {
		233	rfftw_executor_simple(n, in + s * idist,
		234	out + s * odist,
		235	p, istride, ostride);
		236	}
		237	}
		238	}
		239	}
		240
		241	static void rexecutor_many_inplace(int n, fftw_real *in,
		242	fftw_real *out,
		243	fftw_plan_node *p,
		244	int istride,
		245	int howmany, int idist)
		246	{
		247	switch (p->type) {
		248	case FFTW_REAL2HC:
		249	{
		250	fftw_real2hc_codelet *codelet = p->nodeu.real2hc.codelet;
		251	int s;
		252
		253	HACK_ALIGN_STACK_ODD();
		254	for (s = 0; s < howmany; ++s)
		255	codelet(in + s * idist, in + s * idist,
		256	in + n * istride + s * idist,
		257	istride, istride, -istride);
		258
		259	break;
		260	}
		261
		262	case FFTW_HC2REAL:
		263	{
		264	fftw_hc2real_codelet *codelet = p->nodeu.hc2real.codelet;
		265	int s;
		266
		267	HACK_ALIGN_STACK_ODD();
		268	for (s = 0; s < howmany; ++s)
		269	codelet(in + s * idist, in + n * istride + s * idist,
		270	in + s * idist,
		271	istride, -istride, istride);
		272
		273	break;
		274	}
		275
		276	default:
		277	{
		278	int s;
		279	fftw_real *tmp;
		280	if (out)
		281	tmp = out;
		282	else
		283	tmp = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		284
		285	for (s = 0; s < howmany; ++s) {
		286	rfftw_executor_simple(n,
		287	in + s * idist,
		288	tmp,
		289	p, istride, 1);
		290	rfftw_strided_copy(n, tmp, istride, in + s * idist);
		291	}
		292
		293	if (!out)
		294	fftw_free(tmp);
		295	}
		296	}
		297	}
		298
		299	/* user interface */
		300	void rfftw(fftw_plan plan, int howmany, fftw_real *in, int istride,
		301	int idist, fftw_real *out, int ostride, int odist)
		302	{
		303	int n = plan->n;
		304
		305	if (plan->flags & FFTW_IN_PLACE) {
		306	if (howmany == 1) {
		307	rexecutor_simple_inplace(n, in, out, plan->root, istride);
		308	} else {
		309	rexecutor_many_inplace(n, in, out, plan->root, istride, howmany,
		310	idist);
		311	}
		312	} else {
		313	if (howmany == 1) {
		314	rfftw_executor_simple(n, in, out, plan->root, istride, ostride);
		315	} else {
		316	rexecutor_many(n, in, out, plan->root, istride, ostride,
		317	howmany, idist, odist);
		318	}
		319	}
		320	}
		321
		322	void rfftw_one(fftw_plan plan, fftw_real in, fftw_real out)
		323	{
		324	int n = plan->n;
		325
		326	if (plan->flags & FFTW_IN_PLACE)
		327	rexecutor_simple_inplace(n, in, out, plan->root, 1);
		328	else
		329	rfftw_executor_simple(n, in, out, plan->root, 1, 1);
		330	}

Subversion Repositories shark

(root)/shark/branches/pj/ports/fftw/rfftw/rexec.c - Rev 2