numeric-linalg

 #include <stdio.h>
 #include <stdlib.h>
 
 #define __USE_GNU
 #include <sched.h>
 #include <pthread.h>
 
 #include <stdint.h>
 #include <stdbool.h>
 
 #include <string.h>
 #include <errno.h>
 
 #include "progress-bar.h"
 #include "perf.h"
 #include "config.h"
 
 #define HISTOGRAM_SIZE              (MAX_N/STEP)
 #define LOGISTICS_INITIAL_CONDITION (-800.)
 
 // X-macros for defining variables for each perf event
 #define EVENTS  EVENT(loads) EVENT(misses)
 
 #define EVENT(name) uint64_t histogram_##name[HISTOGRAM_SIZE];
 EVENTS
 #undef  EVENT
 
 void (*getrf)(int *m, int *n, double *A, int *lda, int *ipiv, int *info);
 
 ProgressBar progress = {
   .total = HISTOGRAM_SIZE,
   .count = 0,
   .mutex = PTHREAD_MUTEX_INITIALIZER,
 };
 
 // .data is a pointer because it should be allocated in the heap
 // (.data DOES NOT fit in the stack 🤡)
 typedef struct {
   double       *ref_data;
 
   // Input parameters for dgetrf
   double       *data;
   int32_t      ipiv[MAX_N];
   int          n, m, lda, info;
 
   size_t       id;
   PerfRecorder recorder;
 } Thread;
 
 typedef struct {
   Thread threads[N_THREADS];
 } Benchmarker;
 
 /*
  * LAPACK functions
  */
 extern void dgetrf_(int *m, int *n, double *A, int *lda, int *ipiv, int *info);
 extern void dgetrf2_(int *m, int *n, double *A, int *lda, int *ipiv, int *info);
 
 PerfResult thread_run_benchmark(Thread *thread, size_t n)
 {
   // Reinitializes the values in the .data array: this avoids progressively
   // moving larger values to the beginning of the array, which would decrise
   // decrease the number of row interchanges required for computations
   memcpy(thread->data, thread->ref_data, sizeof(double)*n*n);
   thread->n = n; thread->m = n; thread->lda = n;
 
   perf_start_recording(&thread->recorder, thread->id);
   getrf(&thread->n, &thread->m, thread->data, &thread->lda,
         thread->ipiv, &thread->info);
   perf_stop_recording(&thread->recorder);
 
   return thread->recorder.result;
 }
 
 void *thread_benchmark(void *arg)
 {
   Thread *thread = (Thread*) arg;
 
  // We need to lock the running thread to some CPU so that we can mask
  // perf_event_open with this specific CPU
   cpu_set_t set;
   CPU_SET(thread->id, &set);
   if (sched_setaffinity(0, sizeof(set), &set) == -1) {
     fprintf(stderr, "ERROR: Could not lock CPU affinity for thread %lu\n",
             thread->id);
   }
 
   // Computations are distributed evenly across threads
   for (size_t n = 1 + thread->id*STEP; n < MAX_N; n += STEP*N_THREADS) {
     PerfResult result = thread_run_benchmark(thread, n);
     size_t i = (n - 1)/STEP;
 
     histogram_loads[i]  = result.cache_loads;
     histogram_misses[i] = result.cache_misses;
 
     progress_bar_inc(&progress);
   }
 
   return NULL;
 }
 
 Benchmarker benchmarker_new(double *ref_data)
 {
   Benchmarker bench = {0};
 
   // This array will live for the entire duration of the program,
   // so we might as well leak it 🤡
   double *data = malloc(sizeof(double)*MAX_N*MAX_N*N_THREADS);
   if (data == NULL) {
     fprintf(stderr, "ERROR: Buy more RAM!\n");
     exit(EXIT_FAILURE);
   }
 
   for (size_t i = 0; i < N_THREADS; i++) {
     bench.threads[i].data     = data + i*MAX_N*MAX_N;
     bench.threads[i].ref_data = ref_data;
     bench.threads[i].id       = i;
   }
 
   return bench;
 }
 
 void benchmarker_run(Benchmarker *bench)
 {
   pthread_t handles[N_THREADS] = {0};
   for (size_t i = 0; i < N_THREADS; i++) {
     if (pthread_create(&handles[i], NULL,
                        thread_benchmark,
                        &bench->threads[i]) != 0) {
       fprintf(stderr, "ERROR: Failed to spawn thread %lu!\n", i);
       exit(EXIT_FAILURE);
     }
   }
 
   for (size_t i = 0; i < N_THREADS; i++) {
     pthread_join(handles[i], NULL);
   }
 }
 
 
 int main(int argc, char **argv)
 {
 #define EVENT(name) const char *output_path_##name;
   EVENTS
 #undef  EVENT
 
   if (argc < 2 || strcmp(argv[1], "standard") == 0) {
     getrf = dgetrf_;
 #define EVENT(name) \
     output_path_##name = OUTPUT_DIR "cache-" #name ".bin";
     EVENTS
 #undef EVENT
   } else if (strcmp(argv[1], "unblocked") == 0) {
     getrf = dgetrf2_;
 #define EVENT(name) \
     output_path_##name = OUTPUT_DIR "unblocked-cache-" #name ".bin";
     EVENTS
 #undef EVENT
   } else {
     fprintf(stderr, "ERROR: unknown command \"%s\"\n", argv[1]);
     fprintf(stderr, "USAGE: %s [standard|naive|unblocked]\n", argv[0]);
     exit(EXIT_FAILURE);
   }
 
   // ========================================================================
   printf("INFO: Initializing random input data... ");
 
   // This array will live for the entire duration of the program,
   // so we might as well leak it 🤡
   double *ref_data = malloc(sizeof(double)*MAX_N*MAX_N);
   if (ref_data == NULL) {
     fprintf(stderr, "ERROR: Buy more RAM!\n");
     exit(EXIT_FAILURE);
   }
 
   // Pseudorandom data given by the logistic map
   double acc = LOGISTICS_INITIAL_CONDITION;
   for (size_t i = 0; i < MAX_N*MAX_N; i++) {
     ref_data[i] = acc;
     acc = 1000. - acc*acc/500.;
   }
 
   printf("done!\n");
 
   // ========================================================================
   printf("INFO: Benchmarking the cache miss ratio of dgetrf on %ux%u matrices... (using %u threads)\n",
          MAX_N, MAX_N, N_THREADS);
 
   Benchmarker bench = benchmarker_new(ref_data);
   benchmarker_run(&bench);
 
   // ========================================================================
   printf("INFO: Done with benchmarking! Saving histograms to disk...\n");
 
 #define EVENT(name) \
   do { \
     FILE *output = fopen(output_path_##name, "w"); \
     if (output == NULL) { \
       fprintf(stderr, \
               "ERROR: Coundn't open output file \"%s\": %s", \
               output_path_##name, strerror(errno)); \
       return EXIT_FAILURE; \
     } \
     size_t written = fwrite(&histogram_##name, sizeof(uint64_t), HISTOGRAM_SIZE, output); \
     if (written < HISTOGRAM_SIZE) { \
       fprintf(stderr, \
               "ERROR: Failed to write histogram to output file \"%s\": %s\n", \
               output_path_##name, strerror(errno)); \
       return EXIT_FAILURE; \
     } \
   } while(0);
   EVENTS
 #undef EVENT
 
   printf("INFO: Done!\n");
 
   return EXIT_SUCCESS;
 }