rename ExtendedMul functions to a unified unsafe_mul! interface

[Omar-Elrefaei]

I opened the PR as a draft because one test does not pass yet, the QuasiUpperTriangular one.

Everything else works delightfully.

originating from the kron_mul_elem! call in runtests.jl:273, it errors when it gets from the second (abbreviated) A*B' method, to the full one, to the BLAS Ccall with the following:

ERROR: conversion to pointer not defined for QuasiUpperTriangular{Float64, Matrix{Float64}}
  Stacktrace:
    [1] error(s::String)
      @ Base .\error.jl:33
    [2] unsafe_convert(#unused#::Type{Ptr{Float64}}, a::QuasiUpperTriangular{Float64, Matrix{Float64}})
      @ Base .\pointer.jl:67
    [3] unsafe_convert(#unused#::Type{Ptr{Float64}}, A::Adjoint{Float64, QuasiUpperTriangular{Float64, Matrix{Float64}}})
      @ LinearAlgebra C:\Users\elom\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\adjtrans.jl:197
    [4] pointer
      @ .\abstractarray.jl:1167 [inlined]
    [5] unsafe_convert(P::Type{Ptr{Float64}}, b::Base.RefArray{Float64, Adjoint{Float64, QuasiUpperTriangular{Float64, Matrix{Float64}}}, Nothing})
      @ Base .\refpointer.jl:119
    [6] _mul!(c::Vector{Float64}, offset_c::Int64, a::Vector{Float64}, offset_a::Int64, ma::Int64, na::Int64, 
              b::Adjoint{Float64, QuasiUpperTriangular{Float64, Matrix{Float64}}}, offset_b::Int64, nb::Int64)
      @ KroneckerTools C:\Users\elom\.julia\dev\KroneckerTools\src\ExtendedMul.jl:75

I can confirm from a REPL:

[MichelJuillard]

YES! We need to do similar changes in QuasiUpperTriangular.jl, I had forgotten about it. We can't use the methods in ExtendedMul.jl for QuasiUpperTriangular.jl: we need specialized code. To be consistent with the idea that direct multiplications with sub-arrays are only done in KroneckerProductTools.jl, we need to move and transform all the A_mul_B methods with offset... and so on from QuasiUpperTriangular.jl to ExtendedMul.jl. Then, change the names A_mul_B to mul! for the remaining methods (without offset... and so on) in QuasiUpperTriangular.jl In QuasiUpperTriangular.jl, we also use alpha and beta (that are always 1.0 and 0.0 in ExtendedMul.jl). This makes me think that we should add methods for alpha and beta in ExtendedMul.jl

[Omar-Elrefaei]

We can't use the methods in ExtendedMul.jl for QuasiUpperTriangular.jl: we need specialized code.

Turns out we can! QuasiTriangular just was not implementing the necessary interface to be compatible with BLAS operations. Once I defined the functions in DynareJulia/QuasiTriangular.jl#1, all the tests passed flawlessly without needing any of the A_mul_B functions from QuasiTriangular. Of course we can still port them as appropriate becase they might have a more specialized or efficient algorithm, but now the generic _mul!/Ccall OpenBLAS and LinearAlgebra.BLAS work as they should.

ie. The DynareJulia/QuasiTriangular.jl/#1 PR should make this current PR ready to merge and makes our job with QuasiTriangular much easier.

source: https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-strided-arrays more context:

• https://docs.julialang.org/en/v1/base/arrays/#Base.StridedArray
• https://discourse.julialang.org/t/more-specific-definition-for-densearray/33591/6
• https://discourse.julialang.org/t/more-specific-definition-for-densearray/33591/11

[MichelJuillard]

This is because the tests in test_quasi_upper_triangular.jl only use matrices whose storage (X.data) is quasi upper triangular as they are created by a Schur decomposition. In the tests we should use matrices that have garbage in their lower triangular but not contiguous nonzero on the first sub-diagonal

[Omar-Elrefaei]

I've added tests for a few branches and deleted the corresponding #FIXME: UNTESTED labels. I was not able to quickly construct valid tests for a_mul_b_kron_c! and at_mul_b_kron_c! though.

[MichelJuillard]

I added the missing test cases for a_mul_b_kron_c!