1	/*
     2	simd.c (SIMD)
     3	
     4	VC++ : cl.exe /O2 simd.c
     5	gcc  : gcc -O2 -mavx simd.c -o simd
     6	
     7	Usage:
     8	> simd <1|2> <num> <loop> [sse|avx]
     9	*/
    10	
    11	#include <stdlib.h>
    12	#include <stdio.h>
    13	#include <string.h>
    14	#include <time.h>
    15	#include <immintrin.h>
    16	
    17	static void vadd(int simd, int n, const float a[], const float b[], float c[])
    18	{
    19		int    ne = 0;
    20	
    21		if      (simd == 1) {
    22			// SSE
    23			ne = (n / 4) * 4;
    24			for (int i = 0; i < ne; i += 4) {
    25				_mm_storeu_ps(c + i,
    26					_mm_add_ps(
    27						_mm_loadu_ps(a + i),
    28						_mm_loadu_ps(b + i)));
    29			}
    30		}
    31		else if (simd == 2) {
    32			// AVX
    33			ne = (n / 8) * 8;
    34			for (int i = 0; i < ne; i += 8) {
    35				_mm256_storeu_ps(c + i,
    36					_mm256_add_ps(
    37						_mm256_loadu_ps(a + i),
    38						_mm256_loadu_ps(b + i)));
    39			}
    40		}
    41	
    42		for (int i = ne; i < n; i++) {
    43			c[i] = a[i] + b[i];
    44		}
    45	}
    46	
    47	static float sdot(int simd, int n, const float a[], const float b[])
    48	{
    49		int    ne = 0;
    50		float  sum = 0;
    51	
    52		if      (simd == 1) {
    53			// SSE
    54			float fsum[4];
    55			__m128 vsum;
    56			ne = (n / 4) * 4;
    57			vsum = _mm_setzero_ps();
    58			for (int i = 0; i < ne; i += 4) {
    59				vsum = _mm_add_ps(vsum,
    60				       _mm_mul_ps(_mm_loadu_ps(a + i),
    61				                  _mm_loadu_ps(b + i)));
    62			}
    63			_mm_storeu_ps(fsum, vsum);
    64			sum = fsum[0] + fsum[1] + fsum[2] + fsum[3];
    65		}
    66		else if (simd == 2) {
    67			// AVX
    68			float fsum[8];
    69			__m256 vsum;
    70			ne = (n / 8) * 8;
    71			vsum = _mm256_setzero_ps();
    72			for (int i = 0; i < ne; i += 8) {
    73				vsum = _mm256_add_ps(vsum,
    74				       _mm256_mul_ps(_mm256_loadu_ps(a + i),
    75				                     _mm256_loadu_ps(b + i)));
    76			}
    77			_mm256_storeu_ps(fsum, vsum);
    78			sum = fsum[0] + fsum[1] + fsum[2] + fsum[3]
    79			    + fsum[4] + fsum[5] + fsum[6] + fsum[7];
    80		}
    81	
    82		for (int i = ne; i < n; i++) {
    83			sum += a[i] * b[i];
    84		}
    85	
    86		return sum;
    87	}
    88	
    89	int main(int argc, char **argv)
    90	{
    91		int    calc = 1;
    92		int    simd = 0;
    93		int    n = 1000;
    94		int    nloop = 1000;
    95		float *a = NULL, *b = NULL, *c = NULL;
    96		char *str[3] = {"none", "SSE", "AVX"};
    97	
    98		// arguments
    99		if (argc >= 4) {
   100			calc = atoi(argv[1]);
   101			n = atoi(argv[2]);
   102			nloop = atoi(argv[3]);
   103		}
   104		if (argc >= 5) {
   105			if      (!strcmp(argv[4], "sse")) {
   106				simd = 1;
   107			}
   108			else if (!strcmp(argv[4], "avx")) {
   109				simd = 2;
   110			}
   111		}
   112	
   113		// alloc
   114		size_t size = n * sizeof(float);
   115		if      (calc == 1) {
   116			a = (float *)malloc(size);
   117			b = (float *)malloc(size);
   118			c = (float *)malloc(size);
   119		}
   120		else if (calc == 2) {
   121			a = (float *)malloc(size);
   122			b = (float *)malloc(size);
   123		}
   124	
   125		// setup problem
   126		for (int i = 0; i < n; i++) {
   127			a[i] = i + 1.0f;
   128			b[i] = i + 1.0f;
   129		}
   130	
   131		// timer
   132		clock_t t0 = clock();
   133	
   134		// calculation
   135		double result = 0, exact = 0;
   136		if      (calc == 1) {
   137			// vadd
   138			for (int loop = 0; loop < nloop; loop++) {
   139				vadd(simd, n, a, b, c);
   140			}
   141			result = 0;
   142			for (int i = 0; i < n; i++) {
   143				result += c[i];
   144			}
   145			exact = (double)n * (n + 1);
   146		}
   147		else if (calc == 2) {
   148			// sdot
   149			double sum = 0;
   150			for (int loop = 0; loop < nloop; loop++) {
   151				sum += sdot(simd, n, a, b);
   152			}
   153			result = sum / nloop;
   154			exact = (double)n * (n + 1) * (2 * n + 1) / 6.0;
   155		}
   156	
   157		// timer
   158		clock_t t1 = clock();
   159		double cpu = (double)(t1 - t0) / CLOCKS_PER_SEC;
   160	
   161		// output
   162		printf("%d N=%d L=%d %.6e(%.6e) err=%.1e %.3f[sec] %s\n",
   163			calc, n, nloop, result, exact, fabs(result - exact) / exact, cpu, str[simd]);
   164	
   165		// free
   166		if      (calc == 1) {
   167			free(a);
   168			free(b);
   169			free(c);
   170		}
   171		else if (calc == 2) {
   172			free(a);
   173			free(b);
   174		}
   175	
   176		return 0;
   177	}