erg7.c

#include <stdio.h>
#include <stdlib.h>

#include <sys/time.h>

// this is the plain version (pointer walking) of dot product of two NxN float matrices
// compile with:
// gcc -Wall -O2 matmul-float-simple2.c -o matmul-float-simple2 -DN=1000

// matrix dims N rows x N columns: use -DN=.. to define on compilation



void get_walltime(double *wct) {
  struct timeval tp;
  gettimeofday(&tp,NULL);
  *wct = (double)(tp.tv_sec+tp.tv_usec/1000000.0);
}


int main() {
double ts,te;

float *a,*b,*c;	// matrices A,B,C C=AxB, B is transposed
__m128 *pa,*pb;

int k;

k = posix_memalign((void **)&a,16,N*N*sizeof(float));
if (k!=0) { printf("alloc error\n"); exit(1); }
k = posix_memalign((void **)&b,16,N*N*sizeof(float));
if (k!=0) { free(a); printf("alloc error\n"); exit(1); }
k = posix_memalign((void **)&c,16,N*N*sizeof(float));
if (k!=0) { free(a); free(b); printf("alloc error\n"); exit(1); }

  // init input and output matrices
  for (int i=0;i<N*N;i++) {
      a[i] = 2.0;
      b[i] = 3.0;
      c[i] = 0.0;
  }

  // get starting time (double, seconds)
  get_walltime(&ts);

  // load, matrix multiplication
  for (int i=0;i<N;i++) {	// for all rows of A,C

    pb = (__m128 *)b;

    for (int j=0;j<N;j++) {	// for all "columns" (rows) of B

      pa = (__m128 *)a+N*i;

      __m128 sum = _mm_setzero_ps();

      for (int k=0; k<N; k+=4) {	// for each element of selected A row and B "column"
        sum = _mm_add_ps(sum, _mm_mul_ps(*pa, *pb));
        pa++; pb++;
      }
      _mm_store_ps(&c[i*N+j], sum);

    }

  }

  // get ending time
  get_walltime(&te);

  // print computation time
  printf("Computation time = %f sec\n",(te-ts));

  // test result
  for (int i=0;i<N*N;i++) {
      if (c[i]!=6.0*N) { printf("Error! %d (%f)\n",i,c[i]); break; }
  }

  free(c);
  free(b);
  free(a);

  return 0;
}