Skip to content

Aliased output matrix in LinearAlgebra.mul! does not throw error for matrix wrapper types #925

@DNF2

Description

@DNF2

This came up here: https://www.reddit.com/r/Julia/comments/ugm4l7/working_with_linearalgebramul/

It looks like a bug to me, but I'm unable to test on 1.8. Perhaps someone can verify?

According to the docstring of LinearAlgebra.mul!:

Note that C must not be aliased with either A or B.

However:

julia> A = rand(5,5);

julia> B_ = rand(5,5);

julia> B = B_ + B_';

julia> C = similar(A); 

julia> mul!(C, A, Symmetric(B))  # works as expected
5×5 Matrix{Float64}:
 1.10149  1.93057  1.67467  1.39908  1.59191
 1.37769  1.07008  2.50417  1.21931  2.02994
 1.82134  2.73737  3.66532  1.5658   3.55056
 1.45685  1.44589  2.73083  1.17807  2.39641
 1.79866  2.58796  2.48849  1.33673  2.46281

julia> mul!(A, A, B)  # fails as expected
ERROR: ArgumentError: output matrix must not be aliased with input matrix
Stacktrace:
 [1] gemm_wrapper!(C::Matrix{Float64}, tA::Char, tB::Char, A::Matrix{Float64}, B::Matrix{Float64}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
   @ LinearAlgebra C:\Users\dnf\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\LinearAlgebra\src\matmul.jl:647
 [2] mul!
   @ C:\Users\dnf\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\LinearAlgebra\src\matmul.jl:169 [inlined]
 [3] mul!(C::Matrix{Float64}, A::Matrix{Float64}, B::Matrix{Float64})
   @ LinearAlgebra C:\Users\dnf\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\LinearAlgebra\src\matmul.jl:275
 [4] top-level scope
   @ REPL[91]:1

julia> mul!(A, A, Symmetric(B))  # does not fail as it should, and produces wrong result
5×5 Matrix{Float64}:
 0.201062  0.334617  0.478842  0.102553  1.21272
 0.32496   0.540814  0.773914  0.165748  1.96002
 0.575235  0.957335  1.36996   0.293402  3.46957
 0.362894  0.603946  0.864256  0.185096  2.18882
 0.439629  0.731652  1.04701   0.224235  2.65165

The same thing happens for other matrix wrapper types, like Diagonal etc.

julia> versioninfo()
Julia Version 1.7.2
Commit bf53498635 (2022-02-06 15:21 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-9750H CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = code

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions