Segfault with longer SIMD vectors

[fredrikekre]

The following:

using SIMD
A = Vec{N, Float64}((rand(N)...))
B = A + A

works for N = (1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 16, 17, 18, 20, 24, 32, 33, 34, 36, 40, 48, 64, 65, 66, 68, 72, 80, 96) (tested up to N = 100) on v0.5. On master the above code segfaults for all N above 32 (e.g. N = 33, 34, 36, 40, 48, 64, 65, 66, 68, 72, 80, 96). The generated code look okay(?). Example with N = 36:

julia> @code_llvm A + A

define void @"julia_+_68162"(%Vec* noalias nocapture sret, %Vec*, %Vec*) #0 !dbg !5 {
top:
  %3 = getelementptr inbounds %Vec, %Vec* %1, i64 0, i32 0
  %4 = load <36 x double>, <36 x double>* %3, align 512
  %5 = getelementptr inbounds %Vec, %Vec* %2, i64 0, i32 0
  %6 = load <36 x double>, <36 x double>* %5, align 512
  %res.i = fadd <36 x double> %4, %6
  %.sroa.0.0..sroa_idx = getelementptr inbounds %Vec, %Vec* %0, i64 0, i32 0
  store <36 x double> %res.i, <36 x double>* %.sroa.0.0..sroa_idx, align 512
  ret void
}

julia> @code_native A + A
    .text
Filename: SIMD.jl
    pushq   %rbp
    movq    %rsp, %rbp
Source line: 567
    vmovapd 256(%rsi), %ymm0
    vmovapd 224(%rsi), %ymm1
    vmovapd 192(%rsi), %ymm2
    vmovapd 160(%rsi), %ymm3
    vmovapd 128(%rsi), %ymm4
    vmovapd (%rsi), %ymm5
    vmovapd 32(%rsi), %ymm6
    vmovapd 64(%rsi), %ymm7
    vmovapd 96(%rsi), %ymm8
    vaddpd  (%rdx), %ymm5, %ymm5
    vaddpd  32(%rdx), %ymm6, %ymm6
    vaddpd  64(%rdx), %ymm7, %ymm7
    vaddpd  96(%rdx), %ymm8, %ymm8
    vaddpd  128(%rdx), %ymm4, %ymm4
    vaddpd  160(%rdx), %ymm3, %ymm3
    vaddpd  192(%rdx), %ymm2, %ymm2
    vaddpd  224(%rdx), %ymm1, %ymm1
    vaddpd  256(%rdx), %ymm0, %ymm0
Source line: 993
    vmovapd %ymm0, 256(%rdi)
    vmovapd %ymm1, 224(%rdi)
    vmovapd %ymm2, 192(%rdi)
    vmovapd %ymm3, 160(%rdi)
    vmovapd %ymm4, 128(%rdi)
    vmovapd %ymm8, 96(%rdi)
    vmovapd %ymm7, 64(%rdi)
    vmovapd %ymm6, 32(%rdi)
    vmovapd %ymm5, (%rdi)
    movq    %rdi, %rax
    popq    %rbp
    vzeroupper
    retq
    nopl    (%rax)

julia> A + A

signal (11): Segmentation fault
while loading no file, in expression starting on line 0
+ at /home/fredrik/.julia/v0.6/SIMD/src/SIMD.jl:0
unknown function (ip: 0x7f17757ef5e1)
jl_call_fptr_internal at /home/fredrik/julia-master/src/julia_internal.h:326 [inlined]
jl_call_method_internal at /home/fredrik/julia-master/src/julia_internal.h:345 [inlined]
jl_apply_generic at /home/fredrik/julia-master/src/gf.c:2225
do_call at /home/fredrik/julia-master/src/interpreter.c:75
eval at /home/fredrik/julia-master/src/interpreter.c:230
jl_interpret_toplevel_expr at /home/fredrik/julia-master/src/interpreter.c:34
jl_toplevel_eval_flex at /home/fredrik/julia-master/src/toplevel.c:577
jl_toplevel_eval_in at /home/fredrik/julia-master/src/builtins.c:484
eval at ./boot.jl:235
unknown function (ip: 0x7f1986e6a9bf)
jl_call_fptr_internal at /home/fredrik/julia-master/src/julia_internal.h:326 [inlined]
jl_call_method_internal at /home/fredrik/julia-master/src/julia_internal.h:345 [inlined]
jl_apply_generic at /home/fredrik/julia-master/src/gf.c:2225
eval_user_input at ./REPL.jl:66
unknown function (ip: 0x7f17757e30c6)
jl_call_fptr_internal at /home/fredrik/julia-master/src/julia_internal.h:326 [inlined]
jl_call_method_internal at /home/fredrik/julia-master/src/julia_internal.h:345 [inlined]
jl_apply_generic at /home/fredrik/julia-master/src/gf.c:2225
macro expansion at ./REPL.jl:97 [inlined]
#1 at ./event.jl:73
unknown function (ip: 0x7f17757dad2f)
jl_call_fptr_internal at /home/fredrik/julia-master/src/julia_internal.h:326 [inlined]
jl_call_method_internal at /home/fredrik/julia-master/src/julia_internal.h:345 [inlined]
jl_apply_generic at /home/fredrik/julia-master/src/gf.c:2225
jl_apply at /home/fredrik/julia-master/src/julia.h:1410 [inlined]
start_task at /home/fredrik/julia-master/src/task.c:261
unknown function (ip: 0xffffffffffffffff)
Allocations: 3349005 (Pool: 3347456; Big: 1549); GC: 4
Segmentation fault (core dumped)

Actually, doing anything with the Vec, like showing, also segfaults:

julia> A = Vec{36, Float64}((rand(36)...))
36-element SIMD.Vec{36,Float64}:

signal (11): Segmentation fault
while loading no file, in expression starting on line 0
#showarray#256 at ./show.jl:1697
unknown function (ip: 0x7ff0ba4ff36d)
[...]

[KristofferC]

On 0.5 we say that the vectors should be 8 byte aligned:

 %4 = load <36 x double>, <36 x double>* %3, align 8

On 0.6, vectors above length 32 are said to be 512 byte aligned (?):

32:

%4 = load <32 x double>, <32 x double>* %3, align 8

36:

%4 = load <36 x double>, <36 x double>* %3, align 512

[KristofferC]

I also thought data should be 16 byte aligned for SIMD (32 for AVX).

[RalphAS]

Could this be an LLVM codegen problem? FWIW I get a segfault in the constructor, as one might expect from the following (note the large offsets):

julia> x=rand(36);
julia> @code_native(Vec{36,Float64}((x...)))
	.text
Filename: SIMD.jl
	pushq	%rbp
	movq	%rsp, %rbp
	andq	$-32, %rsp
Source line: 121
	subq	$256, %rsp              # imm = 0x100
	vmovups	(%rdx), %ymm0
	vmovups	32(%rdx), %ymm1
	vmovups	64(%rdx), %ymm2
	vmovups	96(%rdx), %ymm3
	vmovups	128(%rdx), %ymm4
	vmovups	160(%rdx), %ymm5
	vmovups	192(%rdx), %ymm6
	vmovups	224(%rdx), %ymm7
	vmovups	256(%rdx), %ymm8
Source line: 119
	vmovaps	%ymm8, 256(%rdi)
	vmovaps	%ymm7, 224(%rdi)
	vmovaps	%ymm6, 192(%rdi)
	vmovaps	%ymm5, 160(%rdi)
	vmovaps	%ymm4, 128(%rdi)
	vmovaps	%ymm3, 96(%rdi)
	vmovaps	%ymm2, 64(%rdi)
	vmovaps	%ymm1, 32(%rdi)
	vmovaps	%ymm0, (%rdi)
	vmovaps	192(%rsp), %ymm0
	vmovaps	%ymm0, 480(%rdi)
	vmovaps	160(%rsp), %ymm0
	vmovaps	%ymm0, 448(%rdi)
	vmovaps	128(%rsp), %ymm0
	vmovaps	%ymm0, 416(%rdi)
	vmovaps	(%rsp), %ymm0
	vmovaps	32(%rsp), %ymm1
	vmovaps	64(%rsp), %ymm2
	vmovaps	96(%rsp), %ymm3
	vmovaps	%ymm3, 384(%rdi)
	vmovaps	%ymm2, 352(%rdi)
	vmovaps	%ymm1, 320(%rdi)
	vmovaps	%ymm0, 288(%rdi)
	movq	%rdi, %rax
	movq	%rbp, %rsp
	popq	%rbp
	vzeroupper
	retq
	nopw	(%rax,%rax)

The N=32 version looks fine (and seems to work). Output from code_llvm appears benign for both. I use a system image for Haswell, i.e. AVX2.

[yuyichao]

No, the issue is already identified AFAICT. Assembly code isn't very helpful here.

[m-j-w]

I stumbled over a similar issue, thus took a look at type sizes and alignment. The following table shows the sizes and alignments for tuples of N VecElements, with N matching the row number, for integers of size 8 to 128 bits. (the table extends to the right). Ran it on today's Julia master.

Note, that begining with N>32, some base types lead to a Tuple alignment of 0. Also note, most fall back to the alignment and size of the pure NTuple{N,T}, rather than NTuple{N,VecElement{T}}.
As far as I know, loading SIMD vectors with an alignment of <16 bytes is an illegal instruction.
Also, the Int128 is as far as I know not supported at all by SIMD instructions. So giving a warning there might be a solution.

My thought on that would be that the alignment should be the tuple size rounded up to the next power of two, but limited to the cache line size of 64 bytes, and a minimum of the necessary 16 bytes.
And the size of the tuple should probably also be rounded up to the next valid SIMD vector size, whatever that is on the executing machine.
Also, Intel occasionally recommends to always to use 64 byte, see e.g. (https://software.intel.com/en-us/articles/data-alignment-to-assist-vectorization).

For the Intel® Many Integrated Core Architecture (Intel® MIC Architecture) such as the Intel® Xeon Phi™ Coprocessor, memory movement is optimal when the data starting address lies on 64 byte boundaries. Thus, it is desired to force the compiler to create data objects with starting addresses that are modulo 64 bytes.

However, there also seems to be an issue in LLVM codegen itself. I found this possibly related LLVM bug.
I saw a similar error when trying to use a vector size of 7 with SIMD.jl, yielding an llvmcall error.

Unfortunately, I have no clue where to change the alignment in the depth of the Julia codebase. Hope this helps anyway.

Row	Int8_sizeof	Int16_sizeof	Int32_sizeof	Int64_sizeof	Int128_sizeof	Int8_alignm	Int16_alignm	Int32_alignm	Int64_alignm	Int128_alignm
1	1	2	4	8	16	1	2	4	8	8
2	2	4	8	16	32	2	4	8	16	8
3	4	8	16	32	48	4	8	16	32	8
4	4	8	16	32	64	4	8	16	32	8
5	8	16	32	64	80	8	16	32	64	8
6	8	16	32	64	96	8	16	32	64	8
7	7	14	28	56	112	1	2	4	8	8
8	8	16	32	64	128	8	16	32	64	8
9	16	32	64	128	144	16	32	64	128	8
10	16	32	64	128	160	16	32	64	128	8
11	11	22	44	88	176	1	2	4	8	8
12	16	32	64	128	192	16	32	64	128	8
13	13	26	52	104	208	1	2	4	8	8
14	14	28	56	112	224	1	2	4	8	8
15	15	30	60	120	240	1	2	4	8	8
16	16	32	64	128	256	16	32	64	128	8
17	32	64	128	256	272	32	64	128	256	8
18	32	64	128	256	288	32	64	128	256	8
19	19	38	76	152	304	1	2	4	8	8
20	32	64	128	256	320	32	64	128	256	8
21	21	42	84	168	336	1	2	4	8	8
22	22	44	88	176	352	1	2	4	8	8
23	23	46	92	184	368	1	2	4	8	8
24	32	64	128	256	384	32	64	128	256	8
25	25	50	100	200	400	1	2	4	8	8
26	26	52	104	208	416	1	2	4	8	8
27	27	54	108	216	432	1	2	4	8	8
28	28	56	112	224	448	1	2	4	8	8
29	29	58	116	232	464	1	2	4	8	8
30	30	60	120	240	480	1	2	4	8	8
31	31	62	124	248	496	1	2	4	8	8
32	32	64	128	256	512	32	64	128	256	8
33	64	128	256	512	528	64	128	256	0	8
34	64	128	256	512	544	64	128	256	0	8
35	35	70	140	280	560	1	2	4	8	8
36	64	128	256	512	576	64	128	256	0	8
37	37	74	148	296	592	1	2	4	8	8
38	38	76	152	304	608	1	2	4	8	8
39	39	78	156	312	624	1	2	4	8	8
40	64	128	256	512	640	64	128	256	0	8
41	41	82	164	328	656	1	2	4	8	8
42	42	84	168	336	672	1	2	4	8	8
43	43	86	172	344	688	1	2	4	8	8
44	44	88	176	352	704	1	2	4	8	8
45	45	90	180	360	720	1	2	4	8	8
46	46	92	184	368	736	1	2	4	8	8
47	47	94	188	376	752	1	2	4	8	8
48	64	128	256	512	768	64	128	256	0	8
49	49	98	196	392	784	1	2	4	8	8
50	50	100	200	400	800	1	2	4	8	8
51	51	102	204	408	816	1	2	4	8	8
52	52	104	208	416	832	1	2	4	8	8
53	53	106	212	424	848	1	2	4	8	8
54	54	108	216	432	864	1	2	4	8	8
55	55	110	220	440	880	1	2	4	8	8
56	56	112	224	448	896	1	2	4	8	8
57	57	114	228	456	912	1	2	4	8	8
58	58	116	232	464	928	1	2	4	8	8
59	59	118	236	472	944	1	2	4	8	8
60	60	120	240	480	960	1	2	4	8	8
61	61	122	244	488	976	1	2	4	8	8
62	62	124	248	496	992	1	2	4	8	8
63	63	126	252	504	1008	1	2	4	8	8
64	64	128	256	512	1024	64	128	256	0	8
65	128	256	512	1024	1040	128	256	0	0	8
66	128	256	512	1024	1056	128	256	0	0	8
67	67	134	268	536	1072	1	2	4	8	8
68	128	256	512	1024	1088	128	256	0	0	8
69	69	138	276	552	1104	1	2	4	8	8
70	70	140	280	560	1120	1	2	4	8	8
71	71	142	284	568	1136	1	2	4	8	8
72	128	256	512	1024	1152	128	256	0	0	8
73	73	146	292	584	1168	1	2	4	8	8
74	74	148	296	592	1184	1	2	4	8	8
75	75	150	300	600	1200	1	2	4	8	8
76	76	152	304	608	1216	1	2	4	8	8
77	77	154	308	616	1232	1	2	4	8	8
78	78	156	312	624	1248	1	2	4	8	8
79	79	158	316	632	1264	1	2	4	8	8
80	128	256	512	1024	1280	128	256	0	0	8
81	81	162	324	648	1296	1	2	4	8	8

Here's the code to reproduce the above. (Not really a piece of art):

using DataFrames
df = DataFrame(n=1:81)
ttype(t,a) = NTuple{a,VecElement{t}}
label(t,s) = Symbol(string(t),'_',s)
for t in [Int8, Int16, Int32, Int64, Int128]
   types = ttype.(t,df[:n])
   df[label(t,"sizeof")] = sizeof.(types)
end
for t in [Int8, Int16, Int32, Int64, Int128]
   types = ttype.(t,df[:n])
   df[label(t,"alignm")] = Base.datatype_alignment.(types)
end
delete!(df, :n)
showall(df)

[KristofferC]

This is where it is decided if the tuple should be passed as a LLVM Vector or Array:

julia/src/datatype.c

Line 174 in 8341444

// LLVM 3.7 and 3.8 either crash or generate wrong code for many

. The required alignment is computed below. The comment:

// LLVM's alignment rule for vectors seems to be to round up to
// a power of two, even if that's overkill for the target hardware.

seems a bit iffy to me.

Also, note that as long as the SIMD vector does not hit global scope, things seem to work ok:

using SIMD

@generated function tosimd{N, T}(a::NTuple{N, T})
    t = Expr(:tuple, [:(a[$i]) for i in 1:N]...)
    ex = :(SIMD.Vec{N, T}($t))
    return quote
        $(Expr(:meta, :inline))
        $ex
    end
end

function sum_simd{N,T}(t::NTuple{N, T})
    v = tosimd(t)
    return sum(v)
end

sum_simd((rand(36)...))

evaluated cleanly.

[fredrikekre]

It wont work with N other than the ones I listed up top because of:

julia/src/datatype.c

Lines 174 to 185 in 64409a0

	// LLVM 3.7 and 3.8 either crash or generate wrong code for many
	// SIMD vector sizes N. It seems the rule is that N can have at
	// most 2 non-zero bits. (This is true at least for N<=100.) See
	// also <https://llvm.org/bugs/show_bug.cgi?id=27708>.
	size_t mask = nfields;
	// See e.g.
	// <https://graphics.stanford.edu/%7Eseander/bithacks.html> for an
	// explanation of this bit-counting algorithm.
	mask &= mask-1; // clear least-significant 1 if present
	mask &= mask-1; // clear another 1
	if (mask)
	return 0; // nfields has more than two 1s

Edit: You are always ahead of me @KristofferC ...

[m-j-w]

Regardless of the comment you linked refering to LLVM 3.7 & 3.8 crashing, returning the wrong alignment in

julia/src/datatype.c

Line 267 in 64409a0

assert(al % alignm == 0);

will certainly always crash. (also I thought C asserts are disabled in release builds...)
So probably the C part of Julia should still try to emit the correct code adhering to basic SIMD requirements (16bytes), and the constructors on top of that should warn when an incompatible combination is chosen?

[KristofferC]

It seems that the alignment computed in jl_special_vector_alignment is only sometimes used. For a length 32 Float32 vector, the computation would give a 256 byte alignment (verified with printf) but still a 8 byte alignment is used by LLVM. For some reason this change when the vector becomes larger than 32 and what is computed in jl_special_vector_alignment is used..

[m-j-w]

@KristofferC Here's one line that applies the special alignment only conditionally if 'not zero' is returned

julia/src/datatype.c

Line 268 in 64409a0

if (al)

.

And here's another one:

julia/src/ccall.cpp

Line 361 in 080dcf4

return jl_special_vector_alignment(n, ft0) != 0;

(and above that line some further restrictions are made when it actually is a simd type)

[KristofferC]

Yes I know but for a 32 length Float64 vector, 256 is returned, not 0.

[m-j-w]

Have you also seen the definition of MAX_ALIGN at

julia/src/julia.h

Lines 19 to 33 in 080dcf4

	#include <setjmp.h>
	#ifndef _OS_WINDOWS_
	# define jl_jmp_buf sigjmp_buf
	# if defined(_CPU_ARM_) || defined(_CPU_PPC_)
	# define MAX_ALIGN 8
	# elif defined(_CPU_AARCH64_)
	// int128 is 16 bytes aligned on aarch64
	# define MAX_ALIGN 16
	# else
	# define MAX_ALIGN sizeof(void*)
	# endif
	#else
	# define jl_jmp_buf jmp_buf
	# include <malloc.h> //for _resetstkoflw
	# define MAX_ALIGN 8

used in julia_alignment

julia/src/cgutils.cpp

Line 1117 in 080dcf4

static unsigned julia_alignment(Value* /*ptr*/, jl_value_t *jltype, unsigned alignment)

only if no alignment is defined yet?

static unsigned julia_alignment(Value* /*ptr*/, jl_value_t *jltype, unsigned alignment)
{
    if (!alignment && ((jl_datatype_t*)jltype)->layout->alignment > MAX_ALIGN) {
        // Type's natural alignment exceeds strictest alignment promised in heap, so return the heap alignment.
        return MAX_ALIGN;
    }
    return alignment;
}

[KristofferC]

That might explain why the vectors are only 8 byte (64 bit) aligned?