serial.c

#include <assert.h>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/times.h>
#include <sys/time.h>
#include <time.h>
#define N 8

typedef struct {float r; float i;} complex;
static complex ctmp;

#define C_SWAP(a,b) {ctmp=(a);(a)=(b);(b)=ctmp;}

void print_array(complex *array)
{
  int row, col;

  if (N < 17)
  {
    printf("\nA =\n\t");
    for (row = 0; row < N; row++)
    {
      for (col = 0; col < N; col++)
      {
        printf("%5.2f%s", array[row*N+col].r, (col < N - 1) ? ", " : ";\n\t");
      }
    }
    printf("\n");
  }
}

void readFile(char fileName[15], complex *array){
    FILE *fp = fopen(fileName, "r");

    int i, j, result;

    for (i=0;i<N;i++){
        for (j=0;j<N;j++){
            result = fscanf(fp,"%g",&array[i*N+j].r);
            array[i*N+j].i = 0.00;
        }
    }

    fclose(fp);
}

void writeToFile(char fileName[15], complex *array){
    FILE *fp = fopen(fileName, "w");
    int i, j;
    for (i=0;i<N;i++) {
        for (j=0;j<N;j++){
            fprintf(fp,"   %.7e",array[i*N+j].r);
        }
        fprintf(fp,"\n");
    }
    fclose(fp);
}

void transpose(complex *r, int n ) {
    for ( int i = 0; i < (n-1); i++ ) {
        for ( int j =i+1; j < n; j++ ) {
            C_SWAP(r[i*n+j],r[j*n+i]);
        }
    }
}

void c_fft1d(complex *r, int n, int isign)
{
   int     m,i,i1,j,k,i2,l,l1,l2;
   float   c1,c2,z;
   complex t, u;

   if (isign == 0) return;

   /* Do the bit reversal */
   i2 = n >> 1;
   j = 0;
   for (i=0;i<n/2;i+=2) {
      if (i < j)
         C_SWAP(r[i], r[j]);
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   /* m = (int) log2((double)n); */
   for (i=n,m=0; i>1; m++,i/=2);

   /* Compute the FFT */
   c1 = -1.0;
   c2 =  0.0;
   l2 =  1;
   for (l=0;l<m;l++) {
      l1   = l2;
      l2 <<= 1;
      u.r = 1.0;
      u.i = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<n;i+=l2) {
            i1 = i + l1;

            /* t = u * r[i1] */
            t.r = u.r * r[i1].r - u.i * r[i1].i;
            t.i = u.r * r[i1].i + u.i * r[i1].r;

            /* r[i1] = r[i] - t */
            r[i1].r = r[i].r - t.r;
            r[i1].i = r[i].i - t.i;

            /* r[i] = r[i] + t */
            r[i].r += t.r;
            r[i].i += t.i;
         }
         z =  u.r * c1 - u.i * c2;

         u.i = u.r * c2 + u.i * c1;
         u.r = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (isign == -1) /* FWD FFT */
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   /* Scaling for inverse transform */
   if (isign == 1) {       /* IFFT*/
      for (i=0;i<n;i++) {
         r[i].r /= n;
         r[i].i /= n;
      }
   }
}

void c_fft2d(complex *r, int n, int isign, int print ) {
    int i;
    
    for (i = 0; i < N ; i++) {
        c_fft1d(&r[i*N], N, isign);
    }

    if ( print == 1 ) {
        print_array(r);
    }

    transpose(r, N);
    if ( print == 1 ) {
        print_array(r);
    }

    for (i = 0; i < N ; i++) {
        c_fft1d(&r[i*N], N, isign);
    }
    if ( print == 1 ) {
        print_array(r);
    }
}

void mm(complex *result, complex *a, complex *b, int n) {
    int i,j;
    for (i = 0; i < n ; i++) {
        for(j = 0; j < n ; j++) {
            int index = i*n+j;
            result[index].r = (a[index].r * b[index].r) - (a[index].i * b[index].i);
            result[index].i = (a[index].r * b[index].i) + (a[index].i * b[index].r);
        }
    }
}

int main(int argc, char **argv){
    /* Timing variables */
    struct timeval etstart, etstop;  /* Elapsed times using gettimeofday() */
    struct timezone tzdummy;
    clock_t etstart2, etstop2;  /* Elapsed times using times() */
    unsigned long long usecstart, usecstop;
    struct tms cputstart, cputstop;  /* CPU times for my processes */

    complex *A, *B, *C;

    A = (complex *)malloc(sizeof(complex)*N*N);
    B = (complex *)malloc(sizeof(complex)*N*N);
    C = (complex *)malloc(sizeof(complex)*N*N);

    char inputA[15] = "sample/3_im1";
    char inputB[15] = "sample/3_im2";
    char output[15] = "out_test";

    readFile(inputA, A);
    readFile(inputB, B);
    

    /* Start Clock */
    printf("\nStarting clock.\n");
    gettimeofday(&etstart, &tzdummy);
    etstart2 = times(&cputstart);

    // printf("\nStep 1 - Starting 2D FFT for A.\n");
    c_fft2d(A, N, -1, 1);
    c_fft2d(B, N, -1, 0);

    // printf("\nStep 2 - Matrix Mutiplication.\n");
    mm(C, A, B, N);
    print_array(C);
    
    // printf("\nStep 3 - Starting 2D FFT for C.\n");
    c_fft2d(C, N, 1, 1);


    /* Stop Clock */
    gettimeofday(&etstop, &tzdummy);
    etstop2 = times(&cputstop);
    printf("Stopped clock.\n");
    usecstart = (unsigned long long)etstart.tv_sec * 1000000 + etstart.tv_usec;
    usecstop = (unsigned long long)etstop.tv_sec * 1000000 + etstop.tv_usec;


    writeToFile(output, C);

    /* Display timing results */
    printf("\nElapsed time = %g s.\n",
    (float)(usecstop - usecstart)/(float)1000000);

    printf("(CPU times are accurate to the nearest %g ms)\n",
    1.0/(float)CLOCKS_PER_SEC * 1000.0);
    printf("My total CPU time for parent = %g ms.\n",
    (float)( (cputstop.tms_utime + cputstop.tms_stime) -
         (cputstart.tms_utime + cputstart.tms_stime) ) /
    (float)CLOCKS_PER_SEC * 1000);
    printf("My system CPU time for parent = %g ms.\n",
    (float)(cputstop.tms_stime - cputstart.tms_stime) /
    (float)CLOCKS_PER_SEC * 1000);
    printf("My total CPU time for child processes = %g ms.\n",
    (float)( (cputstop.tms_cutime + cputstop.tms_cstime) -
         (cputstart.tms_cutime + cputstart.tms_cstime) ) /
    (float)CLOCKS_PER_SEC * 1000);
        /* Contrary to the man pages, this appears not to include the parent */
    printf("--------------------------------------------\n");

    return 0;
}