wolf.c

shpf 1.9一个并行编译器

#include "stdio.h"

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

#include "fft2d.h"


void read_image(char* filename, Array2<int> image,
                int* nx, int* ny, int* maxin) ;

void write_image(char* filename, Array2<int> image,
                 int nx, int ny, int maxin) ;

void fft_slice(SArray1<complex> a, int N, int isgn) ;


// _____________________________________________________________________
//|                                                                     |
//| Driver for 2 dimensional FFT program, in HPF.                       |
//|                                                                     |
//| H W Yau.  25th of January, 1994.                                    |
//| Northeast Parallel Architectures Center.                            |
//| Syracuse University.                                                |
//|_____________________________________________________________________|
//|                                                                     |
//| Edit record:                                                        |
//| 25/Jan/1995: First edit, code modified from Bhaven Avalani's.       |
//| 31/Jan/1995: Rewrote, to get rid of superfluous call to fft1d.      |
//| 02/Feb/1995: Included HPF constructs and directives.                |
//| 10/Feb/1995: Rewrote so fft_slice only sees a 1-D array.            |
//| 15/Feb/1995: Plugged in routines to FFT in a .pgm image.            |
//| 16/Feb/1995: Corrected Fourier mode removal code, now removes the   |
//|              central portion of the Fourier `image'.                |
//| 20/Feb/1995: Separated codes.  This is now FFT of a constant image. |
//| 17/Jun/1996: Added timer routine.                                   |
//| 10/Jul/1996: Added communications & computation times.              |
//| 20/Aug/1996: Added PURE attribute to fft_slice for PGI compiler to  |
//|              do INDEPENDENT parallelisation.                        |
//| 23/Sep/1996: Converted to C++/Adlib.  DBC.                          |
//|_____________________________________________________________________|


void main(int argc, char** argv) {
  AdlibInit(argc, argv) ;
  int nprocs = AdlibNProcs() ;

  Procs1 p(nprocs) ;

  on(p) {
    Range x(NN) ;
    Range y(NN, p.dim(0), BLK) ;

    Array2<int> image(x, y) ;
    Array2<complex> a(x, y) ;
    Array2<complex> tmp1(x, y) ;

    int sign, itrunc ;
    int nx, ny, maxin ;
    float norm ;

    float ti, tf ;
    float tstart, tend ;
    float tcomms, tcomp ;

    int n ;

// Read in from standard in, the number of Fourier modes to keep (in each
// dimension.

    gprintf("Please enter how many Fourier modes to remove"
	    " in each dimension.\n") ;
    gscanf("%d", &itrunc) ;

// Read in the image, and copy into array to be FFT-ed.

    read_image("wolf.pgm", image, &nx, &ny, &maxin) ;

    where(x) {
      where(y) {
	a(x, y) = complex(image(x, y)) ;
      } erewh(y) ;
    } erewh(x) ;

    norm = 1.0 / (nx * ny) ;

// Set parameter to tell fft_slice to do FFT transform.

    sign = 1 ;
 

    where(y)
      fft_slice(a.sect(x, y.operator Subscript&()), NN, sign) ;
      //fft_slice(SArray1<complex>(a.sect(x, Subscript(y))), NN, sign) ;
    erewh(y) ;

//______________________________________________________________________
//
//.. Transpose
//______________________________________________________________________

    transpose(tmp1, a) ;
    copy(a, tmp1) ;

//______________________________________________________________________
//

    where(y)
      fft_slice(a.sect(x, y.operator Subscript&()), NN, sign) ;
      //fft_slice(SArray1<complex> (a.sect(x, Subscript(y))), NN, sign) ;
    erewh(y) ;

//______________________________________________________________________
//
// Scale magnitude.

    where(x) {
      where(y) {
	a(x, y) *= norm ;
      } erewh(y) ;
    } erewh(x) ;

// Throw away modes in the middle of the frequency range.

    int ioff = (nx - itrunc) / 2 ;

    Range xm = x.sub(itrunc, ioff) ;
    Range ym = y.sub(itrunc, ioff) ;
    where(xm) {
      where(ym) {
	a(xm, ym) = 0.0 ;
      } erewh(ym) ;
    } erewh(xm) ;

//// Throw away Low frequency mode...
//
//    Range xl = x.sub(itrunc, 0) ;
//    Range yl = y.sub(itrunc, 0) ;
//    where(xl) {
//      where(yl) {
//	a(xl, yl) = 0.0 ;
//      } erewh(yl) ;
//    } erewh(xl) ;
//
//// Throw away High frequency modes
//
//    Range yh = y.sub(itrunc, ny - itrunc + 1) ;
//    where(x) {
//      where(yh) {
//	a(x, yh) = 0.0 ;
//      } erewh(yh) ;
//    } erewh(x) ;
//
//    Range xh = x.sub(itrunc, nx - itrunc + 1) ;
//    where(xh) {
//      where(y) {
//	a(xh, y) = 0.0 ;
//      } erewh(y) ;
//    } erewh(xh) ;


//_____________________________________________________________________
//
// Do Inverse transform.
//_____________________________________________________________________

    sign = -1 ;

    where(y)
      fft_slice(a.sect(x, y.operator Subscript&()), NN, sign) ;
      //fft_slice(SArray1<complex>(a.sect(x, Subscript(y))), NN, sign) ;
    erewh(y) ;

//______________________________________________________________________
//
//.. Transpose
//______________________________________________________________________

    transpose(tmp1, a) ;
    copy(a, tmp1) ;

//______________________________________________________________________
//

    where(y)
      fft_slice(a.sect(x, y.operator Subscript&()), NN, sign) ;
      //fft_slice(SArray1<complex>(a.sect(x, Subscript(y))), NN, sign) ;
    erewh(y) ;

//______________________________________________________________________
//
// Output the image.
//______________________________________________________________________


    where(x) {
      where(y) {
	image(x, y) = (int) a(x, y).abs() ;
      } erewh(y) ;
    } erewh(x) ;
    maxin = maxval(image) ;


    write_image("fft_out.pgm", image, nx, ny, maxin) ;
//
//    gprintf("Time taken is %12.3f seconds\n", tf) ;
//    gprintf("Number of Processors = %4i\n"
//            "Problem size = %5i\n"
//            "Communications = %9.3f\n"
//            "Compute        = %9.3f\n"
//            "Others         = %9.3f\n"
//            "Total time     = %9.3f\n",
//            nprocs, NN, tcomms, tcomp, (tf - tcomms - tcomp), tf) ;
  } no(p) ;

  AdlibFinalize() ;
}