planner_wisdom()，fftw_wisdom_lookup()，fftw_measure_runtime()，init_test_array()函数代码分析

代码分析以fftw2.15为例，原代码在fftw/planner.c中

planner_wisdom(）函数是fftw为了运行效率提出的wisdom机制，主要思想是通过查找之前相似数据（结构、大小等相似）傅里✌变换时已经计算好的plan来节省时间。

static fftw_plan planner_wisdom(fftw_plan *table, int n,
				fftw_direction dir, int flags,
				int vector_size,
				fftw_complex *in, int istride,
				fftw_complex *out, int ostride)
{
     fftw_plan best = (fftw_plan) 0;
     fftw_plan_node *node;
     int have_wisdom;
     enum fftw_node_type wisdom_type;
     int wisdom_signature;
     fftw_recurse_kind wisdom_recurse_kind;

     /* see if we remember any wisdom for this case */
     have_wisdom = fftw_wisdom_lookup(n, flags, dir, FFTW_WISDOM,
				      istride, ostride,
				      &wisdom_type, &wisdom_signature,
				      &wisdom_recurse_kind, 0);
     //如果没有相关配置方案
     if (!have_wisdom)
	  return best;
     //FFTW_NOTW和FFTW_TWIDDLE模式下的方案生成
     if (wisdom_type == FFTW_NOTW) {
	  FOR_ALL_CODELETS(p) {
	       if (p->dir == dir && p->type == wisdom_type) {
		    /* see if wisdom applies */
		    if (wisdom_signature == p->signature &&
			p->size == n) {
			 //创建对应方案的节点
			 node = fftw_make_node_notw(n, p);
			 //创建执行方案
			 best = fftw_make_plan(n, dir, node, flags,
					       p->type, p->signature,
					       FFTW_NORMAL_RECURSE,
					       vector_size);
			 fftw_use_plan(best);
			 run_plan_hooks(best);
			 return best;
		    }
	       }
	  }
     }
     if (wisdom_type == FFTW_TWIDDLE) {
	  FOR_ALL_CODELETS(p) {
	       if (p->dir == dir && p->type == wisdom_type) {

		    /* see if wisdom applies */
		    if (wisdom_signature == p->signature &&
			p->size > 1 &&
			(n % p->size) == 0) {
			 fftw_plan r = planner(table, n / p->size, dir, 
					       flags | FFTW_NO_VECTOR_RECURSE,
					       wisdom_recurse_kind ==
					       FFTW_VECTOR_RECURSE ?
					       p->size : vector_size,
					       in, istride, out, ostride);
			 node = fftw_make_node_twiddle(n, p,
						       r->root, flags);
			 best = fftw_make_plan(n, dir, node, flags,
					       p->type, p->signature,
					       wisdom_recurse_kind, 
					       vector_size);
			 fftw_use_plan(best);
			 run_plan_hooks(best);
			 fftw_destroy_plan_internal(r);
			 return best;
		    }
	       }
	  }
     }
     /* 
      * BUG (or: TODO)  Can we have generic wisdom? This is probably
      * an academic question
      */

     return best;
}

int fftw_wisdom_lookup(int n, int flags, fftw_direction dir,
		       enum fftw_wisdom_category category,
		       int istride, int ostride,
		       enum fftw_node_type *type,
		       int *signature, fftw_recurse_kind *recurse_kind,
		       int replacep)
{
     struct wisdom *p;

     if (!(flags & FFTW_USE_WISDOM))
	  return 0;		/* simply ignore if wisdom is disabled */

     flags |= FFTW_MEASURE;	/* 
				 * always use (only) wisdom from
				 * measurements 
				 */
     //遍历链表，查找是否有相应的wisdom
     for (p = wisdom_list; p; p = p->next) {
	  if (p->n == n && p->flags == flags && p->dir == dir &&
	      p->istride == istride && p->ostride == ostride &&
	      p->category == category) {
	       /* found wisdom */
	       if (replacep) {
		    /* replace old wisdom with new */
		    p->type = *type;
		    p->signature = *signature;
		    p->recurse_kind = *recurse_kind;
	       } else {
		    *type = p->type;
		    *signature = p->signature;
		    *recurse_kind = p->recurse_kind;
	       }
	       return 1;
	  }
     }

     return 0;
}

 

fftw_measure_runtime()用来评估一种变换方案的执行时间。

static double fftw_measure_runtime(fftw_plan plan,
				   fftw_complex *in, int istride,
				   fftw_complex *out, int ostride)
{
     fftw_time begin, end, start;
     double t, tmax, tmin;
     int i, iter;
     int n;
     int repeat;
     int howmany = plan->vector_size;

     n = plan->n;

     iter = 1;
     /*
     下面是计算耗时的程序
     具体计算方式为：对每一个iter重复FFTW_TIME_REPEAT次，如果计算时间过短，
     可能不够准确，这时通过将iter扩大二倍再进行计算，当tmin大于某一阈值（FFTW_TIME_MIN）时，
     将此时的结果除以iter后返回
     */
     for (;;) {
	  tmin = 1.0E10;
	  tmax = -1.0E10;
	  init_test_array(in, istride, n * howmany);

	  start = fftw_get_time();
	  /* repeat the measurement FFTW_TIME_REPEAT times */
	  //运行FFTW_TIME_REPEAT次
	  for (repeat = 0; repeat < FFTW_TIME_REPEAT; ++repeat) {
	       begin = fftw_get_time();
	       for (i = 0; i < iter; ++i) {
		    fftw(plan, howmany, in, istride, istride,
			 out, ostride, ostride);
	       }
	       end = fftw_get_time();

	       t = fftw_time_to_sec(fftw_time_diff(end, begin));
	       if (t < tmin)
		    tmin = t;
	       if (t > tmax)
		    tmax = t;

	       /* do not run for too long */
		   //超时退出
	       t = fftw_time_to_sec(fftw_time_diff(end, start));
	       if (t > FFTW_TIME_LIMIT)
		    break;
	  }
	  //
	  if (tmin >= FFTW_TIME_MIN)
	       break;

	  iter *= 2;
     }

     tmin /= (double) iter;
     tmax /= (double) iter;

     return tmin;
}

init_test_array()函数初始化测试变换数据，这里拿出来只是因为其代表了C语言初始化数组的一种经典方式。

static void init_test_array(fftw_complex *arr, int stride, int n)
{
     int j;
     //初始化测试数组全部为零
     for (j = 0; j < n; ++j) {
	  c_re(arr[stride * j]) = 0.0;
	  c_im(arr[stride * j]) = 0.0;
     }
}