GPU friendly BLAS calls

[shashi]

@ChrisRackauckas found out that using CuArrays fails at points like https://github.com/JuliaMath/IterativeSolvers.jl/blob/17ef2613dff1d80724886797160ee01bd1bd2c71/src/simple.jl#L38

where the call to BLAS.axpy! explicitly assumes that the inputs are OK with OpenBLAS. The right alternative would be to use LinearAlgebra.axpy! (which as it turns out is a different function) and let it naturally dispatch to BLAS iff possible.

[pkofod]

see also #245 (I actually believe this is a duplicate of that)

[simonbyrne]

There are also several hard-coded Matrix types, e.g.
https://github.com/JuliaMath/IterativeSolvers.jl/blob/17ef2613dff1d80724886797160ee01bd1bd2c71/src/gmres.jl#L7

[pkofod]

Since I need this, I might be able to go over some of these. If anyone has a usecase/ something that doens't currently work, please post it so I can take it into account.

[simonbyrne]

We have some of our own code implemented:
https://github.com/climate-machine/CLIMA/tree/master/src/LinearSolvers
feel free to take/adapt what you need.

See also this comment for our needs: CliMA/ClimateMachine.jl#330 (comment)