#37881

Committed 25 Aug 2024 03:35AM UTC coverage: 85.137% (-2.1%) from 87.259%

Build # #37881

Build Type

push

local

Committed by

web-flow

Commit Message

Fix self-recursion in generic triangular (l/r)mul! (#55547)

This fixes a stack-overflow in the following:
```julia
julia> using LinearAlgebra

julia> struct MyTriangularWithoutLRMul{T, A<:LinearAlgebra.AbstractTriangular{T}} <: LinearAlgebra.AbstractTriangular{T}
           data :: A
       end

julia> Base.size(A::MyTriangularWithoutLRMul) = size(A.data)

julia> Base.getindex(A::MyTriangularWithoutLRMul, i::Int, j::Int) = A.data[i,j]

julia> M = MyTriangularWithoutLRMul(UpperTriangular(rand(4,4)));

julia> A = rand(4,4);

julia> lmul!(M, A)
Warning: detected a stack overflow; program state may be corrupted, so further execution might be unreliable.
ERROR: StackOverflowError:
Stacktrace:
     [1] unsafe_copyto!
       @ ./genericmemory.jl:122 [inlined]
     [2] _copyto_impl!
       @ ./array.jl:308 [inlined]
     [3] copyto!
       @ ./array.jl:299 [inlined]
     [4] copyto!
       @ ./array.jl:322 [inlined]
     [5] _trimul!(C::Matrix{Float64}, A::MyTriangularWithoutLRMul{Float64, UpperTriangular{Float64, Matrix{Float64}}}, B::Matrix{Float64})
       @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:961
     [6] lmul!(A::MyTriangularWithoutLRMul{Float64, UpperTriangular{Float64, Matrix{Float64}}}, B::Matrix{Float64})
       @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:982--- the above 2 lines are repeated 39990 more times ---
```
This is done by rerouting the generic `lmul!`/`rmul!` methods to those
for `UpperTriangular` or `LowerTriangular`, depending on which
triangular half is populated.

A similar issue with `ldiv!`/`rdiv!` is also resolved.

Run Details

16 of 16 new or added lines in 1 file covered. (100.0%)

2397 existing lines in 73 files now uncovered.

75732 of 88953 relevant lines covered (85.14%)

14910158.09 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

88.89

/base/simdloop.jl

# This file is a part of Julia. License is MIT: https://julialang.org/license

# Support for @simd for

module SimdLoop

export @simd, simd_outer_range, simd_inner_length, simd_index

# Error thrown from ill-formed uses of @simd
struct SimdError <: Exception
    msg::String
end

# Parse iteration space expression
#       symbol '=' range
#       symbol 'in' range
function parse_iteration_space(x)
    (isa(x, Expr) && (x.head === :(=) || x.head === :in)) || throw(SimdError("= or in expected"))
    length(x.args) == 2 || throw(SimdError("simd range syntax is wrong"))
    isa(x.args[1], Symbol) || throw(SimdError("simd loop index must be a symbol"))
    x.args # symbol, range
end

# reject invalid control flow statements in @simd loop body
function check_body!(x::Expr)
    if x.head === :break || x.head === :continue
        throw(SimdError("$(x.head) is not allowed inside a @simd loop body"))
    elseif x.head === :macrocall && x.args[1] === Symbol("@goto")
        throw(SimdError("@goto is not allowed inside a @simd loop body"))
    end
    for arg in x.args
        check_body!(arg)
    end
    return true
end
check_body!(x::QuoteNode) = check_body!(x.value)
check_body!(x) = true

# @simd splits a for loop into two loops: an outer scalar loop and
# an inner loop marked with :loopinfo. The simd_... functions define
# the splitting.
# Custom iterators that do not support random access cannot support
# vectorization. In order to be compatible with `@simd` annotated loops,
#they should override `simd_inner_length(v::MyIter, j) = 1`,
#`simd_outer_range(v::MyIter) = v`, and `simd_index(v::MyIter, j, i) = j`.

# Get range for outer loop.
simd_outer_range(r) = 0:0

# Get trip count for inner loop.
@inline simd_inner_length(r, j) = Base.length(r)

# Construct user-level element from original range, outer loop index j, and inner loop index i.
@inline simd_index(r, j, i) = (@inbounds ret = r[i+firstindex(r)]; ret)

# Compile Expr x in context of @simd.
function compile(x, ivdep)
    (isa(x, Expr) && x.head === :for) || throw(SimdError("for loop expected"))
    length(x.args) == 2 || throw(SimdError("1D for loop expected"))
    check_body!(x)

    var,range = parse_iteration_space(x.args[1])
    r = gensym("r") # Range value
    j = gensym("i") # Iteration variable for outer loop
    n = gensym("n") # Trip count for inner loop
    i = gensym("i") # Trip index for inner loop
    quote
        # Evaluate range value once, to enhance type and data flow analysis by optimizers.
        let $r = $range
            for $j in Base.simd_outer_range($r)
                let $n = Base.simd_inner_length($r,$j)
                    if zero($n) < $n
                        # Lower loop in way that seems to work best for LLVM 3.3 vectorizer.
                        let $i = zero($n)
                            while $i < $n
                                local $var = Base.simd_index($r,$j,$i)
                                $(x.args[2])        # Body of loop
                                $i += 1
                                $(Expr(:loopinfo, Symbol("julia.simdloop"), ivdep))  # Mark loop as SIMD loop
                            end
                        end
                    end
                end
            end
        end
        nothing
    end
end

"""
    @simd

Annotate a `for` loop to allow the compiler to take extra liberties to allow loop re-ordering

!!! warning
    This feature is experimental and could change or disappear in future versions of Julia.
    Incorrect use of the `@simd` macro may cause unexpected results.

The object iterated over in a `@simd for` loop should be a one-dimensional range.
By using `@simd`, you are asserting several properties of the loop:

* It is safe to execute iterations in arbitrary or overlapping order, with special consideration for reduction variables.
* Floating-point operations on reduction variables can be reordered or contracted, possibly causing different results than without `@simd`.

In many cases, Julia is able to automatically vectorize inner for loops without the use of `@simd`.
Using `@simd` gives the compiler a little extra leeway to make it possible in more situations. In
either case, your inner loop should have the following properties to allow vectorization:

* The loop must be an innermost loop
* The loop body must be straight-line code. Therefore, [`@inbounds`](@ref) is
    currently needed for all array accesses. The compiler can sometimes turn
    short `&&`, `||`, and `?:` expressions into straight-line code if it is safe
    to evaluate all operands unconditionally. Consider using the [`ifelse`](@ref)
    function instead of `?:` in the loop if it is safe to do so.
* Accesses must have a stride pattern and cannot be "gathers" (random-index
    reads) or "scatters" (random-index writes).
* The stride should be unit stride.

!!! note
    The `@simd` does not assert by default that the loop is completely free of loop-carried
    memory dependencies, which is an assumption that can easily be violated in generic code.
    If you are writing non-generic code, you can use `@simd ivdep for ... end` to also assert that:

* There exists no loop-carried memory dependencies
* No iteration ever waits on a previous iteration to make forward progress.
"""
macro simd(forloop)
    esc(compile(forloop, nothing))
end

macro simd(ivdep, forloop)
    if ivdep === :ivdep
        esc(compile(forloop, Symbol("julia.ivdep")))
    else
        throw(SimdError("Only ivdep is valid as the first argument to @simd"))
    end
end

end # module SimdLoop

1	# This file is a part of Julia. License is MIT: https://julialang.org/license
2
3	# Support for @simd for
4
5	module SimdLoop
6
7	export @simd, simd_outer_range, simd_inner_length, simd_index
8
9	# Error thrown from ill-formed uses of @simd
10	struct SimdError <: Exception
UNCOV 11	msg::String	×
12	end
13
14	# Parse iteration space expression
15	# symbol '=' range
16	# symbol 'in' range
17	function parse_iteration_space(x)	3✔
18	(isa(x, Expr) && (x.head === :(=) \|\| x.head === :in)) \|\| throw(SimdError("= or in expected"))	3✔
19	length(x.args) == 2 \|\| throw(SimdError("simd range syntax is wrong"))	3✔
20	isa(x.args[1], Symbol) \|\| throw(SimdError("simd loop index must be a symbol"))	3✔
21	x.args # symbol, range	3✔
22	end
23
24	# reject invalid control flow statements in @simd loop body
25	function check_body!(x::Expr)	15✔
26	if x.head === :break \|\| x.head === :continue	30✔
UNCOV 27	throw(SimdError("$(x.head) is not allowed inside a @simd loop body"))	×
28	elseif x.head === :macrocall && x.args[1] === Symbol("@goto")	15✔
UNCOV 29	throw(SimdError("@goto is not allowed inside a @simd loop body"))	×
30	end
31	for arg in x.args	15✔
32	check_body!(arg)	72✔
33	end	36✔
34	return true	15✔
35	end
36	check_body!(x::QuoteNode) = check_body!(x.value)	×
37	check_body!(x) = true	×
38
39	# @simd splits a for loop into two loops: an outer scalar loop and
40	# an inner loop marked with :loopinfo. The simd_... functions define
41	# the splitting.
42	# Custom iterators that do not support random access cannot support
43	# vectorization. In order to be compatible with `@simd` annotated loops,
44	#they should override `simd_inner_length(v::MyIter, j) = 1`,
45	#`simd_outer_range(v::MyIter) = v`, and `simd_index(v::MyIter, j, i) = j`.
46
47	# Get range for outer loop.
48	simd_outer_range(r) = 0:0	150✔
49
50	# Get trip count for inner loop.
51	@inline simd_inner_length(r, j) = Base.length(r)	487,872✔
52
53	# Construct user-level element from original range, outer loop index j, and inner loop index i.
54	@inline simd_index(r, j, i) = (@inbounds ret = r[i+firstindex(r)]; ret)	111,727,101✔
55
56	# Compile Expr x in context of @simd.
57	function compile(x, ivdep)	3✔
58	(isa(x, Expr) && x.head === :for) \|\| throw(SimdError("for loop expected"))	3✔
59	length(x.args) == 2 \|\| throw(SimdError("1D for loop expected"))	3✔
60	check_body!(x)	3✔
61
62	var,range = parse_iteration_space(x.args[1])	3✔
63	r = gensym("r") # Range value	3✔
64	j = gensym("i") # Iteration variable for outer loop	3✔
65	n = gensym("n") # Trip count for inner loop	3✔
66	i = gensym("i") # Trip index for inner loop	3✔
67	quote	3✔
68	# Evaluate range value once, to enhance type and data flow analysis by optimizers.
69	let $r = $range	5,027,516✔
70	for $j in Base.simd_outer_range($r)	4,967,950✔
71	let $n = Base.simd_inner_length($r,$j)	5,807,368✔
72	if zero($n) < $n	5,854,272✔
73	# Lower loop in way that seems to work best for LLVM 3.3 vectorizer.
74	let $i = zero($n)	5,567,929✔
75	while $i < $n	137,638,018✔
76	local $var = Base.simd_index($r,$j,$i)	131,891,097✔
77	$(x.args[2]) # Body of loop	131,949,400✔
78	$i += 1	131,891,073✔
79	$(Expr(:loopinfo, Symbol("julia.simdloop"), ivdep)) # Mark loop as SIMD loop	131,891,073✔
80	end	131,891,073✔
81	end
82	end
83	end
84	end	5,939,025✔
85	end
86	nothing	4,836,598✔
87	end
88	end
89
90	"""
91	@simd
92
93	Annotate a `for` loop to allow the compiler to take extra liberties to allow loop re-ordering
94
95	!!! warning
96	This feature is experimental and could change or disappear in future versions of Julia.
97	Incorrect use of the `@simd` macro may cause unexpected results.
98
99	The object iterated over in a `@simd for` loop should be a one-dimensional range.
100	By using `@simd`, you are asserting several properties of the loop:
101
102	* It is safe to execute iterations in arbitrary or overlapping order, with special consideration for reduction variables.
103	* Floating-point operations on reduction variables can be reordered or contracted, possibly causing different results than without `@simd`.
104
105	In many cases, Julia is able to automatically vectorize inner for loops without the use of `@simd`.
106	Using `@simd` gives the compiler a little extra leeway to make it possible in more situations. In
107	either case, your inner loop should have the following properties to allow vectorization:
108
109	* The loop must be an innermost loop
110	* The loop body must be straight-line code. Therefore, [`@inbounds`](@ref) is
111	currently needed for all array accesses. The compiler can sometimes turn
112	short `&&`, `\|\|`, and `?:` expressions into straight-line code if it is safe
113	to evaluate all operands unconditionally. Consider using the [`ifelse`](@ref)
114	function instead of `?:` in the loop if it is safe to do so.
115	* Accesses must have a stride pattern and cannot be "gathers" (random-index
116	reads) or "scatters" (random-index writes).
117	* The stride should be unit stride.
118
119	!!! note
120	The `@simd` does not assert by default that the loop is completely free of loop-carried
121	memory dependencies, which is an assumption that can easily be violated in generic code.
122	If you are writing non-generic code, you can use `@simd ivdep for ... end` to also assert that:
123
124	* There exists no loop-carried memory dependencies
125	* No iteration ever waits on a previous iteration to make forward progress.
126	"""
127	macro simd(forloop)	3✔
128	esc(compile(forloop, nothing))	3✔
129	end
130
131	macro simd(ivdep, forloop)
132	if ivdep === :ivdep
133	esc(compile(forloop, Symbol("julia.ivdep")))
134	else
135	throw(SimdError("Only ivdep is valid as the first argument to @simd"))
136	end
137	end
138
139	end # module SimdLoop

JuliaLang / julia / #37881

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous