WebSVN - shark - Blame - Rev 39 - /shark/branches/pj/ports/fftw/fftw/putils.c

Rev	Author	Line No.	Line
2	pj	1
		2	/*
		3	* Copyright (c) 1997-1999 Massachusetts Institute of Technology
		4	*
		5	* This program is free software; you can redistribute it and/or modify
		6	* it under the terms of the GNU General Public License as published by
		7	* the Free Software Foundation; either version 2 of the License, or
		8	* (at your option) any later version.
		9	*
		10	* This program is distributed in the hope that it will be useful,
		11	* but WITHOUT ANY WARRANTY; without even the implied warranty of
		12	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
		13	* GNU General Public License for more details.
		14	*
		15	* You should have received a copy of the GNU General Public License
		16	* along with this program; if not, write to the Free Software
		17	* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
		18	*
		19	*/
		20
		21	/*
		22	* putils.c -- plan utilities shared by planner.c and rplanner.c
		23	*/
		24
		25	/* $Id: putils.c,v 1.1.1.1 2002-03-29 14:12:57 pj Exp $ */
		26	#ifdef FFTW_USING_CILK
		27	#include <ports/cilk.h>
		28	#include <ports/cilk-compat.h>
		29	#endif
		30
		31	#include <ports/fftw-int.h>
		32	#include <stdlib.h>
		33	//#include <ports/stdio.h>
		34
		35	int fftw_node_cnt = 0;
		36	int fftw_plan_cnt = 0;
		37
		38	/*
		39	* These two constants are used for the FFTW_ESTIMATE flag to help
		40	* create a heuristic plan. They don't affect FFTW_MEASURE.
		41	*/
		42	#define NOTW_OPTIMAL_SIZE 32
		43	#define TWIDDLE_OPTIMAL_SIZE 12
		44
		45	#define IS_POWER_OF_TWO(n) ((n & (n - 1)) == 0)
		46
		47	/* constructors --- I wish I had ML */
		48	fftw_plan_node *fftw_make_node(void)
		49	{
		50	fftw_plan_node p = (fftw_plan_node )
		51	fftw_malloc(sizeof(fftw_plan_node));
		52	p->refcnt = 0;
		53	fftw_node_cnt++;
		54	return p;
		55	}
		56
		57	void fftw_use_node(fftw_plan_node *p)
		58	{
		59	++p->refcnt;
		60	}
		61
		62	fftw_plan_node fftw_make_node_notw(int size, const fftw_codelet_desc config)
		63	{
		64	fftw_plan_node *p = fftw_make_node();
		65
		66	p->type = config->type;
		67	p->nodeu.notw.size = size;
		68	p->nodeu.notw.codelet = (fftw_notw_codelet *) config->codelet;
		69	p->nodeu.notw.codelet_desc = config;
		70	return p;
		71	}
		72
		73	fftw_plan_node *fftw_make_node_real2hc(int size,
		74	const fftw_codelet_desc *config)
		75	{
		76	fftw_plan_node *p = fftw_make_node();
		77
		78	p->type = config->type;
		79	p->nodeu.real2hc.size = size;
		80	p->nodeu.real2hc.codelet = (fftw_real2hc_codelet *) config->codelet;
		81	p->nodeu.real2hc.codelet_desc = config;
		82	return p;
		83	}
		84
		85	fftw_plan_node *fftw_make_node_hc2real(int size,
		86	const fftw_codelet_desc *config)
		87	{
		88	fftw_plan_node *p = fftw_make_node();
		89
		90	p->type = config->type;
		91	p->nodeu.hc2real.size = size;
		92	p->nodeu.hc2real.codelet = (fftw_hc2real_codelet *) config->codelet;
		93	p->nodeu.hc2real.codelet_desc = config;
		94	return p;
		95	}
		96
		97	fftw_plan_node *fftw_make_node_twiddle(int n,
		98	const fftw_codelet_desc *config,
		99	fftw_plan_node *recurse,
		100	int flags)
		101	{
		102	fftw_plan_node *p = fftw_make_node();
		103
		104	p->type = config->type;
		105	p->nodeu.twiddle.size = config->size;
		106	p->nodeu.twiddle.codelet = (fftw_twiddle_codelet *) config->codelet;
		107	p->nodeu.twiddle.recurse = recurse;
		108	p->nodeu.twiddle.codelet_desc = config;
		109	fftw_use_node(recurse);
		110	if (flags & FFTW_MEASURE)
		111	p->nodeu.twiddle.tw = fftw_create_twiddle(n, config);
		112	else
		113	p->nodeu.twiddle.tw = 0;
		114	return p;
		115	}
		116
		117	fftw_plan_node *fftw_make_node_hc2hc(int n, fftw_direction dir,
		118	const fftw_codelet_desc *config,
		119	fftw_plan_node *recurse,
		120	int flags)
		121	{
		122	fftw_plan_node *p = fftw_make_node();
		123
		124	p->type = config->type;
		125	p->nodeu.hc2hc.size = config->size;
		126	p->nodeu.hc2hc.dir = dir;
		127	p->nodeu.hc2hc.codelet = (fftw_hc2hc_codelet *) config->codelet;
		128	p->nodeu.hc2hc.recurse = recurse;
		129	p->nodeu.hc2hc.codelet_desc = config;
		130	fftw_use_node(recurse);
		131	if (flags & FFTW_MEASURE)
		132	p->nodeu.hc2hc.tw = fftw_create_twiddle(n, config);
		133	else
		134	p->nodeu.hc2hc.tw = 0;
		135	return p;
		136	}
		137
		138	fftw_plan_node *fftw_make_node_generic(int n, int size,
		139	fftw_generic_codelet *codelet,
		140	fftw_plan_node *recurse,
		141	int flags)
		142	{
		143	fftw_plan_node *p = fftw_make_node();
		144
		145	p->type = FFTW_GENERIC;
		146	p->nodeu.generic.size = size;
		147	p->nodeu.generic.codelet = codelet;
		148	p->nodeu.generic.recurse = recurse;
		149	fftw_use_node(recurse);
		150
		151	if (flags & FFTW_MEASURE)
		152	p->nodeu.generic.tw = fftw_create_twiddle(n,
		153	(const fftw_codelet_desc *) 0);
		154	else
		155	p->nodeu.generic.tw = 0;
		156	return p;
		157	}
		158
		159	fftw_plan_node *fftw_make_node_rgeneric(int n, int size,
		160	fftw_direction dir,
		161	fftw_rgeneric_codelet *codelet,
		162	fftw_plan_node *recurse,
		163	int flags)
		164	{
		165	fftw_plan_node *p = fftw_make_node();
		166
		167	if (size % 2 == 0 \|\| (n / size) % 2 == 0)
		168	fftw_die("invalid size for rgeneric codelet\n");
		169
		170	p->type = FFTW_RGENERIC;
		171	p->nodeu.rgeneric.size = size;
		172	p->nodeu.rgeneric.dir = dir;
		173	p->nodeu.rgeneric.codelet = codelet;
		174	p->nodeu.rgeneric.recurse = recurse;
		175	fftw_use_node(recurse);
		176
		177	if (flags & FFTW_MEASURE)
		178	p->nodeu.rgeneric.tw = fftw_create_twiddle(n,
		179	(const fftw_codelet_desc *) 0);
		180	else
		181	p->nodeu.rgeneric.tw = 0;
		182	return p;
		183	}
		184
		185	/*
		186	* Note that these two Rader-related things must go here, rather than
		187	* in rader.c, in order that putils.c (and rplanner.c) won't depend
		188	* upon rader.c.
		189	*/
		190
		191	fftw_rader_data *fftw_rader_top = NULL;
		192
		193	static void fftw_destroy_rader(fftw_rader_data * d)
		194	{
		195	if (d) {
		196	d->refcount--;
		197	if (d->refcount <= 0) {
		198	fftw_rader_data cur = fftw_rader_top, prev = NULL;
		199
		200	while (cur && cur != d) {
		201	prev = cur;
		202	cur = cur->next;
		203	}
		204	if (!cur)
		205	fftw_die("invalid Rader data pointer\n");
		206
		207	if (prev)
		208	prev->next = d->next;
		209	else
		210	fftw_rader_top = d->next;
		211
		212	fftw_destroy_plan_internal(d->plan);
		213	fftw_free(d->omega);
		214	fftw_free(d->cdesc);
		215	fftw_free(d);
		216	}
		217	}
		218	}
		219
		220	static void destroy_tree(fftw_plan_node *p)
		221	{
		222	if (p) {
		223	--p->refcnt;
		224	if (p->refcnt == 0) {
		225	switch (p->type) {
		226	case FFTW_NOTW:
		227	case FFTW_REAL2HC:
		228	case FFTW_HC2REAL:
		229	break;
		230
		231	case FFTW_TWIDDLE:
		232	if (p->nodeu.twiddle.tw)
		233	fftw_destroy_twiddle(p->nodeu.twiddle.tw);
		234	destroy_tree(p->nodeu.twiddle.recurse);
		235	break;
		236
		237	case FFTW_HC2HC:
		238	if (p->nodeu.hc2hc.tw)
		239	fftw_destroy_twiddle(p->nodeu.hc2hc.tw);
		240	destroy_tree(p->nodeu.hc2hc.recurse);
		241	break;
		242
		243	case FFTW_GENERIC:
		244	if (p->nodeu.generic.tw)
		245	fftw_destroy_twiddle(p->nodeu.generic.tw);
		246	destroy_tree(p->nodeu.generic.recurse);
		247	break;
		248
		249	case FFTW_RADER:
		250	if (p->nodeu.rader.tw)
		251	fftw_destroy_twiddle(p->nodeu.rader.tw);
		252	if (p->nodeu.rader.rader_data)
		253	fftw_destroy_rader(p->nodeu.rader.rader_data);
		254	destroy_tree(p->nodeu.rader.recurse);
		255	break;
		256
		257	case FFTW_RGENERIC:
		258	if (p->nodeu.rgeneric.tw)
		259	fftw_destroy_twiddle(p->nodeu.rgeneric.tw);
		260	destroy_tree(p->nodeu.rgeneric.recurse);
		261	break;
		262	}
		263
		264	fftw_free(p);
		265	fftw_node_cnt--;
		266	}
		267	}
		268	}
		269
		270	/* create a plan with twiddle factors, and other bells and whistles */
		271	fftw_plan fftw_make_plan(int n, fftw_direction dir,
		272	fftw_plan_node *root, int flags,
		273	enum fftw_node_type wisdom_type,
		274	int wisdom_signature)
		275	{
		276	fftw_plan p = (fftw_plan) fftw_malloc(sizeof(struct fftw_plan_struct));
		277
		278	p->n = n;
		279	p->dir = dir;
		280	p->flags = flags;
		281	fftw_use_node(root);
		282	p->root = root;
		283	p->cost = 0.0;
		284	p->wisdom_type = wisdom_type;
		285	p->wisdom_signature = wisdom_signature;
		286	p->next = (fftw_plan) 0;
		287	p->refcnt = 0;
		288	fftw_plan_cnt++;
		289	return p;
		290	}
		291
		292	/*
		293	* complete with twiddle factors (because nodes don't have
		294	* them when FFTW_ESTIMATE is set)
		295	*/
		296	void fftw_complete_twiddle(fftw_plan_node *p, int n)
		297	{
		298	int r;
		299	switch (p->type) {
		300	case FFTW_NOTW:
		301	case FFTW_REAL2HC:
		302	case FFTW_HC2REAL:
		303	break;
		304
		305	case FFTW_TWIDDLE:
		306	r = p->nodeu.twiddle.size;
		307	if (!p->nodeu.twiddle.tw)
		308	p->nodeu.twiddle.tw =
		309	fftw_create_twiddle(n, p->nodeu.twiddle.codelet_desc);
		310	fftw_complete_twiddle(p->nodeu.twiddle.recurse, n / r);
		311	break;
		312
		313	case FFTW_HC2HC:
		314	r = p->nodeu.hc2hc.size;
		315	if (!p->nodeu.hc2hc.tw)
		316	p->nodeu.hc2hc.tw =
		317	fftw_create_twiddle(n, p->nodeu.hc2hc.codelet_desc);
		318	fftw_complete_twiddle(p->nodeu.hc2hc.recurse, n / r);
		319	break;
		320
		321	case FFTW_GENERIC:
		322	r = p->nodeu.generic.size;
		323	if (!p->nodeu.generic.tw)
		324	p->nodeu.generic.tw =
		325	fftw_create_twiddle(n, (const fftw_codelet_desc *) 0);
		326	fftw_complete_twiddle(p->nodeu.generic.recurse, n / r);
		327	break;
		328
		329	case FFTW_RADER:
		330	r = p->nodeu.rader.size;
		331	if (!p->nodeu.rader.tw)
		332	p->nodeu.rader.tw =
		333	fftw_create_twiddle(n, p->nodeu.rader.rader_data->cdesc);
		334	fftw_complete_twiddle(p->nodeu.rader.recurse, n / r);
		335	break;
		336
		337	case FFTW_RGENERIC:
		338	r = p->nodeu.rgeneric.size;
		339	if (!p->nodeu.rgeneric.tw)
		340	p->nodeu.rgeneric.tw =
		341	fftw_create_twiddle(n, (const fftw_codelet_desc *) 0);
		342	fftw_complete_twiddle(p->nodeu.rgeneric.recurse, n / r);
		343	break;
		344
		345	}
		346	}
		347
		348	void fftw_use_plan(fftw_plan p)
		349	{
		350	++p->refcnt;
		351	}
		352
		353	void fftw_destroy_plan_internal(fftw_plan p)
		354	{
		355	--p->refcnt;
		356
		357	if (p->refcnt == 0) {
		358	destroy_tree(p->root);
		359	fftw_plan_cnt--;
		360	fftw_free(p);
		361	}
		362	}
		363
		364	/* end of constructors */
		365
		366	/* management of plan tables */
		367	void fftw_make_empty_table(fftw_plan *table)
		368	{
		369	*table = (fftw_plan) 0;
		370	}
		371
		372	void fftw_insert(fftw_plan *table, fftw_plan this_plan, int n)
		373	{
		374	fftw_use_plan(this_plan);
		375	this_plan->n = n;
		376	this_plan->next = *table;
		377	*table = this_plan;
		378	}
		379
		380	fftw_plan fftw_lookup(fftw_plan *table, int n, int flags)
		381	{
		382	fftw_plan p;
		383
		384	for (p = *table; p &&
		385	((p->n != n) \|\| (p->flags != flags)); p = p->next);
		386
		387	return p;
		388	}
		389
		390	void fftw_destroy_table(fftw_plan *table)
		391	{
		392	fftw_plan p, q;
		393
		394	for (p = *table; p; p = q) {
		395	q = p->next;
		396	fftw_destroy_plan_internal(p);
		397	}
		398	}
		399
		400	double fftw_estimate_node(fftw_plan_node *p)
		401	{
		402	int k;
		403
		404	switch (p->type) {
		405	case FFTW_NOTW:
		406	k = p->nodeu.notw.size;
		407	goto common1;
		408
		409	case FFTW_REAL2HC:
		410	k = p->nodeu.real2hc.size;
		411	goto common1;
		412
		413	case FFTW_HC2REAL:
		414	k = p->nodeu.hc2real.size;
		415	common1:
		416	return 1.0 + 0.1 * (k - NOTW_OPTIMAL_SIZE) *
		417	(k - NOTW_OPTIMAL_SIZE);
		418
		419	case FFTW_TWIDDLE:
		420	k = p->nodeu.twiddle.size;
		421	return 1.0 + 0.1 * (k - TWIDDLE_OPTIMAL_SIZE) *
		422	(k - TWIDDLE_OPTIMAL_SIZE)
		423	+ fftw_estimate_node(p->nodeu.twiddle.recurse);
		424
		425	case FFTW_HC2HC:
		426	k = p->nodeu.hc2hc.size;
		427	return 1.0 + 0.1 * (k - TWIDDLE_OPTIMAL_SIZE) *
		428	(k - TWIDDLE_OPTIMAL_SIZE)
		429	+ fftw_estimate_node(p->nodeu.hc2hc.recurse);
		430
		431	case FFTW_GENERIC:
		432	k = p->nodeu.generic.size;
		433	return 10.0 + k * k
		434	+ fftw_estimate_node(p->nodeu.generic.recurse);
		435
		436	case FFTW_RADER:
		437	k = p->nodeu.rader.size;
		438	return 10.0 + 10 * k
		439	+ fftw_estimate_node(p->nodeu.rader.recurse);
		440
		441	case FFTW_RGENERIC:
		442	k = p->nodeu.rgeneric.size;
		443	return 10.0 + k * k
		444	+ fftw_estimate_node(p->nodeu.rgeneric.recurse);
		445	}
		446	return 1.0E20;
		447	}
		448
		449	/* pick the better of two plans and destroy the other one. */
		450	fftw_plan fftw_pick_better(fftw_plan p1, fftw_plan p2)
		451	{
		452	if (!p1)
		453	return p2;
		454
		455	if (!p2)
		456	return p1;
		457
		458	if (p1->cost > p2->cost) {
		459	fftw_destroy_plan_internal(p1);
		460	return p2;
		461	} else {
		462	fftw_destroy_plan_internal(p2);
		463	return p1;
		464	}
		465	}
		466
		467	/* find the smallest prime factor of n */
		468	int fftw_factor(int n)
		469	{
		470	int r;
		471
		472	/* try 2 */
		473	if ((n & 1) == 0)
		474	return 2;
		475
		476	/* try odd numbers up to sqrt(n) */
		477	for (r = 3; r * r <= n; r += 2)
		478	if (n % r == 0)
		479	return r;
		480
		481	/* n is prime */
		482	return n;
		483	}
		484	/*
		485	static void print_node(FILE f, fftw_plan_node p, int indent)
		486	{
		487	if (p) {
		488	switch (p->type) {
		489	case FFTW_NOTW:
		490	fprintf(f, "%*sFFTW_NOTW %d\n", indent, "",
		491	p->nodeu.notw.size);
		492	break;
		493	case FFTW_REAL2HC:
		494	fprintf(f, "%*sFFTW_REAL2HC %d\n", indent, "",
		495	p->nodeu.real2hc.size);
		496	break;
		497	case FFTW_HC2REAL:
		498	fprintf(f, "%*sFFTW_HC2REAL %d\n", indent, "",
		499	p->nodeu.hc2real.size);
		500	break;
		501	case FFTW_TWIDDLE:
		502	fprintf(f, "%*sFFTW_TWIDDLE %d\n", indent, "",
		503	p->nodeu.twiddle.size);
		504	print_node(f, p->nodeu.twiddle.recurse, indent);
		505	break;
		506	case FFTW_HC2HC:
		507	fprintf(f, "%*sFFTW_HC2HC %d\n", indent, "",
		508	p->nodeu.hc2hc.size);
		509	print_node(f, p->nodeu.hc2hc.recurse, indent);
		510	break;
		511	case FFTW_GENERIC:
		512	fprintf(f, "%*sFFTW_GENERIC %d\n", indent, "",
		513	p->nodeu.generic.size);
		514	print_node(f, p->nodeu.generic.recurse, indent);
		515	break;
		516	case FFTW_RADER:
		517	fprintf(f, "%*sFFTW_RADER %d\n", indent, "",
		518	p->nodeu.rader.size);
		519
		520	fprintf(f, "%*splan for size %d convolution:\n",
		521	indent + 4, "", p->nodeu.rader.size - 1);
		522	print_node(f, p->nodeu.rader.rader_data->plan->root,
		523	indent + 6);
		524
		525	print_node(f, p->nodeu.rader.recurse, indent);
		526	break;
		527	case FFTW_RGENERIC:
		528	fprintf(f, "%*sFFTW_RGENERIC %d\n", indent, "",
		529	p->nodeu.rgeneric.size);
		530	print_node(f, p->nodeu.rgeneric.recurse, indent);
		531	break;
		532	}
		533	}
		534	}
		535
		536	void fftw_fprint_plan(FILE *f, fftw_plan p)
		537	{
		538
		539	fprintf(f, "plan: (cost = %e)\n", p->cost);
		540	print_node(f, p->root, 0);
		541	}
		542
		543	void fftw_print_plan(fftw_plan p)
		544	{
		545	fftw_fprint_plan(stdout, p);
		546	}
		547	*/
		548	size_t fftw_sizeof_fftw_real(void)
		549	{
		550	return(sizeof(fftw_real));
		551	}

Subversion Repositories shark

(root)/shark/branches/pj/ports/fftw/fftw/putils.c - Rev 39