fft_slice.c

shpf 1.9一个并行编译器

#include "ad++.h"
#include "admacros.h"
#include "adftn.h"

int ipow(int i, int j) ;
int ilog2(int n) ;



void fft_slice(SArray1<complex> a, int N, int isgn)
// a.rng(0) == N
{

// _____________________________________________________________________
//|                                                                     |
//| Subroutine, to perform FFT over one dimension of a square matrix.   |
//| Namely, the leading dimension of the array `a', column `icol'.      |
//| isgn is an integer which determines whether to do an FFT or the     |
//| inverse operation.                                                  |
//|                                                                     |
//| H W Yau.  26th of January, 1995.                                    |
//| Northeast Parallel Architectures Center.                            |
//| Syracuse University.                                                |
//|_____________________________________________________________________|
//|                                                                     |
//| Edit record:                                                        |
//| 31/Jan/1995: First edit, code modified from Numerical Recipes and   |
//|              Bhaven Avalani.                                        |
//| 02/Feb/1995: Included HPF constructs and directives.                |
//| 06/Feb/1995: Minor change to allow the inverse transform.           |
//| 10/Feb/1995: Rewrote using explict array slices, ie minus HPF array |
//|              decomposition directives.             	                |
//| 16/Feb/1995: Removed PURE declaration of ilog2 & ipow functions,    |
//|              PGI does not like it.                                  |
//| 20/Aug/1996: Added PURE attribute to `fft_slice' & extrinsic        |
//|              functions for PGI compiler to do INDEPENDENT           |
//|              parallelisation.                                       |
//|_____________________________________________________________________|

  int m, ln2, nv2, nm1 ;
  int index, jndex ;
  int nrows, ncols, le, lev2 ;
  int ilevel, ii, jj ;
  complex tmp, u, w ;
  const float pi = 3.141592654 ;

  nrows = N ;
  ncols = N ;
  nv2   = nrows / 2 ;
  nm1   = nrows - 1 ;
  ln2   = ilog2(N)  ;  // Base to the 2 log of the leading dimension.

// Sort by bit-reversed ordering.

  jndex = 0 ;
  for(index = 0 ; index < nm1 ; index++) {
    if(jndex > index) {
      tmp      = a(jndex) ;  // Swap entries
      a(jndex) = a(index) ;
      a(index) = tmp ;
    }
    m = nv2 ;
    while ((m >= 2) && (jndex >= m)) {
      jndex = jndex - m ;
      m = m / 2 ;
    }
    jndex = jndex + m ;
  }

// Danielson-Lanczos algorithm for FFT.

  for(ilevel = 1 ; ilevel <= ln2 ; ilevel++) {
    le   = ipow(2,ilevel) ;
    lev2 = le / 2 ;
    u = complex(1.0, 0.0) ;
    w = complex(cos(isgn * pi / lev2), sin(isgn * pi / lev2)) ;
    for(jj = 0 ; jj < lev2 ; jj++) {
      for(ii = jj ; ii < nrows ; ii += le) {
        jndex = ii + lev2 ;
        index = ii ;
        tmp      = u * a(jndex) ;
        a(jndex) = a(index) - tmp ;
        a(index) = a(index) + tmp ;
      }
      tmp = u * w ;
      u   = tmp ;
    }
  }
}

//======================================================================