BLAS threads

[filoferra]

Draft, for discussion.

Since PETSc_jll 3.25, all PETSc builds call BLAS through Julia's libblastrampoline, so PETSc's BLAS runs in Julia's own OpenBLAS thread pool. Under MPI, every rank has such a pool, and its idle threads keep spinning. With several ranks on one node, they fight the ranks for the cores. A 4-rank DMStag Stokes solve ran 25× slower (2.0 s instead of 0.07 s per step), with the same results and no warning.

PETSc's -blas_num_threads doesn't help: PETSc can't see that it's OpenBLAS, so it only stores the number.

• Commit 1 (fix): initialize passes -blas_num_threads on to BLAS.set_num_threads. It also stops initialize from changing the caller's options vector. Test: test/blas_threads.jl.
• Commit 2 (to discuss): when several ranks share a node, initialize uses 1 BLAS thread, unless -blas_num_threads, OPENBLAS_NUM_THREADS or OMP_NUM_THREADS is set. Serial runs and one rank per node are unchanged. This changes a global Julia setting, so we can also drop this commit and just document the option. Test: test/mpi_blas_threads.jl, run on 4 ranks.

[filoferra]

Commit 3 (fix, found casually by this PR's new test on Windows): on Windows, no options given to initialize ever reached PETSc, because PETSc there can't see PETSC_OPTIONS once Julia sets it. initialize now passes them as PETSc's command line.

[boriskaus]

yes, in LaMEM.jl I do this and start julia with julia -t 1 to avoid this issue. We cannot have multithreading and MPI at the same time, they bite each other.

[filoferra]

Okay then now the default is one BLAS thread per rank when ranks share a node.
The user can still set -blas_num_threads, OPENBLAS_NUM_THREADS or OMP_NUM_THREADS, or calls BLAS.set_num_threads after initialize if wished otherwise