Есть планы написать виртуальную машину для исполнения шейдеров? :-)
В ближайшее время таких планов точно нет, возможно в будущем. Было бы интересно над этим поработать
Объясните пожалуйста, кто отломал от вашего компилятора циклы?
matrix4x4 matrix4x4::operator*(matrix4x4 const &m) const
{
return matrix4x4{
values[0][0] * m.values[0][0] + values[0][1] * m.values[1][0] + values[0][2] * m.values[2][0] + values[0][3] * m.values[3][0],
values[0][0] * m.values[0][1] + values[0][1] * m.values[1][1] + values[0][2] * m.values[2][1] + values[0][3] * m.values[3][1],
values[0][0] * m.values[0][2] + values[0][1] * m.values[1][2] + values[0][2] * m.values[2][2] + values[0][3] * m.values[3][2],
values[0][0] * m.values[0][3] + values[0][1] * m.values[1][3] + values[0][2] * m.values[2][3] + values[0][3] * m.values[3][3],
values[1][0] * m.values[0][0] + values[1][1] * m.values[1][0] + values[1][2] * m.values[2][0] + values[1][3] * m.values[3][0],
values[1][0] * m.values[0][1] + values[1][1] * m.values[1][1] + values[1][2] * m.values[2][1] + values[1][3] * m.values[3][1],
values[1][0] * m.values[0][2] + values[1][1] * m.values[1][2] + values[1][2] * m.values[2][2] + values[1][3] * m.values[3][2],
values[1][0] * m.values[0][3] + values[1][1] * m.values[1][3] + values[1][2] * m.values[2][3] + values[1][3] * m.values[3][3],
values[2][0] * m.values[0][0] + values[2][1] * m.values[1][0] + values[2][2] * m.values[2][0] + values[2][3] * m.values[3][0],
values[2][0] * m.values[0][1] + values[2][1] * m.values[1][1] + values[2][2] * m.values[2][1] + values[2][3] * m.values[3][1],
values[2][0] * m.values[0][2] + values[2][1] * m.values[1][2] + values[2][2] * m.values[2][2] + values[2][3] * m.values[3][2],
values[2][0] * m.values[0][3] + values[2][1] * m.values[1][3] + values[2][2] * m.values[2][3] + values[2][3] * m.values[3][3],
values[3][0] * m.values[0][0] + values[3][1] * m.values[1][0] + values[3][2] * m.values[2][0] + values[3][3] * m.values[3][0],
values[3][0] * m.values[0][1] + values[3][1] * m.values[1][1] + values[3][2] * m.values[2][1] + values[3][3] * m.values[3][1],
values[3][0] * m.values[0][2] + values[3][1] * m.values[1][2] + values[3][2] * m.values[2][2] + values[3][3] * m.values[3][2],
values[3][0] * m.values[0][3] + values[3][1] * m.values[1][3] + values[3][2] * m.values[2][3] + values[3][3] * m.values[3][3]};
}

Если вы так сделали, думая что компилятор — идиот и не умеет разворачивать циклы, а затем реализовывать их в векторных инструкциях в глубоко ошибаетесь.
Задам вопрос иначе — а менее уродливого способа работы с матрицами разве нет? Разве красивое использование синтаксиса в дальнейшем (всего-то две операции будут выглядить красиво — сложение и умножение, для остальных придется лепить методы) стоит появления в исходнике таких кирпичей?
И в чем проблема молча все поправить, оттестить и отправить коммит на гитхабе?
Статья отличная!
Как говорил один сантехник, которого потом правда расстреляли, — «Тут систему менять надо» (в смысле — все переписывать).
Что касается контента статьи — ну учебник Немнюгина содержит все то же самое, правда издан 14 лет назад, в 2000 году.

Что касается избавления от ужаса, на данный момент я заметил следующее:
Авторский кирпич транслируется в 204 инструкции:
Скрытый текст
.LFB3:
	.cfi_startproc
	vmovss	56(%rsi), %xmm7
	movq	%rdi, %rax
	vmovss	60(%rsi), %xmm6
	vmovss	60(%rdx), %xmm0
	vmovss	8(%rdx), %xmm12
	vmovss	24(%rdx), %xmm13
	vmovss	48(%rsi), %xmm2
	vmovss	12(%rdx), %xmm3
	vmovss	52(%rsi), %xmm4
	vmovss	28(%rdx), %xmm5
	vmovss	44(%rdx), %xmm1
	vmovss	40(%rdx), %xmm10
	vmovss	56(%rdx), %xmm11
	vmovss	%xmm7, -8(%rsp)
	vmovss	%xmm6, -4(%rsp)
	vmovss	4(%rdx), %xmm7
	vmovss	20(%rdx), %xmm6
	vmovss	%xmm0, -28(%rsp)
	vmovss	%xmm12, -48(%rsp)
	vmovss	36(%rdx), %xmm0
	vmovss	52(%rdx), %xmm12
	vmovss	%xmm13, -44(%rsp)
	vmovss	32(%rdx), %xmm13
	vmovss	%xmm2, -24(%rsp)
	vmovss	%xmm3, -20(%rsp)
	vmovss	%xmm4, -16(%rsp)
	vmovss	%xmm5, -12(%rsp)
	vmovss	%xmm1, -32(%rsp)
	vmovss	%xmm10, -72(%rsp)
	vmovss	(%rdx), %xmm1
	vmovss	%xmm11, -68(%rsp)
	vmovss	%xmm7, -64(%rsp)
	vmovss	%xmm6, -60(%rsp)
	vmovss	%xmm0, -56(%rsp)
	vmovss	%xmm12, -52(%rsp)
	vmovss	16(%rdx), %xmm0
	vmovss	%xmm13, -40(%rsp)
	vmovss	(%rsi), %xmm13
	vmovss	4(%rsi), %xmm12
	vmulss	%xmm13, %xmm1, %xmm15
	vmovss	48(%rdx), %xmm14
	vmovss	%xmm14, -36(%rsp)
	vmulss	%xmm12, %xmm0, %xmm14
	vmovss	8(%rsi), %xmm11
	vmovss	12(%rsi), %xmm10
	vmovss	16(%rsi), %xmm9
	vaddss	%xmm14, %xmm15, %xmm15
	vmulss	-40(%rsp), %xmm11, %xmm14
	vmovss	20(%rsi), %xmm8
	vmovss	24(%rsi), %xmm7
	vmovss	28(%rsi), %xmm6
	vmovss	32(%rsi), %xmm5
	vmovss	36(%rsi), %xmm4
	vaddss	%xmm14, %xmm15, %xmm15
	vmulss	-36(%rsp), %xmm10, %xmm14
	vmovss	40(%rsi), %xmm3
	vmovss	44(%rsi), %xmm2
	vaddss	%xmm14, %xmm15, %xmm14
	vmulss	-64(%rsp), %xmm13, %xmm15
	vmovss	%xmm14, (%rdi)
	vmulss	-60(%rsp), %xmm12, %xmm14
	vaddss	%xmm14, %xmm15, %xmm15
	vmulss	-56(%rsp), %xmm11, %xmm14
	vaddss	%xmm14, %xmm15, %xmm15
	vmulss	-52(%rsp), %xmm10, %xmm14
	vaddss	%xmm14, %xmm15, %xmm14
	vmulss	-48(%rsp), %xmm13, %xmm15
	vmovss	%xmm14, 4(%rdi)
	vmulss	-44(%rsp), %xmm12, %xmm14
	vaddss	%xmm14, %xmm15, %xmm15
	vmulss	-72(%rsp), %xmm11, %xmm14
	vmulss	-32(%rsp), %xmm11, %xmm11
	vaddss	%xmm14, %xmm15, %xmm15
	vmulss	-68(%rsp), %xmm10, %xmm14
	vmulss	-28(%rsp), %xmm10, %xmm10
	vaddss	%xmm14, %xmm15, %xmm14
	vmovss	-20(%rsp), %xmm15
	vmulss	%xmm15, %xmm13, %xmm13
	vmovss	%xmm14, 8(%rdi)
	vmovss	-12(%rsp), %xmm14
	vmulss	%xmm14, %xmm12, %xmm12
	vaddss	%xmm12, %xmm13, %xmm13
	vmovaps	%xmm15, %xmm12
	vaddss	%xmm11, %xmm13, %xmm13
	vmulss	%xmm9, %xmm1, %xmm11
	vaddss	%xmm10, %xmm13, %xmm13
	vmulss	%xmm8, %xmm0, %xmm10
	vaddss	%xmm10, %xmm11, %xmm11
	vmulss	-40(%rsp), %xmm7, %xmm10
	vmovss	%xmm13, 12(%rdi)
	vmovaps	%xmm14, %xmm13
	vaddss	%xmm10, %xmm11, %xmm11
	vmulss	-36(%rsp), %xmm6, %xmm10
	vaddss	%xmm10, %xmm11, %xmm10
	vmovss	%xmm10, 16(%rdi)
	vmulss	-60(%rsp), %xmm8, %xmm10
	vmulss	-64(%rsp), %xmm9, %xmm11
	vaddss	%xmm10, %xmm11, %xmm11
	vmulss	-56(%rsp), %xmm7, %xmm10
	vaddss	%xmm10, %xmm11, %xmm11
	vmulss	-52(%rsp), %xmm6, %xmm10
	vaddss	%xmm10, %xmm11, %xmm10
	vmulss	-48(%rsp), %xmm9, %xmm11
	vmulss	%xmm15, %xmm9, %xmm9
	vmovss	-32(%rsp), %xmm15
	vmovss	%xmm10, 20(%rdi)
	vmulss	-44(%rsp), %xmm8, %xmm10
	vmulss	%xmm14, %xmm8, %xmm8
	vmovss	-48(%rsp), %xmm14
	vaddss	%xmm8, %xmm9, %xmm9
	vmovss	-40(%rsp), %xmm8
	vaddss	%xmm10, %xmm11, %xmm11
	vmulss	-72(%rsp), %xmm7, %xmm10
	vmulss	%xmm15, %xmm7, %xmm7
	vaddss	%xmm7, %xmm9, %xmm9
	vmulss	%xmm5, %xmm1, %xmm7
	vaddss	%xmm10, %xmm11, %xmm11
	vmulss	-68(%rsp), %xmm6, %xmm10
	vaddss	%xmm10, %xmm11, %xmm10
	vmovss	-28(%rsp), %xmm11
	vmulss	%xmm11, %xmm6, %xmm6
	vmovss	%xmm10, 24(%rdi)
	vaddss	%xmm6, %xmm9, %xmm9
	vmulss	%xmm4, %xmm0, %xmm6
	vmovss	-44(%rsp), %xmm10
	vaddss	%xmm6, %xmm7, %xmm7
	vmulss	%xmm8, %xmm3, %xmm6
	vmovss	%xmm9, 28(%rdi)
	vmovss	-36(%rsp), %xmm9
	vaddss	%xmm6, %xmm7, %xmm7
	vmulss	%xmm9, %xmm2, %xmm6
	vaddss	%xmm6, %xmm7, %xmm6
	vmulss	-64(%rsp), %xmm5, %xmm7
	vmovss	%xmm6, 32(%rdi)
	vmulss	-60(%rsp), %xmm4, %xmm6
	vaddss	%xmm6, %xmm7, %xmm7
	vmulss	-56(%rsp), %xmm3, %xmm6
	vaddss	%xmm6, %xmm7, %xmm7
	vmulss	-52(%rsp), %xmm2, %xmm6
	vaddss	%xmm6, %xmm7, %xmm6
	vmulss	%xmm14, %xmm5, %xmm7
	vmulss	%xmm12, %xmm5, %xmm5
	vmovss	%xmm6, 36(%rdi)
	vmulss	%xmm10, %xmm4, %xmm6
	vmulss	%xmm13, %xmm4, %xmm4
	vaddss	%xmm6, %xmm7, %xmm7
	vmulss	-72(%rsp), %xmm3, %xmm6
	vaddss	%xmm4, %xmm5, %xmm5
	vmulss	%xmm15, %xmm3, %xmm3
	vmovss	-16(%rsp), %xmm4
	vmulss	%xmm4, %xmm0, %xmm0
	vaddss	%xmm3, %xmm5, %xmm5
	vmovss	-64(%rsp), %xmm3
	vaddss	%xmm6, %xmm7, %xmm7
	vmulss	-68(%rsp), %xmm2, %xmm6
	vmulss	%xmm11, %xmm2, %xmm2
	vaddss	%xmm2, %xmm5, %xmm5
	vmovss	-24(%rsp), %xmm2
	vmulss	%xmm2, %xmm1, %xmm1
	vaddss	%xmm6, %xmm7, %xmm6
	vmovaps	%xmm12, %xmm7
	vmovaps	%xmm15, %xmm12
	vmovaps	%xmm8, %xmm15
	vmovss	-8(%rsp), %xmm8
	vaddss	%xmm0, %xmm1, %xmm1
	vmulss	%xmm8, %xmm15, %xmm0
	vmovss	%xmm5, 44(%rdi)
	vmovaps	%xmm9, %xmm15
	vmovss	%xmm6, 40(%rdi)
	vmovss	-4(%rsp), %xmm9
	vmovss	-60(%rsp), %xmm5
	vmovaps	%xmm13, %xmm6
	vaddss	%xmm0, %xmm1, %xmm1
	vmulss	%xmm9, %xmm15, %xmm0
	vmovaps	%xmm11, %xmm13
	vaddss	%xmm0, %xmm1, %xmm1
	vmulss	%xmm2, %xmm3, %xmm0
	vmovss	-56(%rsp), %xmm3
	vmovss	%xmm1, 48(%rdi)
	vmulss	%xmm4, %xmm5, %xmm1
	vmovss	-52(%rsp), %xmm5
	vaddss	%xmm1, %xmm0, %xmm0
	vmulss	%xmm8, %xmm3, %xmm1
	vaddss	%xmm1, %xmm0, %xmm0
	vmulss	%xmm9, %xmm5, %xmm1
	vaddss	%xmm1, %xmm0, %xmm0
	vmulss	%xmm2, %xmm14, %xmm1
	vmovss	%xmm0, 52(%rdi)
	vmulss	%xmm4, %xmm10, %xmm0
	vmovss	-72(%rsp), %xmm10
	vmovss	-68(%rsp), %xmm11
	vaddss	%xmm0, %xmm1, %xmm0
	vmulss	%xmm8, %xmm10, %xmm1
	vaddss	%xmm1, %xmm0, %xmm0
	vmulss	%xmm9, %xmm11, %xmm1
	vaddss	%xmm1, %xmm0, %xmm0
	vmulss	%xmm2, %xmm7, %xmm1
	vmovss	%xmm0, 56(%rdi)
	vmulss	%xmm4, %xmm6, %xmm0
	vaddss	%xmm0, %xmm1, %xmm0
	vmulss	%xmm8, %xmm12, %xmm1
	vaddss	%xmm1, %xmm0, %xmm0
	vmulss	%xmm9, %xmm13, %xmm1
	vaddss	%xmm1, %xmm0, %xmm0
	vmovss	%xmm0, 60(%rdi)
	ret
	.cfi_endproc


Банальное и очевидное решение в три цикла:
Скрытый текст
struct matint
{
    float data[4][4];
};

class matrix4x4
{
public:
    matrix4x4();
    matrix4x4(const matint& data);
    matrix4x4 operator*(matrix4x4 const& m) const;
    void randInit();
    void print() const;
private:
    matint values;

};


matrix4x4 matrix4x4::operator*(matrix4x4 const &m) const
{
    matint out;
    for(size_t i=0;i<4;i++)
    {
        for(size_t j=0;j<4;j++)
        {
            out.data[j][i]=0.0;
            for(size_t k=0;k<4;k++)
            {
                out.data[j][i]+=values.data[k][i]*m.values.data[j][k];
            }
        }
    }
    return{out};
}



Транслируется в 104 инструкции:

Скрытый текст
	.cfi_startproc
	pushq	%rbp
	.cfi_def_cfa_offset 16
	.cfi_offset 6, -16
	movq	%rdi, %rax
	movq	%rsp, %rbp
	.cfi_def_cfa_register 6
	andq	$-32, %rsp
	addq	$16, %rsp
	vmovss	60(%rdx), %xmm4
	vmovss	16(%rdx), %xmm12
	vmovss	32(%rdx), %xmm8
	vshufps	$0, %xmm12, %xmm12, %xmm12
	vmovss	%xmm4, -112(%rsp)
	vshufps	$0, %xmm8, %xmm8, %xmm8
	vmovups	(%rsi), %xmm4
	vbroadcastss	(%rdx), %xmm14
	vmulps	%xmm12, %xmm4, %xmm12
	vmovss	56(%rdx), %xmm3
	vmulps	%xmm14, %xmm4, %xmm14
	vmovss	%xmm3, -108(%rsp)
	vmovups	16(%rsi), %xmm2
	vxorps	%xmm3, %xmm3, %xmm3
	vmulps	%xmm8, %xmm4, %xmm8
	vmovss	20(%rdx), %xmm11
	vmovss	36(%rdx), %xmm7
	vbroadcastss	4(%rdx), %xmm15
	vshufps	$0, %xmm11, %xmm11, %xmm11
	vshufps	$0, %xmm7, %xmm7, %xmm7
	vmulps	%xmm15, %xmm2, %xmm15
	vaddps	%xmm3, %xmm14, %xmm14
	vmovss	52(%rdx), %xmm1
	vaddps	%xmm3, %xmm12, %xmm12
	vmulps	%xmm2, %xmm11, %xmm11
	vmovss	24(%rdx), %xmm10
	vaddps	%xmm3, %xmm8, %xmm8
	vmulps	%xmm2, %xmm7, %xmm7
	vmovss	40(%rdx), %xmm6
	vmovss	%xmm1, -104(%rsp)
	vshufps	$0, %xmm10, %xmm10, %xmm10
	vmovups	32(%rsi), %xmm1
	vshufps	$0, %xmm6, %xmm6, %xmm6
	vaddps	%xmm14, %xmm15, %xmm14
	vbroadcastss	8(%rdx), %xmm13
	vmulps	%xmm10, %xmm1, %xmm10
	vaddps	%xmm7, %xmm8, %xmm7
	vmovss	28(%rdx), %xmm9
	vmulps	%xmm13, %xmm1, %xmm13
	vaddps	%xmm11, %xmm12, %xmm11
	vmovss	44(%rdx), %xmm5
	vmulps	%xmm6, %xmm1, %xmm6
	vmovss	48(%rdx), %xmm0
	vmovss	%xmm0, -100(%rsp)
	vshufps	$0, %xmm9, %xmm9, %xmm9
	vmovups	48(%rsi), %xmm0
	vshufps	$0, %xmm5, %xmm5, %xmm5
	vaddps	%xmm14, %xmm13, %xmm13
	vbroadcastss	12(%rdx), %xmm14
	vmulps	%xmm9, %xmm0, %xmm9
	vaddps	%xmm6, %xmm7, %xmm6
	vmulps	%xmm0, %xmm14, %xmm14
	vaddps	%xmm10, %xmm11, %xmm10
	vmulps	%xmm5, %xmm0, %xmm5
	vaddps	%xmm14, %xmm13, %xmm15
	vaddps	%xmm5, %xmm6, %xmm5
	vaddps	%xmm9, %xmm10, %xmm9
	vmovaps	%xmm15, -96(%rsp)
	movq	-96(%rsp), %rdi
	vmovaps	%xmm15, -80(%rsp)
	vmovaps	%xmm9, -96(%rsp)
	movq	-96(%rsp), %rsi
	vmovaps	%xmm9, -64(%rsp)
	movq	%rdi, (%rax)
	movq	-72(%rsp), %rdi
	vmovaps	%xmm5, -96(%rsp)
	vmovaps	%xmm5, -48(%rsp)
	vbroadcastss	-100(%rsp), %xmm5
	movq	-96(%rsp), %rcx
	movq	%rsi, 16(%rax)
	vmulps	%xmm5, %xmm4, %xmm4
	movq	-56(%rsp), %rsi
	movq	%rdi, 8(%rax)
	movq	%rcx, 32(%rax)
	movq	-40(%rsp), %rcx
	movq	%rsi, 24(%rax)
	vaddps	%xmm3, %xmm4, %xmm3
	movq	%rcx, 40(%rax)
	vbroadcastss	-104(%rsp), %xmm4
	vmulps	%xmm4, %xmm2, %xmm2
	vaddps	%xmm2, %xmm3, %xmm2
	vbroadcastss	-108(%rsp), %xmm3
	vmulps	%xmm3, %xmm1, %xmm1
	vaddps	%xmm1, %xmm2, %xmm1
	vbroadcastss	-112(%rsp), %xmm2
	vmulps	%xmm2, %xmm0, %xmm0
	vaddps	%xmm0, %xmm1, %xmm7
	vmovaps	%xmm7, -96(%rsp)
	movq	-96(%rsp), %rdx
	vmovaps	%xmm7, -32(%rsp)
	movq	%rdx, 48(%rax)
	movq	-24(%rsp), %rdx
	movq	%rdx, 56(%rax)
	leave
	.cfi_def_cfa 7, 8
	ret
	.cfi_endproc



Что я упускаю?
А что показывает секундомер?
/me ушел писать тесты.
Мы когда-то в школе ray tracing в Borland C 3.1 писали. Вектора и матрицы множили такими же конструкциями, но у нас значения лежали не в массиве, а в мемберах структуры вида m11, m12, m13… Но это было давно и неправда, и скорее всего тогда компилятор был сильно глупее.

Кстати, если заменить массив на именованные мемберы, сколько инструкций будет и что на секундомере?
У меня до сих пор хранится книжка Шикина и Борескова 97-го года издания, вся исписанная на полях карандашными комментариями. До сих пор восхищаюсь изяществом, с каким там были написаны классы матриц, векторов, 3D-объектов, перегрузка операторов для работы с ними.
image
Упускаете бессмысленность сравнения производительности двух вариантов по количеству инструкций, а не по таймингу.
Прогнал тесты. Итого:
Скрытый текст
#define itemCount (20000000)
 void testMy()
{
    //Сброс ГСЧ, чтобы числовые значения всегда повторялись
    srand(0);

    high_resolution_clock::time_point t1 = high_resolution_clock::now();
    matrix4x4* out=new matrix4x4[itemCount];
    matrix4x4 a;
    out[0].randInit();
    a.Init9999999();
    for(size_t i=1;i<itemCount;i++)
    {
        out[i]=a*out[i-1];
    }
    high_resolution_clock::time_point t2 = high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>( t2 - t1 ).count();
    std::cout<< duration <<std::endl;
    out[itemCount-1].print();
    delete[] out;
}


Дает итог (значения меняются незначительно):
698649  -- время работы кода с циклами в микросекундах
matrix(
     [0.840188, 0.911647, 0.277775, 0.364784]
, [0.394383, 0.197551, 0.55397, 0.513401]
, [0.783099, 0.335223, 0.477397, 0.95223]
, [0.79844, 0.76823, 0.628871, 0.916195]
);
1091429  -- время работы расписанного кода без циклов в микросекундах
matrix(
     [0.840188, 0.911647, 0.277775, 0.364784]
, [0.394383, 0.197551, 0.55397, 0.513401]
, [0.783099, 0.335223, 0.477397, 0.95223]
, [0.79844, 0.76823, 0.628871, 0.916195]
);


Разница почти в полтора-два раза в пользу циклов. Исходники на гитхабе.

ВАЖНО: При компиляции включить третий уровень оптимизации и AVX
Спасибо
А если без AVX? Просто мне тут для игры (Вангеры, там много софт рендринга) пришлось sse2 отключать так как были люди у которых его не было. :(
Так как после тех тестов я пару раз перезагрузил машину и обновил ядро, тест с AVX я переделал:
image
Верхние три строки получены на машине с процессором:
 $ cat /proc/cpuinfo
processor	: 0
vendor_id	: GenuineIntel
cpu family	: 6
model		: 58
model name	: Intel(R) Core(TM) i7-3517U CPU @ 1.90GHz
stepping	: 9
microcode	: 0x1b
cpu MHz		: 849.468
cache size	: 4096 KB
physical id	: 0
siblings	: 4
core id		: 0
cpu cores	: 2
apicid		: 0
initial apicid	: 0
fpu		: yes
fpu_exception	: yes
cpuid level	: 13
wp		: yes
flags		: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx rdtscp lm constant_tsc arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc aperfmperf eagerfpu pni pclmulqdq dtes64 monitor ds_cpl vmx est tm2 ssse3 cx16 xtpr pdcm pcid sse4_1 sse4_2 x2apic popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm ida arat epb xsaveopt pln pts dtherm tpr_shadow vnmi flexpriority ept vpid fsgsbase smep erms
bogomips	: 4788.82
clflush size	: 64
cache_alignment	: 64
address sizes	: 36 bits physical, 48 bits virtual



Значения в самой нижней строке — результат работы программы на настоящем Pentium III/1GHz. В этом тесте умножалось не двадцать миллионов матриц, а всего два миллиона, иначе бы старикану не хватило памяти.
Увидеть характеристики процессора музейного экспоната:
 $ cat /proc/cpuinfo
processor       : 0
vendor_id       : GenuineIntel
cpu family      : 6
model           : 8
model name      : Pentium III (Coppermine)
stepping        : 10
cpu MHz         : 997.535
cache size      : 256 KB
fdiv_bug        : no
hlt_bug         : no
f00f_bug        : no
coma_bug        : no
fpu             : yes
fpu_exception   : yes
cpuid level     : 2
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 sep mtrr pge mca cmov pse36 mmx fxsr sse up
bogomips        : 1995.07
clflush size    : 32
cache_alignment : 32
address sizes   : 36 bits physical, 32 bits virtual


На самом деле ваша вторая строчка это SSE2, поскольку x86-64 по определению включает в себя SSE2.
Спасибо за замечание, я с вами полностью согласен, в ассемблере 64 разрядной версии стоят инструкции SSE.
Хоть я и согласен с вами, но я думаю что вопросы импрувмента кода лучше обсуждать лично, если конечно статья не о написании идеального кода :)

PS. Ну, в ручном режиме можно и до 68-ми инструкций дожать ;) (а если еще и память выровнять, да movaps везде....)
См. под кат
		movups xmm0, [eax]
		movups xmm1, [eax + 16]
		movups xmm2, [eax + 32]
		movups xmm3, [eax + 48]
		movups xmm4, [ebx]
		movaps xmm5, xmm4
		movaps xmm6, xmm4
		movaps xmm7, xmm4
		shufps xmm4, xmm4, 00000000b
		shufps xmm5, xmm5, 01010101b
		shufps xmm6, xmm6, 10101010b
		shufps xmm7, xmm7, 11111111b
		mulps xmm4, xmm0
		mulps xmm5, xmm1
		mulps xmm6, xmm2
		mulps xmm7, xmm3
		addps xmm6, xmm7
		addps xmm4, xmm5
		addps xmm4, xmm6
		movups [ecx], xmm4
		movups xmm4, [ebx + 16]
		movaps xmm5, xmm4
		movaps xmm6, xmm4
		movaps xmm7, xmm4
		shufps xmm4, xmm4, 00000000b
		shufps xmm5, xmm5, 01010101b
		shufps xmm6, xmm6, 10101010b
		shufps xmm7, xmm7, 11111111b
		mulps xmm4, xmm0
		mulps xmm5, xmm1
		mulps xmm6, xmm2
		mulps xmm7, xmm3
		addps xmm6, xmm7
		addps xmm4, xmm5
		addps xmm4, xmm6
		movups [ecx + 16], xmm4
		movups xmm4, [ebx + 32]
		movaps xmm5, xmm4
		movaps xmm6, xmm4
		movaps xmm7, xmm4
		shufps xmm4, xmm4, 00000000b
		shufps xmm5, xmm5, 01010101b
		shufps xmm6, xmm6, 10101010b
		shufps xmm7, xmm7, 11111111b
		mulps xmm4, xmm0
		mulps xmm5, xmm1
		mulps xmm6, xmm2
		mulps xmm7, xmm3
		addps xmm6, xmm7
		addps xmm4, xmm5
		addps xmm4, xmm6
		movups [ecx + 32], xmm4
		movups xmm4, [ebx + 48]
		movaps xmm5, xmm4
		movaps xmm6, xmm4
		movaps xmm7, xmm4
		shufps xmm4, xmm4, 00000000b
		shufps xmm5, xmm5, 01010101b
		shufps xmm6, xmm6, 10101010b
		shufps xmm7, xmm7, 11111111b
		mulps xmm4, xmm0
		mulps xmm5, xmm1
		mulps xmm6, xmm2
		mulps xmm7, xmm3
		addps xmm6, xmm7
		addps xmm4, xmm5
		addps xmm4, xmm6
		movups [ecx + 48], xmm4

Сразу видно код «ручной работы», все так ровненько. Тут надо посмотреть, как он на кэш и конвейер ляжет.
Конкатенация матриц не сильно частая операция, в софтвер-рендере лучше оптимизировать другие части — например трансформацию вершин. Тут руками можно чудес наделать, несколько вершин (4?) за итерацию преобразовывать, да еще и в кеш префетчить на упреждение…
Умножение пары матриц 4х4 — это полная трансформация вершин четырехугольника.
А за что минус?

Я согласен с вами, однако по сравнению другими операциями — не такая уж и частая, согласитесь.
Я написал из своего опыта — что оптимизация конвеерных функций (тот же TnL), которые выполняют одни и теже операции, и часть данных константна — там можно получить сещественный прирост при правильной ручной оптимизации.
Что касается вопросов, связанных с обсуждением кодов — на мой взгляд, на профильном техническом ресурсе обнаружение и обсуждение скользких мест, недостатков, плохих практик и так далее должно вестись наиболее широко и гласно — это говорит о том, что аудитория ресурса откровенный ужас мимо себя просто так не пропустит.
Не стоит, согласен. Нужно поправить. Но этот пост не из хаба «идеальный код», и на него не претендует, поэтому можете смело сделать фикс и отправить пулл реквест, если хотите :)
В комментарии выше я предлагаю вариант исправления кошмарика.
Прошу простить мой резкий тон, был шокирован.
Для примеров лучше использовать код попроще и попонятнее — в три цикла.
Для реальности — нужно мерить, сравнивать и оптимизировать обоснованно.
И уж никак не предлагать такой код в качестве учебных примеров.

P.S. Это все IMHO, забейте — статья отличная! :-)
Можно использовать и готовую библиотеку матриц, например Eigen (там и SIMD встроено).
Немного не в тему, но хочу поделиться ссылкой на интересную реализацию программного рендера: www.ixbt.com/video/hier-tiling.html
Кстати программный рендер используется не только в академических целях. В ZBrush, например, используется полностью программный рендер, при этом количество полигонов при отрисовке может достигать нескольких десятков миллионов.
Простите пожалуйста, не могли бы поделиться источником инфы про софтверный рендер в ZBrush? Дико интересно.
Да собственно официальный сайт. Алгоритм по понятным причинам не доступен, но сам факт того, что там используется программный рендер написано на сайте:
ZBrush is software rendered, meaning that ZBrush itself is doing the rendering rather than the GPU. Your choice of GPU will not matter so long as it supports the recommended monitor resolution.

CPU: P4 or AMD Opteron or Athlon64 Processor (Must have SSE2: Streaming SIMD Extensions 2)

pixologic.com/zbrush/system/
Спасибо, казалось облазил весь сайт, а надо было в системные требования зайти.
<ИМХО> Software Render уже разобран и в статьях и в книгах.
Прошло даже несколько конкурсов на эту тему, например:
www.gamedev.ru/code/forum/?id=126749

Статья ничего нового не дает. </ИМХО>
Хороший учебник с примерами.

Всем начинающим пионерам рекомендую отличную книгу Роджерса, прочитав которую, вы станете великим в деле компьютерной графики до конца дней своих. Единственное, что в ней не хватает — это bumpmapping-а всякого и работы с текстурами.
Сюда же добавлю 2-томник Фоли и ван Дэма.
Большое спасибо за статью! Я как раз уже давно хотел написать софтрендерер, но с матчастью были проблемы (осилил пока только простой Wireframe-рендеринг без корректного отсечения), жду следующих частей статьи.
А потом переписать всё под OpenCL!
Довольно избитая тема. Любопытно будет посмотреть на реализацию растеризации, есть там подводные камни.
Подача материала в разных статьях и соответствующей литературе существенно различается. Всегда хорошо иметь выбор, что почитать на эту тему. Кому-то, возможно, приглянется мой вариант
Хорошо бы, если бы вы писали, какой у Вас FPS получается на каждом этапе, для сравнения с другими подобными реализациями на других языках и оценки оверхеда от каждого этапа конвеера (с указанием процессора, количества вершин в модели и размеров вьюпорта). К примеру, в wireframe какой получился?
Я вообще не занимался никакой оптимизацией на данном этапе. Может быть, позже, это интересная тема. В ближайшее время в статьях будет упор на описание алгоритмов, а не на их оптимизацию в коде. Сейчас, на модели с 27к вершин, 90 FPS (с -О3) с вьюпортом 640х480, i7-3770K
Алгебраично!
У меня есть такой глупый-глупый вопрос: а почему в 3д всегда используют трехточечную проекцию?

Или, иными словами, почему при малом угле наклона вертикальные линии не вертикальны? Я понимаю, что математика, но почему именно такая? Завалы вертикалей — это же плохо, не?
Может, я неверно понял вопрос, но разве кол-во точек схода в итоговом изображении не зависит от ориентации камеры? Если плоскость проекции пересекает все три оси, то мы получим трехточечную проекцию, но мы так же можем получить двух или одно-точечную, в зависимости от угла просмотра
Я к тому, что в рисовании принято избегать трехточечной проекции без серьёзных оснований. Нет ли таких систем рендеринга, которые бы выдерживали это правило? И более того, использовали бы более одного горизонта для помещений, и другие трюки из ассортимента художников. Вон, bloom уже используют…
Так можно легко этого добиться — держать ось вирт камеры всегда горизонтальной, а затем в постобработке кропнуть то что нужно.
Проблема изложена например в сто первом параграфе Ководства А. Лебедева.
Это, может быть, и плохо с «эстетической» точки зрения, зато абсолютно реалистично и соответствует действительно видимой картинке. К примеру, можете сфотографировать высокое здание, стоя не очень далеко от него — если вы поднимете камеру чуть-чуть вверх от горизонта при этом, увидите, что вертикальные стены здания сходятся. Это нормально.

Слегка неформальное объяснение: подъем камеры вверх — то же самое, что опускание всей сцены вниз, и вертикальные линии становятся, например, почти горизонтальными или даже строго горизонтальными (если угол равен 90 градусам). А параллельные горизонтальные (да и вообще какие угодно) линии становятся все ближе и ближе при приближении к горизонту из-за того, что расстояние между ними хоть и неизменно, но становится меньше в пересчете на угловой размер при удалении от наблюдателя.
Не соответствует оно никакой картинке. Проекция изображения в глазу идёт не на плоскость, а на вогнутую полусферу (неровной формы). Дальше на это накладываются алгоритмы компрессии, плюс алгоритм восстановления изображения, включающий в себя восстановление прямых.

Вот, простое вам упражнение: подойдите к столбу, и внимательно на него посмотрите, держа голову (глаза) прямо. Лучше с чем-то типа зубочистки для оценки относительных размеров. Столб будет толще в середине и тоньше к верху/низу (понятно почему — расстояние до центра столба меньше, чем до верха/низа). Таким образом вместо прямой линии будет кривая. И эти кривые — вокруг все. Но мозги их исправляют — и, заодно, исправляют мелкий угол от проекции, представляя прямоугольные объекты прямоугольными.

Ещё более важная вещь в рисовании: использование нескольких горизонтов. Это когда пол сходится к одному, а потолок к другому. Для человеческого глаза этот формат рисунка правильнее, чем формальное построение.
Не соглашусь. Неуместно сравнивать получаемую картинку с проекцией на глаз человека, нужно сравнивать с проекцией на экран, на который мы смотрим. Представьте, что вы смотрите на какую-то трехмерную сцену через «окно». Точки сцены проецируются на плоскость этого окна (= экрана). А уж и то, что получается на настоящем окне, и то, что получается на экране монитора, глаз человека спроецирует одинаково (то есть применительно к данной ситуации — вообще неважно, как именно, хоть вверх ногами в мозг).

Причина того, что отрендеренные картинки (и фотографии!) содержат то, что многие называют перспективными «искажениями» — несоответствие углу обзора. У современных камер горизонтальный угол обзора — что-то в районе 90 градусов, в играх тоже используются похожие значения. А смотрим мы на эту картинку на мониторе потом с такого расстояния, что ее угловой размер — градусов 40. Конечно, перспектива «плывет». В свою очередь, предложу вам свой эксперимент в ответ: возьмите картинку, FOV которой вы точно знаете, и расположитесь относительно нее так, чтобы FOV совпадал. В этом случае, говоря, опять же, «простым» языком, на все те вытянутые и «искаженные» детали ближе к краям картинки вы будете смотреть под таким (правильным) углом, что они сожмутся обратно в естественную с точки зрения вашего взгляда форму.

На ваш пример со столбом могу сказать, что уменьшение деталей с удалением от точки обзора (верх/низ столба) как раз и компенсируется тем, что на те «участки картинки», которые содержат эти точки столба, мы смотрим под таким углом, что, если отойти подальше и уменьшить свой угол обзора, они будут казаться растянутыми (см. предыдущий параграф).

Разумеется, это всё неприменимо к изогнутым экранам. Именно поэтому, в частности, в Oculus Rift картинка такая «выпученная», если ее перенести с рифтовского изогнутого экрана на плоский экран монитора.
В статье слишком много матана и по сути нет рендера.
Не так давно тоже баловался mentalx.org/tmp/metla_2014_09_22.7z 31к полигонов.
Рендер треугольника с освещением в SSE переписал, как оказалось, не зря, т.к. получил свои 10% к FPS по сравнению с SSE оптимизацией компилятора 8)
Суть рендерера — во многом «матан» :) Прежде чем писать, нужно понимать, что писать и откуда берутся такие формулы
Только полноправные пользователи могут оставлять комментарии.
Войдите, пожалуйста.