numeric-linalg

 #include <stdio.h>
 #include <stdlib.h>
 #include <pthread.h>
 
 #include <stdint.h>
 #include <stdbool.h>
 #include <time.h>
 
 #include <string.h>
 #include <errno.h>
 
 #include "progress-bar.h"
 #include "config.h"
 
 #define HISTOGRAM_SIZE              (MAX_N/STEP)
 #define CLOCKS_PER_MILLIS           (CLOCKS_PER_SEC/1000)
 #define LOGISTICS_INITIAL_CONDITION (-800.)
 
 uint32_t histogram[HISTOGRAM_SIZE];
 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  *data;
   double  *ref_data;
   int32_t ipiv[MAX_N];
   size_t  id;
 } 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);
 extern void dgetrfnaive_(int *m, int *n, double *A, int *lda, int *ipiv, int *info);
 
 uint32_t thread_run_benchmark(Thread *thread, size_t n)
 {
   int64_t total_time = 0;
   clock_t start, end;
   int l_n, l_m, lda, info;
 
   for (size_t i = 0; i < N_TESTS; i++) {
     // 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);
     l_n = n; l_m = n; lda = n;
 
     start = clock();
     getrf(&l_m, &l_n, thread->data, &lda, thread->ipiv, &info);
     end = clock();
     total_time += (end - start)/CLOCKS_PER_MILLIS;
   }
 
   return (uint32_t)(total_time/N_TESTS);
 }
 
 void *thread_benchmark(void *arg)
 {
   Thread *thread = (Thread*) arg;
 
   // Computations are distributed evenly across threads
   for (size_t n = 1 + thread->id*STEP; n < MAX_N; n += STEP*N_THREADS) {
     uint32_t duration = thread_run_benchmark(thread, n);
 
     histogram[(n - 1)/STEP] = duration;
     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)
 {
   const char *output_path;
 
   if (argc < 2 || strcmp(argv[1], "standard") == 0) {
     getrf = dgetrf_;
     output_path = OUTPUT_DIR "histogram-dgetrf.bin";
   } else if (strcmp(argv[1], "naive") == 0) {
     getrf = dgetrfnaive_;
     output_path = OUTPUT_DIR "histogram-naive.bin";
   } else if (strcmp(argv[1], "unblocked") == 0) {
     getrf = dgetrf2_;
     output_path = OUTPUT_DIR "histogram-unblocked.bin";
   } 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 dgetrf function 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 histogram to disk...\n");
 
   FILE *output = fopen(output_path, "w");
 
   if (output == NULL) {
     fprintf(stderr,
             "ERROR: Coundn't open output file \"%s\": %s",
             output_path, strerror(errno));
     return EXIT_FAILURE;
   }
   size_t written = fwrite(histogram, sizeof(uint32_t), HISTOGRAM_SIZE, output);
 
   if (written < HISTOGRAM_SIZE) {
     fprintf(stderr,
             "ERROR: Failed to write histogram to output file \"%s\": %s\n",
             output_path, strerror(errno));
     return EXIT_FAILURE;
   }
 
   printf("INFO: Done!\n");
 
   return EXIT_SUCCESS;
 }