support non-BLAS numeric types

[stevengj]

It would be good to support non-BLAS numeric types where possible, which may mean writing a simple reference fallback to some of the BLAS routines you are calling.

On the other hand, the Base LinearAlgebra package does not currently support this for the Hessenberg factorization, so it's not a deal-breaker.

[smataigne]

The issue with mul! and Int32 doesn't occure in hessenberg where the type is casted to float, but for the generic method mul! defined for the skewhermitian "class" and tested in the testset "SkewLinearAlgebra.jl". Therefore matrices initialized with Int32 cannot be multiplied by vectors of that type. But this is not a big deal

[stevengj]

Can you be more specific? It seems to work for me:

julia> A = rand(Int32, 5,5);

julia> mul!(similar(A), A, A)
5×5 Matrix{Int32}:
   872558029  -2043782810   -384672508   -933696300  -2049908152
   318733205    499145450    112310675   -184415062    631224005
 -1714681140  -1813236787  -1986934138  -1833312600  -1243817060
   830155170    -24749663  -1104103074    497580760   1956840138
 -2101335698   1987201211   2016869374    306843582  -1894678931

[smataigne]

When I write these commands with A initialized as a skewhermitian matrix of Int32:

x = rand(T, n)
y = zeros(T, n)
mul!(y, A, x, 2, 0)

I obtain the following error message

[stevengj]

That's because the arguments 2 and 0 are of type Int (== Int64 on 64-bit machines), so the arithmetic is promoted to 64-bit integers, but then it can't store the result back in to the 32-bit array y.

If you do mul!(y, A, x, T(2), T(0)) it works, because then everything is of type T = Int32, and there is never any attempt to truncate a 64-bit result to a 32-bit result. (It still overflows, but purely Int32 arithmetic in Julia has 2s-complement overflow behavior.)

[smataigne]

Thank you very much!