「最適化・並列化コードを生み出す最新コンパイラ「インテル コンパイラー」サンプルコード（3/4） - HPC/並列プログラミングポータル

Up: ［記事全文PDFのダウンロード］

Prev: 「最適化・並列化コードを生み出す最新コンパイラ「インテル コンパイラー」サンプルコード（2/4）

Next: 「最適化・並列化コードを生み出す最新コンパイラ「インテル コンパイラー」サンプルコード（4/4）

リスト3 『インテル謹製の数値演算ライブラリ「MKL」を使ってプログラムを高速化』記事内で使用したFFTWによるFFT処理プログラム

/* fft1.c */

#include stdio.h
#include string.h
#include stdlib.h
#include memory.h
#include math.h
#include time.h

#include "fftw3.h"

/*
  fft1 input file output_file
 */

#define FMODE_WRITE "w"
#define NORM(x,y) (sqrt(x*x+y*y))
#define POWR(x,y) (x*x+y*y)


#define c_re(x) (x[0])
#define c_im(x) (x[1])

/* prototypes */
unsigned char* load_image( char* filename,
                           int* pWidth,
                           int* pHeight,
                           int* pDepth );
void write_image( char* filename,
                  unsigned char* buff,
                  int width,
                  int height,
                  int depth );


unsigned char* load_image( char* filename,
                           int* pWidth,
                           int* pHeight,
                           int* pDepth ) {
    FILE* fp_in;
    unsigned char* buff;
    char szBuf[512];
    int height = 0;
    int width = 0;
    int depth = 0;

    fp_in = fopen( filename, "r" );
    if( fp_in == NULL ) {
        fprintf( stderr, "cannot open file: %s\n", filename );
        return NULL;
    }

    /* read PGM Header */
    fgets( szBuf, 512, fp_in );
    if( strlen( szBuf )  3 ) {
        fprintf( stderr, "invalid image file: %s\n", filename );
        fclose( fp_in );
        return NULL;
    }
    if( strcmp( szBuf, "P5\n" ) != 0 ) {
        fprintf( stderr, "invalid image file: %s\n", filename );
        fclose( fp_in );
        return NULL;
    }

    /* read #comment */
    fgets( szBuf, 512, fp_in );
    while( szBuf[0] == '#' ) {
        while( (szBuf[strlen(szBuf)-1] != '\n')  (!feof(fp_in)) ) {
            fgets( szBuf, 512, fp_in );
        }
        fgets( szBuf, 512, fp_in );
    }

    if(feof(fp_in)) {
        fprintf( stderr, "invalid image file: %s\n", filename );
        fclose( fp_in );
        return NULL;
    }

    /* read image size */
    if (szBuf[strlen(szBuf)-1] != '\n') {
        fprintf( stderr, "image file is strange: %s\n", filename );
        fclose( fp_in );
        return NULL;
    }
    sscanf( szBuf, "%d %d", width, height );
    if( (width == 0) || (height == 0) ) {
        fprintf( stderr, "image file is strange: %s\n", filename );
        fclose( fp_in );
        return NULL;
    }
    fgets( szBuf, 512, fp_in );

    /* read depth */
    if (szBuf[strlen(szBuf)-1] != '\n') {
        fprintf( stderr, "image file is strange: %s\n", filename );
        fclose( fp_in );
        return NULL;
    }
    sscanf( szBuf, "%d", depth );
    if( depth == 0 ) {
        fprintf( stderr, "image file is strange: %s\n", filename );
        fclose( fp_in );
        return NULL;
    }

    printf( "image size: %d x %d, depth: %d\n", width, height, depth );

    /* allocate buffer */
    buff = (unsigned char*)malloc( width * height * sizeof(unsigned char) );
    if( buff == NULL ) {
        fclose( fp_in );
        return NULL;
    }

    fread( buff, 1, width*height, fp_in );
    fclose( fp_in );

    *pWidth = width;
    *pHeight = height;
    *pDepth = depth;

    return buff;
}

void write_image( char* filename,
                  unsigned char* buff,
                  int width,
                  int height,
                  int depth ) {
    FILE* fp_out;

    fp_out = fopen( filename, FMODE_WRITE );
    if( fp_out == NULL ) {
        fprintf( stderr, "cannot open file: %s\n", filename );
        return;
    }

    /* write PGM Header */
    fprintf( fp_out, "P5\n" );

    /* write image size */
    fprintf( fp_out, "%d %d\n", width, height );

    /* write depth */
    fprintf( fp_out, "%d\n", depth );
    fwrite( buff, 1, width*height, fp_out );

    if( ferror(fp_out) ) {
        fprintf( stderr, "file output error: %s\n", filename );
        return;
    }

}

void FFTWExchange(int sizex,int sizey, fftw_complex *src, fftw_complex *dst) {
  int i,j,u,v;
  int imgxd2=sizex/2, imgyd2=sizey/2;
  for(j=0;jsizey;j++){
    for(i=0;isizex;i++){
      u = (iimgxd2)? imgxd2+i: i-imgxd2;
      v = (jimgyd2)? imgyd2+j: j-imgyd2;
      c_re(dst[j*sizex+i]) = c_re(src[u+sizey*v]);
      c_im(dst[j*sizex+i]) = c_im(src[u+sizey*v]);
    }
  }
}

unsigned char* execute_fft( unsigned char* input_image_buf,
                            int width,
                            int height,
                            int depth ) {

  fftw_complex *in, *out;
  fftw_plan p;
  unsigned char* output_image_buf;
  int n;
  int i,j;
  double dMax, dMin, dNorm;
  double* pTmp;
  dMax = 0.0;
  dMin = 0.0;

  pTmp = (double*)malloc( width*height*sizeof(double) );
  if( pTmp == NULL ) {
    return NULL;
  }

  /* allocate memory */
  in  = fftw_malloc(sizeof(fftw_complex) * width * height);
  out = fftw_malloc(sizeof(fftw_complex) * width * height);

  /* copy data to fftw buffer */
  for( n = 0; n  width*height; n++ ) {
    in[n][0] = (float)input_image_buf[n];
    in[n][1] = 0.0;
  }

  /* create plan */
  p = fftw_plan_dft_2d(height, width, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

  /* do FFT */
  fftw_execute(p);

  /* allocate buffer */
  output_image_buf =
    (unsigned char*)malloc( width * height * sizeof(unsigned char) );

  if( output_image_buf != NULL ) {
    /* copy data to fftw buffer */
    for( n = 0; n  width*height; n++ ) {
      pTmp[n] = log(POWR(out[n][0],out[n][1]) + 1);
      if( pTmp[n]  dMin ) {
        dMin = pTmp[n];
      }
      if( dMax  pTmp[n] ) {
        dMax = pTmp[n];
      }
    }

   for( n = 0; n  width*height; n++ ) {
     double dWork;
     dWork = (pTmp[n] * 255.0 / (dMax - dMin));
     output_image_buf[n] = (unsigned char)(dWork);
    }
  }

  /* free memories */
  fftwnd_destroy_plan(p);

  free( pTmp );
  fftw_free(in);
  fftw_free(out);

  return output_image_buf;

}

unsigned char* execute_ifft( unsigned char* input_image_buf,
                            int width,
                            int height,
                            int depth ) {

  fftw_complex *in, *out;
  fftw_plan p;
  unsigned char* output_image_buf;
  int n;

  /* allocate memory and create plan */
  in  = fftw_malloc(sizeof(fftw_complex) * width * height);
  out = fftw_malloc(sizeof(fftw_complex) * width * height);

  p   = fftw_plan_dft_2d(height, width, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);

  /* allocate output buffer */


  /* copy data to fftw buffer */
  for( n = 0; n  width*height; n++ ) {
    in[n][0] = input_image_buf[n]; /* Re */
    in[n][1] = 0; /* Im */
  }

  /* do FFT */
  fftw_execute(p);

  /* allocate buffer */
  output_image_buf =
    (unsigned char*)malloc( width * height * sizeof(unsigned char) );

  if( output_image_buf != NULL ) {
    /* copy data to fftw buffer */
    for( n = 0; n  width*height; n++ ) {
      output_image_buf[n] = out[n][0];
    }
  }

  /* free memories */
  /*fftwnd_destroy_plan(p);*/
  fftw_free(in);
  fftw_free(out);

  return output_image_buf;

}



int main( int argv, char** argc ) {
    char* in_file_path;
    char* out_file_path;
    unsigned char* input_image_buf;
    unsigned char* output_image_buf;
    unsigned char* image_buf;
    int width;
    int height;
    int depth;
    clock_t clock_start, clock_end;

    if( argv  3 ) {
        printf( "%s input_file output_file\n", argc[0] );
        return -1;
    }
    in_file_path = argc[1];
    out_file_path = argc[2];

    input_image_buf = load_image( in_file_path, width, height, depth );
    if( input_image_buf == NULL ) {
        fprintf( stderr, "cannot load image. exit.\n" );
        return -1;
    }

    clock_start = clock();

    output_image_buf =
      execute_fft( input_image_buf, width, height, depth );

    clock_end = clock();

    if( output_image_buf != NULL ) {
      write_image( out_file_path, output_image_buf, width, height, depth );
    } else {
      fprintf( stderr, "FFT operation faild.\n" );
    }

    free( input_image_buf );
    free( output_image_buf );

    printf( "start: %d end: %d elapse: %d.\n", clock_start, clock_end, clock_end - clock_start );
    printf( "elapse time: %lf sec.", ((double)(clock_end - clock_start)) / (double)CLOCKS_PER_SEC );

    /*    free( image_buf ); */
    return 0;
};

Up: ［記事全文PDFのダウンロード］

Prev: 「最適化・並列化コードを生み出す最新コンパイラ「インテル コンパイラー」サンプルコード（2/4）

Next: 「最適化・並列化コードを生み出す最新コンパイラ「インテル コンパイラー」サンプルコード（4/4）