Merkwürdige Effekte durch __restrict



  • Hallo allerseits!

    Ich bin auf ein merkwürdiges Verhalten meines Programmes durch die Benutzung von restrict (bzw. seiner GCC-Erweiterung für C++) gestoßen. Folgender Code führt eine Bereichnung durch, wobei alle 5 beteiligten Arrays garantiert disjunkt sind:

    template <typename T, typename Idx>
    inline void apply(const Idx& idx, 
    		  const T* const __restrict H, 
    		  const T* const __restrict in_re, 
    		  const T* const __restrict in_im, 
    		  T* const __restrict out_re,
    		  T* const __restrict out_im) 
    {
    	for(int j = 1; j <= idx.cols(); ++j) {
    		double re_tmp = 0., im_tmp = 0.;
    		for(int i = idx.lower_idx(j)+1; i <= idx.upper_idx(j); ++i) {
    			const double Hij = H[idx(i, j)];
    			out_re[i-1]	+= Hij * in_im[j-1];
    			out_im[i-1]	-= Hij * in_re[j-1];
    			re_tmp		+= Hij * in_im[i-1];
    			im_tmp		-= Hij * in_re[i-1];
    		}
    		out_re[j-1] += re_tmp;
    		out_im[j-1] += im_tmp;
    
    		const double Hjj = H[idx(j, j)];
    		out_re[j-1] += Hjj * in_im[j-1];
    		out_im[j-1] -= Hjj * in_re[j-1];
    	}
    }
    
    

    Der Code ist absolut performance-kritisch, weshalb ich es für eine gute Idee hielt, __restrict zu verwenden. Tatsächlich führt das dazu, dass andere Zahlen aus der Rechnung resultieren (die sich zwar nur knapp oberhalb vom Maschinenepsilon unterscheiden), als wenn ich kein __restrict nähme. Leider macht das aber im weiteren Verlauf mein Programm gut 100% langsamer.
    Kompiliert wird das ganze bei mir mit -O3 -march=native (insb. mit SSE/AVX2 Instructions). Eigentlich auch noch mit -ffast-math, aber das würde hier nur für Verwirrung sorgen.

    UB kann ich soweit ausschließen, dass es weder Warnings auf dem Level -Wall -Wextra -Wpedantic gibt, noch, dass -fsanitize=address irgendetwas beanstandet hätte. Der Effekt tritt auf, solange auch nur irgendein Funktionsargument mit __restrict kommt, egal welches. Sobald alle __restricts entfernt sind, wird der Code doppelt so schnell, zeigt minimal andere Zahlen und auch der ASM-Output ist bedeutent kürzer:

    	mov	eax, DWORD PTR 4[rdx]
    	test	eax, eax
    	jle	.L224
    	mov	r14d, DWORD PTR 12[rdx]
    	mov	r15d, DWORD PTR 8[rdx]
    	dec	eax
    	lea	ebx, 1[r14+r15]
    	mov	r13d, DWORD PTR [rdx]
    	add	rax, 2
    	movsx	rdx, ebx
    	mov	DWORD PTR -72[rbp], ebx
    	mov	QWORD PTR -80[rbp], rax
    	lea	rbx, 0[0+rdx*8]
    	mov	rdx, QWORD PTR -64[rbp]
    	mov	QWORD PTR -88[rbp], rbx
    	movsx	rbx, r14d
    	lea	r11, [rdx+rbx*8]
    	mov	QWORD PTR -96[rbp], rbx
    	mov	r12, -8
    	xor	ebx, ebx
    	mov	edx, 1
    	.p2align 4,,10
    	.p2align 3
    .L223:
    	mov	eax, edx
    	sub	eax, r14d
    	test	eax, eax
    	mov	r10d, 1
    	cmovle	eax, r10d
    	lea	r9d, [r15+rdx]
    	inc	eax
    	cmp	r9d, r13d
    	cmovg	r9d, r13d
    	cmp	eax, r9d
    	jg	.L267
    	movsx	r10, ebx
    	add	r10, QWORD PTR -96[rbp]
    	lea	r10, [r12+r10*8]
    	cdqe
    	add	r10, QWORD PTR -64[rbp]
    	vmovapd	xmm1, xmm4
    	vmovapd	xmm2, xmm4
    	.p2align 4,,10
    	.p2align 3
    .L222:
    	vmovsd	xmm0, QWORD PTR [r10+rax*8]
    	vmovsd	xmm3, QWORD PTR -8[r8+rdx*8]
    	vfmadd213sd	xmm3, xmm0, QWORD PTR -8[rsi+rax*8]
    	vmovsd	QWORD PTR -8[rsi+rax*8], xmm3
    	vmovsd	xmm3, QWORD PTR -8[rdi+rdx*8]
    	vfnmadd213sd	xmm3, xmm0, QWORD PTR -8[rcx+rax*8]
    	vmovsd	QWORD PTR -8[rcx+rax*8], xmm3
    	vfmadd231sd	xmm2, xmm0, QWORD PTR -8[r8+rax*8]
    	vfnmadd231sd	xmm1, xmm0, QWORD PTR -8[rdi+rax*8]
    	inc	rax
    	cmp	r9d, eax
    	jge	.L222
    .L221:
    	vaddsd	xmm2, xmm2, QWORD PTR -8[rsi+rdx*8]
    	add	ebx, DWORD PTR -72[rbp]
    	sub	r12, 8
    	vmovsd	QWORD PTR -8[rsi+rdx*8], xmm2
    	vaddsd	xmm1, xmm1, QWORD PTR -8[rcx+rdx*8]
    	vmovsd	QWORD PTR -8[rcx+rdx*8], xmm1
    	vmovsd	xmm0, QWORD PTR [r11]
    	vmovsd	xmm1, QWORD PTR -8[r8+rdx*8]
    	add	r11, QWORD PTR -88[rbp]
    	vfmadd213sd	xmm1, xmm0, QWORD PTR -8[rsi+rdx*8]
    	vmovsd	QWORD PTR -8[rsi+rdx*8], xmm1
    	vmovsd	xmm5, QWORD PTR -8[rcx+rdx*8]
    	vfnmadd132sd	xmm0, xmm5, QWORD PTR -8[rdi+rdx*8]
    	vmovsd	QWORD PTR -8[rcx+rdx*8], xmm0
    	inc	rdx
    	cmp	QWORD PTR -80[rbp], rdx
    	jne	.L223
    
    

    gegen mehr als dreimal so viele Instructions mit __restrict:

    	mov	eax, DWORD PTR 4[rdx]
    	test	eax, eax
    	jle	.L228
    	mov	ebx, DWORD PTR [rdx]
    	mov	edi, DWORD PTR 12[rdx]
    	mov	DWORD PTR -240[rbp], ebx
    	mov	ebx, DWORD PTR 8[rdx]
    	mov	DWORD PTR -204[rbp], edi
    	mov	DWORD PTR -248[rbp], ebx
    	lea	ebx, 1[rdi+rbx]
    	movsx	rdx, ebx
    	mov	DWORD PTR -208[rbp], ebx
    	lea	rbx, 0[0+rdx*8]
    	mov	QWORD PTR -256[rbp], rbx
    	dec	eax
    	movsx	rbx, edi
    	mov	rdi, QWORD PTR -176[rbp]
    	add	rax, 2
    	mov	QWORD PTR -224[rbp], rbx
    	mov	QWORD PTR -264[rbp], rax
    	lea	rbx, [rdi+rbx*8]
    	mov	rax, QWORD PTR -336[rbp]
    	mov	rdi, QWORD PTR -232[rbp]
    	mov	QWORD PTR -168[rbp], rbx
    	add	rax, rdi
    	mov	QWORD PTR -280[rbp], rax
    	mov	rax, QWORD PTR -344[rbp]
    	mov	DWORD PTR -132[rbp], 0
    	add	rax, rdi
    	mov	QWORD PTR -288[rbp], rax
    	mov	QWORD PTR -64[rbp], r11
    	mov	r12, r11
    	mov	r11, QWORD PTR -312[rbp]
    	mov	r14, r15
    	mov	r13d, 1
    	mov	r15, r9
    	.p2align 4,,10
    	.p2align 3
    .L227:
    	mov	eax, r13d
    	sub	eax, DWORD PTR -204[rbp]
    	test	eax, eax
    	mov	edi, 1
    	cmovle	eax, edi
    	mov	edi, DWORD PTR -248[rbp]
    	mov	ecx, DWORD PTR -240[rbp]
    	lea	edx, [rdi+r13]
    	cmp	edx, ecx
    	mov	edi, edx
    	cmovg	edi, ecx
    	mov	rdx, QWORD PTR -152[rbp]
    	lea	ebx, 1[rax]
    	mov	DWORD PTR -216[rbp], r13d
    	mov	DWORD PTR -136[rbp], ebx
    	mov	DWORD PTR -160[rbp], edi
    	vmovsd	xmm4, QWORD PTR -8[rdx+r13*8]
    	cmp	ebx, edi
    	jg	.L271
    	sub	edi, eax
    	mov	rax, QWORD PTR -280[rbp]
    	movsx	rcx, ebx
    	lea	rsi, [rax+rcx]
    	lea	r10, 0[0+rsi*8]
    	lea	rbx, -8[r10]
    	mov	rax, QWORD PTR -232[rbp]
    	lea	rdx, [r11+rbx]
    	mov	QWORD PTR -184[rbp], rbx
    	add	r10, 24
    	mov	rbx, QWORD PTR -64[rbp]
    	mov	QWORD PTR -112[rbp], rcx
    	mov	r9d, edi
    	add	rcx, rax
    	lea	r8, [r11+r10]
    	mov	QWORD PTR -88[rbp], r9
    	mov	QWORD PTR -96[rbp], r8
    	add	rbx, 8
    	mov	r8, QWORD PTR -120[rbp]
    	lea	r9, 0[0+rcx*8]
    	mov	QWORD PTR -80[rbp], rbx
    	lea	rbx, -8[r9]
    	add	r9, 24
    	mov	DWORD PTR -72[rbp], edi
    	lea	rdi, [r8+rbx]
    	add	r8, r9
    	mov	QWORD PTR -104[rbp], r8
    	movsx	r8, DWORD PTR -132[rbp]
    	mov	QWORD PTR -192[rbp], rsi
    	mov	rsi, QWORD PTR -224[rbp]
    	mov	QWORD PTR -272[rbp], r8
    	sub	r8, r13
    	mov	rax, QWORD PTR -176[rbp]
    	add	rsi, r8
    	add	rsi, QWORD PTR -112[rbp]
    	sal	rsi, 3
    	mov	r8, rax
    	add	r8, rsi
    	lea	rax, 32[rax+rsi]
    	cmp	rbx, r10
    	lea	rsi, [r11+r9]
    	mov	QWORD PTR -128[rbp], rsi
    	setge	sil
    	cmp	r9, QWORD PTR -184[rbp]
    	setle	r9b
    	or	r9d, esi
    	cmp	QWORD PTR -96[rbp], rdi
    	setbe	sil
    	cmp	rdx, QWORD PTR -104[rbp]
    	setnb	r10b
    	mov	QWORD PTR -200[rbp], rcx
    	lea	rcx, [r11+rbx]
    	mov	ebx, DWORD PTR -72[rbp]
    	or	esi, r10d
    	and	r9d, esi
    	lea	esi, -1[rbx]
    	cmp	esi, 2
    	seta	sil
    	and	r9d, esi
    	cmp	rdi, QWORD PTR -128[rbp]
    	setnb	bl
    	cmp	rcx, QWORD PTR -104[rbp]
    	setnb	sil
    	or	ebx, esi
    	mov	rsi, QWORD PTR -192[rbp]
    	mov	r10, QWORD PTR -80[rbp]
    	add	rsi, QWORD PTR -88[rbp]
    	lea	rsi, -8[r11+rsi*8]
    	and	r9d, ebx
    	cmp	QWORD PTR -64[rbp], rsi
    	setnb	sil
    	cmp	rdx, r10
    	setnb	bl
    	or	esi, ebx
    	and	esi, r9d
    	mov	r9, QWORD PTR -88[rbp]
    	add	r9, QWORD PTR -200[rbp]
    	lea	r9, -8[r11+r9*8]
    	cmp	QWORD PTR -64[rbp], r9
    	setnb	r9b
    	cmp	rcx, r10
    	setnb	bl
    	or	r9d, ebx
    	and	esi, r9d
    	cmp	QWORD PTR -96[rbp], r8
    	setbe	r9b
    	cmp	rdx, rax
    	setnb	r10b
    	or	r9d, r10d
    	test	sil, r9b
    	je	.L222
    	cmp	r8, QWORD PTR -128[rbp]
    	setnb	sil
    	cmp	rcx, rax
    	setnb	r9b
    	or	sil, r9b
    	je	.L222
    	mov	rax, QWORD PTR -112[rbp]
    	mov	rsi, QWORD PTR -120[rbp]
    	add	rax, QWORD PTR -288[rbp]
    	lea	r9, -8[rsi+rax*8]
    	mov	esi, DWORD PTR -72[rbp]
    	vbroadcastsd	ymm5, QWORD PTR -8[r12+r13*8]
    	shr	esi, 2
    	vbroadcastsd	ymm6, xmm4
    	sal	rsi, 5
    	xor	eax, eax
    	vmovapd	xmm1, xmm7
    	vmovapd	xmm3, xmm7
    	.p2align 4,,10
    	.p2align 3
    .L223:
    	vmovupd	ymm0, YMMWORD PTR [r8+rax]
    	vmovapd	ymm2, ymm0
    	vfmadd213pd	ymm2, ymm6, YMMWORD PTR [rcx+rax]
    	vmovupd	YMMWORD PTR [rcx+rax], ymm2
    	vmovapd	ymm2, ymm0
    	vfnmadd213pd	ymm2, ymm5, YMMWORD PTR [rdx+rax]
    	vmovupd	YMMWORD PTR [rdx+rax], ymm2
    	vmulpd	ymm2, ymm0, YMMWORD PTR [r9+rax]
    	vmulpd	ymm0, ymm0, YMMWORD PTR [rdi+rax]
    	add	rax, 32
    	vaddsd	xmm3, xmm3, xmm2
    	vsubsd	xmm1, xmm1, xmm0
    	vunpckhpd	xmm9, xmm2, xmm2
    	vunpckhpd	xmm8, xmm0, xmm0
    	vaddsd	xmm3, xmm3, xmm9
    	vsubsd	xmm1, xmm1, xmm8
    	vextractf128	xmm2, ymm2, 0x1
    	vextractf128	xmm0, ymm0, 0x1
    	vaddsd	xmm3, xmm3, xmm2
    	vsubsd	xmm1, xmm1, xmm0
    	vunpckhpd	xmm2, xmm2, xmm2
    	vunpckhpd	xmm0, xmm0, xmm0
    	vaddsd	xmm3, xmm3, xmm2
    	vsubsd	xmm1, xmm1, xmm0
    	cmp	rax, rsi
    	jne	.L223
    	mov	edi, DWORD PTR -72[rbp]
    	mov	eax, DWORD PTR -136[rbp]
    	mov	edx, edi
    	and	edx, -4
    	add	eax, edx
    	cmp	edi, edx
    	je	.L221
    	mov	ebx, DWORD PTR -216[rbp]
    	mov	r8d, DWORD PTR -204[rbp]
    	mov	edx, eax
    	mov	r10d, DWORD PTR -132[rbp]
    	sub	edx, ebx
    	add	edx, r8d
    	mov	rdi, QWORD PTR -176[rbp]
    	add	edx, r10d
    	movsx	rdx, edx
    	vmovsd	xmm0, QWORD PTR [rdi+rdx*8]
    	movsx	rdx, eax
    	sal	rdx, 3
    	lea	rcx, -8[rdx]
    	lea	rsi, [r15+rcx]
    	vmovapd	xmm2, xmm0
    	vfmadd213sd	xmm2, xmm4, QWORD PTR [rsi]
    	add	rcx, r14
    	vmovsd	QWORD PTR [rsi], xmm2
    	vmovsd	xmm2, QWORD PTR -8[r12+r13*8]
    	vfnmadd213sd	xmm2, xmm0, QWORD PTR [rcx]
    	vmovsd	QWORD PTR [rcx], xmm2
    	mov	rcx, QWORD PTR -152[rbp]
    	vfnmadd231sd	xmm1, xmm0, QWORD PTR -8[r12+rdx]
    	vfmadd231sd	xmm3, xmm0, QWORD PTR -8[rcx+rdx]
    	lea	edx, 1[rax]
    	cmp	DWORD PTR -160[rbp], edx
    	jl	.L221
    	mov	r9d, edx
    	sub	r9d, ebx
    	mov	ecx, r9d
    	add	ecx, r8d
    	add	ecx, r10d
    	movsx	rdx, edx
    	lea	rsi, 0[0+rdx*8]
    	movsx	rcx, ecx
    	vmovsd	xmm0, QWORD PTR [rdi+rcx*8]
    	lea	rcx, -8[rsi]
    	mov	r9, rdi
    	vmovapd	xmm2, xmm4
    	lea	rdi, [r15+rcx]
    	vfmadd213sd	xmm2, xmm0, QWORD PTR [rdi]
    	add	rcx, r14
    	add	eax, 2
    	vmovsd	QWORD PTR [rdi], xmm2
    	vmovsd	xmm2, QWORD PTR -8[r12+r13*8]
    	mov	rdi, QWORD PTR -152[rbp]
    	vfnmadd213sd	xmm2, xmm0, QWORD PTR [rcx]
    	vfmadd231sd	xmm3, xmm0, QWORD PTR -8[rdi+rsi]
    	vmovsd	QWORD PTR [rcx], xmm2
    	vfnmadd231sd	xmm1, xmm0, QWORD PTR -8[r12+rsi]
    	cmp	DWORD PTR -160[rbp], eax
    	jl	.L221
    	sub	eax, ebx
    	add	eax, r8d
    	add	eax, r10d
    	cdqe
    	vmovsd	xmm0, QWORD PTR [r9+rax*8]
    	lea	rcx, [r15+rsi]
    	vmovapd	xmm2, xmm4
    	vfmadd213sd	xmm2, xmm0, QWORD PTR [rcx]
    	mov	rax, rsi
    	add	rax, r14
    	vfmadd231sd	xmm3, xmm0, QWORD PTR [rdi+rdx*8]
    	vmovsd	QWORD PTR [rcx], xmm2
    	vmovsd	xmm2, QWORD PTR -8[r12+r13*8]
    	vfnmadd213sd	xmm2, xmm0, QWORD PTR [rax]
    	vmovsd	QWORD PTR [rax], xmm2
    	vfnmadd231sd	xmm1, xmm0, QWORD PTR [r12+rdx*8]
    .L221:
    	vaddsd	xmm3, xmm3, QWORD PTR -8[r15+r13*8]
    	mov	rax, QWORD PTR -168[rbp]
    	mov	rdi, QWORD PTR -256[rbp]
    	vmovsd	QWORD PTR -8[r15+r13*8], xmm3
    	vaddsd	xmm0, xmm1, QWORD PTR -8[r14+r13*8]
    	mov	ebx, DWORD PTR -208[rbp]
    	vmovsd	QWORD PTR -8[r14+r13*8], xmm0
    	vmovsd	xmm0, QWORD PTR [rax]
    	add	rax, rdi
    	vfmadd213sd	xmm4, xmm0, QWORD PTR -8[r15+r13*8]
    	mov	QWORD PTR -168[rbp], rax
    	mov	rax, QWORD PTR -80[rbp]
    	add	DWORD PTR -132[rbp], ebx
    	mov	QWORD PTR -64[rbp], rax
    	vmovsd	QWORD PTR -8[r15+r13*8], xmm4
    	vmovsd	xmm4, QWORD PTR -8[r14+r13*8]
    	vfnmadd132sd	xmm0, xmm4, QWORD PTR -8[r12+r13*8]
    	vmovsd	QWORD PTR -8[r14+r13*8], xmm0
    	inc	r13
    	cmp	r13, QWORD PTR -264[rbp]
    	jne	.L227
    
    

    Die nahelegende Lösung wäre sicherlich, eben auf __restrict zu pfeifen. Aber ich möchte gerne verstehen, was hier schiefgeht bzw. wie ich das weiter untersuchen könnte. Der ASM-Code sagt mir im Detail jetzt nicht furchtbar viel. Vielen Dank!

    p.s. Ich verwende GCC 9.3.



  • @Jodocus sagte in Merkwürdige Effekte durch __restrict:

    Tatsächlich führt das dazu, dass andere Zahlen aus der Rechnung resultieren (die sich zwar nur knapp oberhalb vom Maschinenepsilon unterscheiden), als wenn ich kein __restrict nähme. Leider macht das aber im weiteren Verlauf mein Programm gut 100% langsamer.

    Hm. Denormals? Die sollten mit -ffast-math eigentlich abgedreht sein. Das funktioniert aber soweit ich weiss nur wenn du das Executable selbst mit -ffast-math baust, und es kann auf jeden Fall passieren dass Denormals nicht deaktiviert sind wenn dein Code als Library in einem anderen Executable läuft.

    Probier vielleicht einfach mal
    _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON) und _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON) vor deinen Berechnungen zu machen. Siehe https://scc.ustc.edu.cn/zlsc/chinagrid/intel/compiler_c/main_cls/GUID-1659EAE1-583E-44EE-BDEA-7C68C46061C7.htm

    Wenns daran liegt, musst du dir aber 'was überlegen. Denn einfach so FTZ und DAZ zu setzen (und nie wieder zu restoren) ist auf jeden Fall "pfui!".



  • Mal eine andere Herangehensweise:
    Wieso übergibst du nicht ein struct und returnst ein struct (dank RTO kein Problem)? Dann dürftest du das Aliasing-Problem nicht haben und der Compiler dürfte besser optimieren können.



  • @Steffo Dir ist schon aufgefallen dass das alles Arras sind, oder?



  • @hustbaer Habe ich mal probiert, aber macht keinen Unterschied. Zur Info: Mein Code läuft nicht als Library, also ich habe volle Kontrolle über alle Linker/Compiler Flags vom Code bis zur Executable.

    Hintergrundinfo: Die Routine wird während einer Adams-Bashforth-Routine zum Lösen einer ODE aufgerufen, mit adaptiver Schrittweite. Je weniger Schritte (= Aufrufe dieser Funktion), desto besser.
    Meine obigen Beobachtungen stehen nach wie vor. Mit restrict sind die Zahlen leicht anders und am Ende werden ca. 100.000 Aufrufe auf die Funktion gemacht. Nehme ich restrict raus, sind es nur noch 50.000 Aufrufe. Die Ergebnisse sind in beiden Fällen bis auf die eingestellte Genauigkeit der Konvergenz des ODE-Solvers korrekt.

    Zu allem Überfluss kommt jetzt hinzu: wenn ich am Code leicht etwas ändere, z.B. die Zeilen mit Hjj am Anfang statt am Ende der äußeren Loop setze (was mathematisch zwar äquivalent ist), also

    	for(int j = 1; j <= idx.cols(); ++j) {
    		const double Hjj = H[idx(j, j)]; // *
    		out_re[j-1] += Hjj * in_im[j-1]; // *
    		out_im[j-1] -= Hjj * in_re[j-1]; // *
    
    		double re_tmp = 0., im_tmp = 0.;
    		for(int i = idx.lower_idx(j)+1; i <= idx.upper_idx(j); ++i) {
    			const double Hij = H[idx(i, j)];
    			out_re[i-1]	+= Hij * in_im[j-1];
    			out_im[i-1]	-= Hij * in_re[j-1];
    			re_tmp		+= Hij * in_im[i-1];
    			im_tmp		-= Hij * in_re[i-1];
    		}
    		out_re[j-1] += re_tmp;
    		out_im[j-1] += im_tmp;
    	}
    
    

    Dann ist das Ergebnis von den __restricts völlig unabhängig (wie man es ja auch erwarten würde)! Mir ist schon klar, dass für die Fließpunktkommazahlen Addition nicht assoziativ/kommutativ sein könnte (die Reihenfolge also technisch, wenn auch nicht mathematisch eine Rolle spielt), aber ich habe explizit -ffast-math ausgeschaltet (und mehrfach überprüft und neugebaut), sodass Optimierungen, die der Compiler durch __restrict machen darf, nicht z.B. an dieser Reihenfolge rütteln dürfte (da sonst nicht mehr IEEE-konform). Warum das in einem Fall doch passiert, riecht für mich irgendwie schon arg nach Compiler-Bug.

    Ganz abgesehen davon, dass die ODE-Routine numerisch alles andere als stabil ist, wenn solche Effekte zu 100% Laufzeitsteigerung führen können...



  • Hm, ich hatte deinen Code mal in Godbolt gekippt und habe diese großen Assember-Unterschiede nicht gesehen. Jeweils 95 Zeilen. Allerdings musste ich dazu natürlich das template und das inline rausnehmen und mir eine Fake-Index-Klasse machen. Vielleicht tritt der Effekt, den du siehst, nur bei bestimmten Klassen für Idx und T auf - vielleicht trifft der Compiler irgendwelche Optimierungen? Oder das Problem kommt von irgendwo anders her? Finde ich schon komisch, dass die Änderung im Epsilon-Bereich gleich doppelt so viele Calls zur Folge haben soll. Dann würde ich denken, dass der Solver nicht korrekt arbeitet - oder dass du irgendeinen Spezialfall triffst, wo das tatsächlich große Relevanz hat.



  • @Jodocus sagte in Merkwürdige Effekte durch __restrict:

    Zu allem Überfluss kommt jetzt hinzu: wenn ich am Code leicht etwas ändere, z.B. die Zeilen mit Hjj am Anfang statt am Ende der äußeren Loop setze (was mathematisch zwar äquivalent ist), also

    Wie groß sind die Felder? Der Loop über i liest vier Felder auf einmal, wenn diese groß sind, ist das schlecht für die Cache Hits.



  • @john-0 sagte in Merkwürdige Effekte durch __restrict:

    @Jodocus sagte in Merkwürdige Effekte durch __restrict:

    Zu allem Überfluss kommt jetzt hinzu: wenn ich am Code leicht etwas ändere, z.B. die Zeilen mit Hjj am Anfang statt am Ende der äußeren Loop setze (was mathematisch zwar äquivalent ist), also

    Wie groß sind die Felder? Der Loop über i liest vier Felder auf einmal, wenn diese groß sind, ist das schlecht für die Cache Hits.

    Das variiert, ich kann da keine Rücksicht auf die Cachegröße nehmen.
    Ich glaube, ich mache hier Schluss. Es ist nicht ausgeschlossen, dass ich noch das eine oder andere an dem Code ändern muss, und da der Effekt so empfindlich ist (und evtl. am Ende gar nicht mehr auftritt), lohnt sich die Nachforschung an der Stelle nicht. Es bleibt nur ein komischer Nachgeschmack. Trotzdem danke an alle für die Empfehlungen!



  • @Jodocus sagte in Merkwürdige Effekte durch __restrict:

    Das variiert, ich kann da keine Rücksicht auf die Cachegröße nehmen.

    Das solltest Du aber, weil das zum Teil drastische Auswirkungen hat. Die Optimierung auf Cache Hits ist bei modernen Computer bei HPC die wichtigste Optimierung, die man machen kann. Es nützt Dir nichts einige wenige Taktzyklen einzusparen, wenn Du am Ende tausende Zyklen auf das Laden aus dem RAM warten musst.



  • @Jodocus
    Im ersten Beitrag hast du geschrieben du baust mit -ffast-math und nun schreibst du du baust mit -ffast-math ausgeschaltet. Jetzt bin ich verwirrt 🙂
    Meinst du du hast beide Varianten ( -ffast-math und -fno-fast-math) probiert?

    Mit restrict sind die Zahlen leicht anders und am Ende werden ca. 100.000 Aufrufe auf die Funktion gemacht. Nehme ich restrict raus, sind es nur noch 50.000 Aufrufe.

    OK. Das heisst die Funktion läuft gar nicht unbedingt langsamer, sie erzeugt nur leicht andere Ergebnisse so dass die Werte langsamer konvergieren und du mehr Iterationen brauchst. Das kann mit -ffast-math natürlich passieren.

    Ich bin mir nicht ganz sicher, aber ich denke dass es auch mit -fno-fast-math Unterschiede geben kann. Und zwar ist bei sowas wie a = x * y + z glaube ich nicht definiert ob das Zwischenergebnis x * y erstmal gerundet wird vor der Addition oder nicht. Der Compiler generiert für deinen Code (mit und ohne restrict) fused-multiply-add Instruktionen (die ganzen vfmaddXXX und vfnmaddXXX Dinger), und diese runden das Zwischenergebnis nicht. Da kann ich mir schon vorstellen das was anderes rauskommt wenn z.B. mehr oder weniger fused-multiply-add generiert werden, abhängig von den Annahmen die der Compiler durch restrict oder nicht-restrict treffen kann.

    Was vielleicht auch noch einen Versucht wert wäre: du könntest std::fma in deinem Code verwenden: https://en.cppreference.com/w/cpp/numeric/math/fma
    Dann müsste der Compiler immer sicherstellen dass die Zwischenergebnisse nicht gerundet werden, da std::fma das garantiert.
    Dann bliebe nur noch die Frage ob std::fma zu mehr oder weniger Iterationen führt, verglichen mit dem ursprünglichen Code ohne restrict.


Anmelden zum Antworten