Linear algebra refactor

[marcelo-alvarez]

The blas / lapack routines dot, axpy, syr, syrk, gemv, gemm, posv and potri are now called directly using the cblas and lapacke C interfaces, instead of boost. The resulting PSF fit is bitwise identical to that in the master branch when run on cori and dirac, with no changes to desispec or the build commands.

Please check that it builds and runs as expected, and merge into master if successful.

[sbailey]

Thanks, it compiles as expected for my environment at NERSC.

Maintainability question: It looks like the call structure is now (picking a specific example):

specex::Legendre1DPol ->
    specex::syr (specex_linalg) ->
        specex_syr (specex_blas) ->
            cblas_dsyr (from either mkl or cblas)

vs. previously:

specex::Legendre1DPol ->
    specex::syr (specex_linalg) ->
        blas::syr (boost::numeric::bindings::blas via specex_linalg)
            --> (presumably to cblas_dsyr but I haven't tracked into boost)

It looks like specex_syr is taking the place of boost::numeric::bindings::blas::syrso the call stack isn't any deeper than it used to be, but there is now an extra layer to maintain within specex itself. i.e. I was surprised by having both specex::syr and specex_syr. Is that extra internal layer specifically needed / advantageous, or is it a leftover?

Put another way, previously:

algorithm code -> specex::syr interface wrapper -> external library

vs now it looks like

algorithm code -> specex::syr interface wrapper -> another specex_syr interface wrapper -> external library

I understand the benefits of having one interface layer to the external libraries, but this smells like two interface layers. Comments?

[marcelo-alvarez]

I understand the benefits of having one interface layer to the external libraries, but this smells like two interface layers. Comments?

I agree. This was done to have the prototype declarations (through vendor supplied include files) be compiled from separate source to that containing any objects with harp/boost dependencies, to avoid name conflicts, e.g.

mkl_lapack.h: error: conflicting declaration of C function 'void zunmbr'
...
In file included from /...include/boost/numeric/.../blas_names.h:

Once boost and harp dependencies are removed, then it will be easy to remove the additional layer and have the calls to cblas/lapacke directly in specex_linalg.cc. I suggest we go ahead and merge as is and address this in the next PR which will explicitly remove remaining boost and harp dependencies in all remaining files other than for program options.

[sbailey]

OK, thanks for the explanation. That makes sense for why it is structured as it is. OK to merge when ready. Given that your last commit was 5 min ago I won't merge until confirming realtime with you so that I don't merge out from under you.