Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

rename ExtendedMul functions to a unified unsafe_mul! interface #4

Merged
merged 5 commits into from
Mar 13, 2023

Conversation

Omar-Elrefaei
Copy link
Member

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:
image

@MichelJuillard
Copy link
Member

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
Copy link
Member Author

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:

@Omar-Elrefaei Omar-Elrefaei marked this pull request as ready for review February 24, 2023 03:28
@MichelJuillard MichelJuillard self-requested a review February 24, 2023 07:11
@MichelJuillard
Copy link
Member

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

- Major reorganization
- Add proper QuasiUpperTriangular support in ExtendedMul directy
@Omar-Elrefaei Omar-Elrefaei changed the title rename ExtendedMul functions to a unified _mul! interface rename ExtendedMul functions to a unified unsafe_mul! interface Mar 8, 2023
@Omar-Elrefaei Omar-Elrefaei marked this pull request as draft March 8, 2023 14:30
@Omar-Elrefaei Omar-Elrefaei marked this pull request as ready for review March 9, 2023 20:54
@Omar-Elrefaei
Copy link
Member Author

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 MichelJuillard merged commit af572e8 into main Mar 13, 2023
@MichelJuillard
Copy link
Member

I added the missing test cases for a_mul_b_kron_c!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants