Rigorous benchmark of matrix utilities for parallelism

[JamesYang007]

Most matrix utility functions are in util.hpp under matrix/. We just spotted a huge cost in omp even when n_threads=1! We should do some benchmarks to figure out how to perform the parallelism properly.

[JamesYang007]

• [x] First, rewrite all parallelizable linear algebra routines in terms of the utility functions in utils.hpp. API should be like xxx(..., n_threads).
• [x] We should add global configurations that determine cutoffs for when parallelization is used. Something like "if n_threads <= 1 || n < config.value, then single-thread, otherwise parallelize".
• [x] Next, we need to benchmark these utility functions to see exactly when parallelization is helpful, i.e. what are good cutoffs? There must be a sweet-spot when parallelization is helpful. This is probably machine-dependent? Anyways, worth experimenting on mac m1 and Sherlock.

[JamesYang007]

Found a critical flag of omp! On GOMP, GOMP_SPINCOUNT and on mac brew's libomp KMP_BLOCKTIME should be set sufficiently high that threads don't sleep too soon. OMP is way too costly when waking up the threads.

[JamesYang007]

OMP_PROC_BIND also helps localize memory since each thread is bound to a single CPU.