Optional BLAS/LAPACK implementations

[TLCFEM]

Will it be a good idea to provide an option to offload linear algerabra related operations to third party implementations such as OpenBLAS and MKL?

The linalg module can be mostly offloaded. With the optimised parallel implementations, it may perform better for large problems.

Alternatively, openmp primitives can be used to parallelise some loops I presume.

[zaikunzhang]

Hi @TLCFEM , thank you for raising this!

First of all, the linalg module should never be removed. We should never assume that users have BLAS etc available and (more importantly) are capable of linking to it properly.

Secondly, the cost of numerical linear algebra is usually not an issue for the problems that PRIMA intends to solve. In these problems, costs are dominated by function evaluations, each of which may take minutes, hours, or even months. On the other hand, the most expensive linear algebra in PRIMA may be solving a linear system of 1000 dimensions, because PRIMA is designed to solve problems of at most hundreds of variables (larger problems should be solved by other methods).

Thirdly, it will be of course quite nice to use BLAS as long as this is done in an elegant and user-friendly way, and makes sure that things fall back automatically and gracefully to the internal linalg when BLAS does not exist or is not linked properly.

Thanks.

[zaikunzhang]

Alternatively, openmp primitives can be used to parallelise some loops I presume

As a reference implementation, PRIMA intends to implement linear algebra using matrix-vector procedures instead of loops, whenever the former is possible. As the one who coded each line of PRIMA, I am pretty sure that it does not contain any parallelizable loops. In case they still exist, tell me and I will rewrite them using matrix-vector operations.

There may exist some performance differences between matrix-vector operations and parallelised loops, but this really does not matter (or even make any difference), recalling that each function evaluation takes maybe hours. However, the former is much easier to read, understand, maintain, and extend.

If we prefer loops, then we can just stay with the original F77 implementation, which coded everything based on loops. We know how the code looks like.

Thank you.

[zaikunzhang]

Ref:

https://github.com/jacobwilliams/slsqp/releases/tag/1.5.0

jacobwilliams/slsqp@dffab84