Быстрое инвертирование матрицы 4х4

← →
@!!ex © (2008-06-27 10:33) [0]

Как?
Посомтрел тут свой код...
Инвертирование матрицы методом Гаусса... Больше 150 строк кода... Вообще не радует....
Возможно есть какое то частное решение, подходящее только для матрицы 4х4, о котором я не знаю?
Как например расчет определителя для матрицы 3х3 методом треугольника.

← →
han_malign © (2008-06-27 10:47) [1]

http://algolist.manual.ru/maths/linalg/index.php

← →
Jeer © (2008-06-27 10:53) [2]

Мы тут как-то проезжали эту тему по поводу быстродействия.
Вот один из моих вариантов на 2800 ticks:

{$APPTYPE CONSOLE}
program inv;

uses Math, TimeUnit;

type
TMat4 = array[0..3, 0..3] of double;
TV4 = array[0..15] of double;
var
AA, BB, CC: TMat4;
i, j, k: Integer;
T1,T2: int64;

procedure Inv4(const S; var T);
var det: double;
A12,A13,A14: double;
A22,A23,A24: double;
A32,A33,A34: double;
A41,A42,A43,A44: double;
B1,B2,B3,B4,B5,B6,B7,B8: double;
begin

A12 := TV4(S)[1];
A13 := TV4(S)[2];
A14 := TV4(S)[3];
//
A22 := TV4(S)[5];
A23 := TV4(S)[6];
A24 := TV4(S)[7];
//
A32 := TV4(S)[9];
A33 := TV4(S)[10];
A34 := TV4(S)[11];
//
A41 := TV4(S)[12];
A42 := TV4(S)[13];
A43 := TV4(S)[14];
A44 := TV4(S)[15];

B1 := A22*A33*A44 + A23*A34*A42 + A24*A32*A43;
B2 := TV4(S)[4]*A34*A43 + A23*TV4(S)[8]*A44 + A24*A33*A41;
B3 := TV4(S)[4]*A32*A44 + A22*A34*A41 + A24*TV4(S)[8]*A42;
B4 := TV4(S)[4]*A33*A42 + A22*TV4(S)[8]*A43 + A23*A32*A41;

B5 := A22*A34*A43 + A23*A32*A44 + A24*A33*A42;
B6 := TV4(S)[4]*A33*A44 + A23*A34*A41 + A24*TV4(S)[8]*A43;
B7 := TV4(S)[4]*A34*A42 + A22*TV4(S)[8]*A44 + A24*A32*A41;
B8 := TV4(S)[4]*A32*A43 + A22*A33*A41 + A23*TV4(S)[8]*A42;

det := 1.0 / (TV4(S)[0]*(B1-B5) + A12*(B2-B6) + A13*(B3-B7) + A14*(B4-B8));

TV4(T)[0] := det*(B1 - B5);
TV4(T)[1] := det*(A12*(A34*A43 - A33*A44) + A13*(A32*A44 - A34*A42) + A14*(A33*A42 - A32*A43));
TV4(T)[2] := det*(A12*A23*A44 + A13*A24*A42 + A14*A22*A43 - A12*A24*A43 - A13*A22*A44 - A14*A23*A42);
TV4(T)[3] := det*(A12*A24*A33 + A13*A22*A34 + A14*A23*A32 - A12*A23*A34 - A13*A24*A32 - A14*A22*A33);

TV4(T)[4] := det*(B2 - B6);
TV4(T)[5] := det*(TV4(S)[0]*(A33*A44 - A34*A43) + A13*(A34*A41 - TV4(S)[8]*A44) + A14*(TV4(S)[8]*A43 - A33*A41));
TV4(T)[6] := det*(TV4(S)[0]*(A24*A43 - A23*A44) + A13*(TV4(S)[4]*A44 - A24*A41) + A14*(A23*A41 - TV4(S)[4]*A43));
TV4(T)[7] := det*(TV4(S)[0]*(A23*A34 - A24*A33) + A13*(A24*TV4(S)[8] - TV4(S)[4]*A34) + A14*(TV4(S)[4]*A33 - A23*TV4(S)[8]));

TV4(T)[8] := det*(B3 - B7);
TV4(T)[9] := det*(TV4(S)[0]*(A34*A42 - A32*A44) + A12*(TV4(S)[8]*A44 - A34*A41) + A14*(A32*A41 - TV4(S)[8]*A42));
TV4(T)[10] := det*(TV4(S)[0]*(A22*A44 - A24*A42) + A12*(A24*A41 - TV4(S)[4]*A44) + A14*(TV4(S)[4]*A42 - A22*A41));
TV4(T)[11] := det*(TV4(S)[0]*(A24*A32 - A22*A34) + A12*(TV4(S)[4]*A34 - A24*TV4(S)[8]) + A14*(A22*TV4(S)[8] - TV4(S)[4]*A32));

TV4(T)[12] := det*(B4 - B8);
TV4(T)[13] := det*(TV4(S)[0]*(A32*A43 - A33*A42) + A12*(A33*A41 - TV4(S)[8]*A43) + A13*(TV4(S)[8]*A42 - A32*A41));
TV4(T)[14] := det*(TV4(S)[0]*(A23*A42 - A22*A43) + A12*(TV4(S)[4]*A43 - A23*A41) + A13*(A22*A41 - TV4(S)[4]*A42));
TV4(T)[15] := det*(TV4(S)[0]*(A22*A33 - A23*A32) + A12*(A23*TV4(S)[8] - TV4(S)[4]*A33) + A13*(TV4(S)[4]*A32 - A22*TV4(S)[8]));

end;

begin

for i := 0 to 3 do
for j := 0 to 3 do
AA[i,j] := 10.0 * random - 5.0;

T1 := GetCPUTicksP_;
Inv4(AA, BB);
T2 := GetCPUTicksP_;

Writeln((T2-T1)," ticks");
WriteLn;

for i := 0 to 3 do begin
for j := 0 to 3 do begin
CC[i,j] := 0.0;
for k := 0 to 3 do
CC[i,j] := CC[i,j] + AA[i,k] * BB[k,j];
Write(CC[i,j]:12:6);
end;
WriteLn
end;
// Должна вывестись единичная матрица
readln;
end.

← →
Jeer © (2008-06-27 11:03) [3]

А это вариант с вычислением подматриц 3*3, но хуже по быстродействию.

{$APPTYPE CONSOLE}
program inv;

uses Math;

type
TMat4 = array[0..3, 0..3] of Double;
TV3 = array[0..8] of Double;
TV4 = array[0..15] of Double;

function GetCPUTicksP_:int64;
// 88 ticks
asm
dw 310Fh // rdtsc
end;

procedure m4_submat(const mr: TV4; var mb: TV3; i, j: integer );
var
di, dj, si, sj: integer;
begin
si := 0; sj := 0;
// loop through 3x3 submatrix
for di := 0 to 2 do
for dj := 0 to 2 do begin
// map 3x3 element (destination) to 4x4 element (source)
if di >= i then si := di + 1
else si := di;
if dj >= j then sj := dj + 1
else sj := dj;
// copy element
mb[di shl 1 + di + dj] := mr[si shl 2 + sj];
end;
end;

function m3_det( const mat: TV3 ): double;
begin
Result := mat[0] * mat[4] * mat[8] - mat[0] * mat[7] * mat[5]
- mat[1] * mat[3] * mat[8] + mat[1] * mat[6] * mat[5]
+ mat[2] * mat[3] * mat[7] - mat[2] * mat[6] * mat[4];
end;

function m4_det( const mr: TV4 ): double;
var
det: double;
i,n: integer;
msub3: TV3;
x: double;
begin
x := 0.0;
i := 1;
for n := 0 to 3 do begin
i := -i;
m4_submat( mr, msub3, 0, n );
det := m3_det( msub3 );
x := x + mr[n] * det * i;
end;
Result := x;
end;

procedure m4_identity(var mr: TV4);
var i: integer;
begin
for i:=0 to 3 do
mr[i*5] := 1.0;
end;

function m4_inverse( var ma, mr ): integer;
var
mdet: double;
mtemp: TV3;
i, j, sgn: integer;
begin

mdet := 1.0 / m4_det( TV4(ma) );
if ( abs( mdet ) < 0.0005 ) then begin
m4_identity( TV4(mr) );
Result := 0;
Exit;
end;
for i := 0 to 3 do
for j := 0 to 3 do begin
sgn := 1 - ( (i + j) mod 2 ) * 2;
m4_submat( TV4(ma), mtemp, i, j );
if sgn >=0 then
TV4(mr)[i + j * 4] := mdet * m3_det( mtemp )
else
TV4(mr)[i + j * 4] := - mdet * m3_det( mtemp )

end;
Result := 1;
end;

var
AA, BB, CC: TMat4;
i, j, k: integer;
T1,T2: int64;
begin

for i := 0 to 3 do
for j:= 0 to 3 do begin
AA[i,j] := 10.0 * random - 5.0;
end;

T1 := GetCPUTicksP_;
m4_inverse(AA, BB);
T2 := GetCPUTicksP_;

Writeln((T2-T1)," ticks");
WriteLn;

for i := 0 to 3 do begin
for j := 0 to 3 do begin
CC[i,j] := 0.0;
for k := 0 to 3 do
CC[i,j] := CC[i,j] + AA[i,k] * BB[k,j];
Write(CC[i,j]:12:6);
end;
WriteLn;
end;
// Должна вывестись единичная матрица
readln;
end.

P.S.

У Intel в апноте AP-928 за 1999 г. есть примеры использования SIMD для инверсии матрицы 4*4.
С использованием правила Крамера там вообще вышло 210 cycles :)

← →
Sapersky (2008-06-27 12:12) [4]

Мы тут как-то проезжали эту тему по поводу быстродействия.

Причём @!!ex же и спрашивал, примерно год назад.
Вот эта ветка (не полностью, в архиве лень ковыряться):
http://narod.ru/disk/1147837000/InverseMatrix.zip

← →
@!!ex © (2008-06-27 12:47) [5]

> [4] Sapersky (27.06.08 12:12)

Какая память однако. :) Я уже и забыл что спрашивал это.

← →
@!!ex © (2008-06-27 12:49) [6]

Вот есть такой вариант:
http://www.everfall.com/paste/id.php?z0neuzcvr9fy

Насколько это в тему?

P.S.
Да, матрица не общая, матрица хранит только вращение.

← →
@!!ex © (2008-06-27 13:02) [7]

Это в тему. Работает.

← →
Zeqfreed © (2008-06-27 13:11) [8]

> Да, матрица не общая, матрица хранит только вращение.

На википедии написано:
In matrix theory, a rotation matrix is a real square matrix whose transpose is its inverse and whose determinant is +1.

http://en.wikipedia.org/wiki/Rotation_matrix

← →
Sapersky (2008-06-27 14:27) [9]

http://www.everfall.com/paste/id.php?z0neuzcvr9fy

Похоже на вариант, который я давал ещё в той ветке, но математически сокращённый (корректность преобразований не проверял). Хотя в моём несколько оптимальнее считалась последняя строка матрицы, с учётом m->_11*l1 = res->_11 и т.д.

> [9] Sapersky (27.06.08 14:27)

Раз уж ты все равно здесь, может поможешь разобраться, почему SV не корректно строится?
http://img80.imageshack.us/img80/6520/swbugvz6.jpg

← →
Sapersky (2008-06-27 17:13) [11]

Сходу не скажу - 3D-графикой давно не занимался, да и год назад только эпизодически поковырялся с SV, так что о типовых глюках имею слабое представление.
Не знаю, имеет ли смысл отсылать сюда:
http://delphimaster.net/view/9-1181032876/
хотя судя по начальному вопросу ветки, смысл вполне может быть :)
Во всяком случае, там пример гарантированно работающего кода, правда, он относится к z-pass, для z-fail нужны ещё передняя/задняя "затычки".

← →
Sapersky (2008-06-27 18:07) [12]

Могу предположить, что моделька кривовата, некоторые треугольники - двухсторонние, и рёбра их "задних" сторон воспринимаются как часть силуэта (алгоритм выделения "силуэтных" ребёр определяет их по одинарному вхождению в список, и для одиночных "задних" cторон это условие соблюдается).

У меня код из http://delphimaster.net/view/9-1181032876/ не пахал нормально.
ПРишлось дублировать грани.
А проблема на скриншоте тупо в том, что я нормаль не правильно считал... :)
Когда писал перемножение векторов, в одном месте вместо x - y написал... и все поломалось.

Быстрое инвертирование матрицы 4х4 Найти похожие ветки