Обратная матрица. Метод Гаусса.

← →
@!!ex_ (2007-06-07 10:14) [0]

Я написал, но получилось слишком много кода и мне кажется, это медленно работает(если вообще работает, пока проверить не могу).
Плиз, посоветуйте, как сделать лучше.

Function Matrix.Inverse:boolean; var Ret:array[0..3,0..3] of single; Gauss:array[0..7,0..3] of single; Procedure Sub(Line1,Line2:integer; K:single); begin Gauss[0,Line1]:=Gauss[0,Line1]-Gauss[0,Line2]*K; Gauss[1,Line1]:=Gauss[1,Line1]-Gauss[1,Line2]*K; Gauss[2,Line1]:=Gauss[2,Line1]-Gauss[2,Line2]*K; Gauss[3,Line1]:=Gauss[3,Line1]-Gauss[3,Line2]*K; Gauss[4,Line1]:=Gauss[4,Line1]-Gauss[4,Line2]*K; Gauss[5,Line1]:=Gauss[5,Line1]-Gauss[5,Line2]*K; Gauss[6,Line1]:=Gauss[6,Line1]-Gauss[6,Line2]*K; Gauss[7,Line1]:=Gauss[7,Line1]-Gauss[7,Line2]*K; end; Procedure ChangeLines(Line1,Line2:integer); var K:single; begin K:=Gauss[0,Line1]; Gauss[0,Line1]:=Gauss[0,Line2]; Gauss[0,Line2]:=K; K:=Gauss[1,Line1]; Gauss[1,Line1]:=Gauss[1,Line2]; Gauss[1,Line2]:=K; K:=Gauss[2,Line1]; Gauss[2,Line1]:=Gauss[2,Line2]; Gauss[2,Line2]:=K; K:=Gauss[3,Line1]; Gauss[3,Line1]:=Gauss[3,Line2]; Gauss[3,Line2]:=K; K:=Gauss[4,Line1]; Gauss[4,Line1]:=Gauss[4,Line2]; Gauss[4,Line2]:=K; K:=Gauss[5,Line1]; Gauss[5,Line1]:=Gauss[5,Line2]; Gauss[5,Line2]:=K; K:=Gauss[6,Line1]; Gauss[6,Line1]:=Gauss[6,Line2]; Gauss[6,Line2]:=K; K:=Gauss[7,Line1]; Gauss[7,Line1]:=Gauss[7,Line2]; Gauss[7,Line2]:=K; end; Procedure Scale(Line:integer; K:single); begin Gauss[0,Line]:=Gauss[0,Line]*K; Gauss[1,Line]:=Gauss[1,Line]*K; Gauss[2,Line]:=Gauss[2,Line]*K; Gauss[3,Line]:=Gauss[3,Line]*K; Gauss[4,Line]:=Gauss[4,Line]*K; Gauss[5,Line]:=Gauss[5,Line]*K; Gauss[6,Line]:=Gauss[6,Line]*K; Gauss[7,Line]:=Gauss[7,Line]*K; end; begin Result:=False; Gauss[0,0]:=m_matrix[0]; Gauss[1,0]:=m_matrix[1]; Gauss[2,0]:=m_matrix[2]; Gauss[3,0]:=m_matrix[3]; Gauss[0,1]:=m_matrix[4]; Gauss[1,1]:=m_matrix[5]; Gauss[2,1]:=m_matrix[6]; Gauss[3,1]:=m_matrix[7]; Gauss[0,2]:=m_matrix[8]; Gauss[1,2]:=m_matrix[9]; Gauss[2,2]:=m_matrix[10]; Gauss[3,2]:=m_matrix[11]; Gauss[0,3]:=m_matrix[12]; Gauss[1,3]:=m_matrix[13]; Gauss[2,3]:=m_matrix[14]; Gauss[3,3]:=m_matrix[15]; Gauss[4,0]:=1; Gauss[5,0]:=0; Gauss[6,0]:=0; Gauss[7,0]:=0; Gauss[4,1]:=0; Gauss[5,1]:=1; Gauss[6,1]:=0; Gauss[7,1]:=0; Gauss[4,2]:=0; Gauss[5,2]:=0; Gauss[6,2]:=1; Gauss[7,2]:=0; Gauss[4,3]:=0; Gauss[5,3]:=0; Gauss[6,3]:=0; Gauss[7,3]:=1; //Gauss = M|E //Прямой ход //Первый элемент первой строки должен быть не 0, иначе не получить 1. if (Gauss[0,0]=0) and (Gauss[0,1]<>0) then ChangeLines(0,1) else if (Gauss[0,0]=0) and (Gauss[0,2]<>0) then ChangeLines(0,2) else if (Gauss[0,0]=0) and (Gauss[0,3]<>0) then ChangeLines(0,3) else Exit; Sub(1,0,Gauss[0,1]/Gauss[0,0]); Sub(2,0,Gauss[0,2]/Gauss[0,0]); Sub(3,0,Gauss[0,3]/Gauss[0,0]); if (Gauss[1,1]=0) and (Gauss[1,2]<>0) then ChangeLines(1,2) else if (Gauss[1,1]=0) and (Gauss[1,3]<>0) then ChangeLines(1,3) else Exit; Sub(2,1,Gauss[1,2]/Gauss[1,1]); Sub(3,1,Gauss[1,3]/Gauss[1,1]); if (Gauss[2,2]=0) and (Gauss[2,3]<>0) then ChangeLines(2,3) else Exit; Sub(3,2,Gauss[2,3]/Gauss[2,2]); //Получаем на диагонали 1 Scale(0,1/Gauss[0,0]); Scale(1,1/Gauss[1,1]); Scale(2,1/Gauss[2,2]); Scale(3,1/Gauss[3,3]); //Обратный ход Sub(2,0,Gauss[3,2]/Gauss[3,3]); Sub(1,0,Gauss[3,1]/Gauss[3,3]); Sub(0,0,Gauss[3,0]/Gauss[3,3]); Sub(1,0,Gauss[2,1]/Gauss[2,2]); Sub(0,0,Gauss[2,0]/Gauss[2,2]); Sub(0,0,Gauss[1,0]/Gauss[1,1]); Result:=true; end;

← →
MBo © (2007-06-07 10:18) [1]

http://alglib.sources.ru/linequations/obsolete/gaussm.php

← →
@!!ex_ (2007-06-07 10:31) [2]

> [1] MBo © (07.06.07 10:18)

Спасибо, конечно.
Там решение СЛУ, универсальное и в дебильном виде...
Смысл мне от него?
Мне надо обратную матрицу найти и не универсального размера, а конкретного. 4х4

← →
Думкин © (2007-06-07 10:34) [3]

> MBo © (07.06.07 10:18) [1]

Да он прикалывается.

← →
palva © (2007-06-07 10:46) [4]

Очень плохая идея проверять плавающее число на равенство.
Gauss[0,0]=0
Обычно ищут "главный элемент" в столбце, т. е. с максимальным модулем. И строку с главным элементом перемещают на первое место.

← →
palva © (2007-06-07 11:12) [5]

А зачем для частного случая 4x4 использовать алгоритм Гаусса? Есть много исследований на такие темы. Вот здесь посмотрите.
http://www.cellperformance.com/articles/2006/06/a_4x4_matrix_inverse_1.html

← →
@!!ex_ (2007-06-07 11:25) [6]

> http://www.cellperformance.com/articles/2006/06/a_4x4_matrix_inverse_1.ht
> ml

Не загрузилась статья.. Все оформление загрузилось, а статья - нет...

← →
Jeer © (2007-06-07 11:35) [7]

Значит так надо.

← →
palva © (2007-06-07 11:52) [8]

У меня грузится. Ну вот там основные формулы. Не знаю, понятно ли отобразится.

Inversion by Partitioning: To inverse a matrix A (size N) by partitioning, the matrix is partitioned into:
| A0 A1 | A = | | with A0 and A3 squared matrix with the respective size | A2 A3 | s0 and s3 following the rule: s0 + s3 = N
The inverse is
| B0 B1 | InvA = | | | B2 B3 |
with:

B0 = Inv(A0 - A1 * InvA3 * A2)
B1 = - B0 * (A1 * InvA3)
B2 = - (InvA3 * A2) * B0
B3 = InvA3 + B2 * (A1 * InvA3)

Здесь, правда, должно быть |A3|<>0
Если это не так, то нужно будет что-то делать переставляя строки или столбцы.

← →
@!!ex_ (2007-06-07 11:56) [9]

Так этож матрица 2х2 или я чего то не понимаю?

← →
palva © (2007-06-07 11:57) [10]

Если матрица 4x4 то разбить на клетки 2x2. Вычисления сведутся к обращению двух матриц 2x2 и матричной арифметике. Наверно подробности о клеточных матрицах можно найти в толстых книгах по матрицам, типа Гантмахера.

← →
Думкин © (2007-06-07 11:57) [11]

> @!!ex_ (07.06.07 11:56) [9]

Нет. Это не элементы - это блоки.

← →
@!!ex_ (2007-06-07 12:05) [12]

> [10] palva © (07.06.07 11:57)

> [11] Думкин © (07.06.07 11:57)

Мда... Не прочитал текст, сорри.

Че то не сказать, чтобы этот алгоритм был ыбстрее Гаусса.
Хотя надо сравнить будет.

← →
palva © (2007-06-07 12:13) [13]

> Че то не сказать, чтобы этот алгоритм был ыбстрее Гаусса.
Тут насколько я понял ищут не быстроту а возможность распараллеливания и использования специфических команд специфического процессора. Возможно там аппаратно реализована работа с матрицами 2x2. Но в самой статье я не разбирался.

← →
Jeer © (2007-06-07 12:46) [14]

В статье описано правило Крамера для инверсии матриц.
Эффективно для N<=6

Пример скорострельности для 4*4, cycles
Гаусс - 1074
Cramer rule - 846
Cramer rule with SIMD - 210

Так, что тут не думать надо, а реализовывать.

← →
Внук © (2007-06-07 15:09) [15]

Обращение матрицы 4x4 медленно работает? Чудеса...

← →
Sapersky (2007-06-07 15:32) [16]

Привет от страшного и ужасного GLScene (VectorGeometry.pas):

function MatrixDetInternal(const a1, a2, a3, b1, b2, b3, c1, c2, c3: Single): Single;
// internal version for the determinant of a 3x3 matrix
begin
Result:= a1 * (b2 * c3 - b3 * c2)
- b1 * (a2 * c3 - a3 * c2)
+ c1 * (a2 * b3 - a3 * b2);
end;

function MatrixDeterminant(const M: TMatrix): Single;
begin
Result:= M[X, X]*MatrixDetInternal(M[Y, Y], M[Z, Y], M[W, Y], M[Y, Z], M[Z, Z], M[W, Z], M[Y, W], M[Z, W], M[W, W])
-M[X, Y]*MatrixDetInternal(M[Y, X], M[Z, X], M[W, X], M[Y, Z], M[Z, Z], M[W, Z], M[Y, W], M[Z, W], M[W, W])
+M[X, Z]*MatrixDetInternal(M[Y, X], M[Z, X], M[W, X], M[Y, Y], M[Z, Y], M[W, Y], M[Y, W], M[Z, W], M[W, W])
-M[X, W]*MatrixDetInternal(M[Y, X], M[Z, X], M[W, X], M[Y, Y], M[Z, Y], M[W, Y], M[Y, Z], M[Z, Z], M[W, Z]);
end;

procedure AdjointMatrix(var M : TMatrix);
var
a1, a2, a3, a4,
b1, b2, b3, b4,
c1, c2, c3, c4,
d1, d2, d3, d4: Single;
begin
a1:= M[X, X]; b1:= M[X, Y];
c1:= M[X, Z]; d1:= M[X, W];
a2:= M[Y, X]; b2:= M[Y, Y];
c2:= M[Y, Z]; d2:= M[Y, W];
a3:= M[Z, X]; b3:= M[Z, Y];
c3:= M[Z, Z]; d3:= M[Z, W];
a4:= M[W, X]; b4:= M[W, Y];
c4:= M[W, Z]; d4:= M[W, W];

// row column labeling reversed since we transpose rows & columns
M[X, X]:= MatrixDetInternal(b2, b3, b4, c2, c3, c4, d2, d3, d4);
M[Y, X]:=-MatrixDetInternal(a2, a3, a4, c2, c3, c4, d2, d3, d4);
M[Z, X]:= MatrixDetInternal(a2, a3, a4, b2, b3, b4, d2, d3, d4);
M[W, X]:=-MatrixDetInternal(a2, a3, a4, b2, b3, b4, c2, c3, c4);

M[X, Y]:=-MatrixDetInternal(b1, b3, b4, c1, c3, c4, d1, d3, d4);
M[Y, Y]:= MatrixDetInternal(a1, a3, a4, c1, c3, c4, d1, d3, d4);
M[Z, Y]:=-MatrixDetInternal(a1, a3, a4, b1, b3, b4, d1, d3, d4);
M[W, Y]:= MatrixDetInternal(a1, a3, a4, b1, b3, b4, c1, c3, c4);

M[X, Z]:= MatrixDetInternal(b1, b2, b4, c1, c2, c4, d1, d2, d4);
M[Y, Z]:=-MatrixDetInternal(a1, a2, a4, c1, c2, c4, d1, d2, d4);
M[Z, Z]:= MatrixDetInternal(a1, a2, a4, b1, b2, b4, d1, d2, d4);
M[W, Z]:=-MatrixDetInternal(a1, a2, a4, b1, b2, b4, c1, c2, c4);

M[X, W]:=-MatrixDetInternal(b1, b2, b3, c1, c2, c3, d1, d2, d3);
M[Y, W]:= MatrixDetInternal(a1, a2, a3, c1, c2, c3, d1, d2, d3);
M[Z, W]:=-MatrixDetInternal(a1, a2, a3, b1, b2, b3, d1, d2, d3);
M[W, W]:= MatrixDetInternal(a1, a2, a3, b1, b2, b3, c1, c2, c3);
end;

procedure ScaleMatrix(var M : TMatrix; const factor : Single);
var
i : Integer;
begin
for i:=0 to 3 do begin
M[I, 0]:=M[I, 0] * Factor;
M[I, 1]:=M[I, 1] * Factor;
M[I, 2]:=M[I, 2] * Factor;
M[I, 3]:=M[I, 3] * Factor;
end;
end;

Const
IdentityHmgMatrix: TMatrix = ((1, 0, 0, 0),
(0, 1, 0, 0),
(0, 0, 1, 0),
(0, 0, 0, 1));

procedure InvertMatrix(var M : TMatrix);
var
det : Single;
begin
det:=MatrixDeterminant(M);
if Abs(Det)<EPSILON then
M:=IdentityHmgMatrix
else begin
AdjointMatrix(M);
ScaleMatrix(M, 1/det);
end;
end;

← →
@!!ex_ (2007-06-07 15:42) [17]

> [16] Sapersky (07.06.07 15:32)

Удажсно медленно.
Нас в ВУЗе учили не пользоваться этим методом в порграммировании, посольку он жутко медленный, и Гаусс будет работать быстрее, поэтому я Гаусс и выбрал.

← →
Jeer © (2007-06-07 15:49) [18]

> методом в пор..граммировании

Чему только сейчас не учат:))
Всему ? :)))

← →
MBo © (2007-06-07 16:11) [19]

>Мне надо обратную матрицу найти и не универсального размера, а конкретного. 4х4

Может, у тебя и матрица не общего вида, а, скажес, аффинного преобразования?

← →
Sapersky (2007-06-07 16:36) [20]

Нашёл у себя более простой вариант. Даже подозрительно простой, но по крайней мере с видовой матрицей он работал (перевод координат курсора мыши в мировые):

function MatrixInvert(Const a: TMatrix; Var q: TMatrix): Boolean;
var
DetInv: Single;
begin
Result := False;
if (abs(a._44 - 1) > 0.001) or
(abs(a._14) > 0.001) or
(abs(a._24) > 0.001) or
(abs(a._34) > 0.001) then Exit;
DetInv := 1 /( a._11 * ( a._22 * a._33 - a._23 * a._32 ) -
a._12 * ( a._21 * a._33 - a._23 * a._31 ) +
a._13 * ( a._21 * a._32 - a._22 * a._31 ) );

q._11 := DetInv * ( a._22 * a._33 - a._23 * a._32 );
q._12 := -DetInv * ( a._12 * a._33 - a._13 * a._32 );
q._13 := DetInv * ( a._12 * a._23 - a._13 * a._22 );
q._14 := 0;

q._21 := -DetInv * ( a._21 * a._33 - a._23 * a._31 );
q._22 := DetInv * ( a._11 * a._33 - a._13 * a._31 );
q._23 := -DetInv * ( a._11 * a._23 - a._13 * a._21 );
q._24 := 0;

q._31 := DetInv * ( a._21 * a._32 - a._22 * a._31 );
q._32 := -DetInv * ( a._11 * a._32 - a._12 * a._31 );
q._33 := DetInv * ( a._11 * a._22 - a._12 * a._21 );
q._34 := 0;

q._41 := -( a._41 * q._11 + a._42 * q._21 + a._43 * q._31 );
q._42 := -( a._41 * q._12 + a._42 * q._22 + a._43 * q._32 );
q._43 := -( a._41 * q._13 + a._42 * q._23 + a._43 * q._33 );
q._44 := 1;

Result := True;
end;

← →
palva © (2007-06-07 17:18) [21]

http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm

← →
Правильный Вася (2007-06-07 17:20) [22]

> Обратная матрица. Метод Гаусса.
ты гоняешься за агентами по главной диагонали?

← →
Jeer © (2007-06-07 18:00) [23]

http://www.cvl.iis.u-tokyo.ac.jp/~miyazaki/tech/teche23.html

← →
palva © (2007-06-07 18:29) [24]

Jeer © (07.06.07 18:00) [23]
Автора по ссылке интересно зовут. Дайсуке Миязаки

← →
Jeer © (2007-06-07 18:51) [25]

> palva © (07.06.07 18:29) [24]

Я тоже повеселился:)

← →
palva © (2007-06-08 11:40) [26]

> palva © (07.06.07 11:52) [8]
> B3 = InvA3 + B2 * (A1 * InvA3)

Вынужден извиниться. По ссылке в формуле ошибка. Нужен минус вместо плюса.
B3 = InvA3 - B2 * (A1 * InvA3)
Сейчас попробую проверить численно.

← →
palva © (2007-06-08 13:21) [27]

Ну вот написал.
Я не стал бы писать этот код,
1) если бы по моей ссылке не стояла ошибочная формула и не требовалось бы загладить мою вину перед автором;
2) если бы этот код не послужил примером программы, которую нельзя (по крайней мере я не умею) написать на паскале эффективно и прозрачно без большого числа комментариев, в то время как на си с использованием препроцессора подобный код был бы прозрачным и комментариев не потребовал. Правда, паскалевский код возможно пропустить через сишный препроцессор, но это уже для очень больших поклонников паскаля.

Я старался писать код максимальный по скорости. Кстати, если уважаемый @!!ex_ перепишет этот код на си, хотя бы только одну главную процедуру, то он получит еще небольшое увеличение производительности.
{$APPTYPE CONSOLE} uses Math; type TMat4 = array[0..3, 0..3] of Double; procedure Inv4(var A, B: TMat4); // type // TMat2 = array[0..1, 0..1] of Double; procedure Inv2(var M1, M2); var Det: Double; begin Det := TMat4(M1)[0,0]*TMat4(M1)[1,1] - TMat4(M1)[1,0]*TMat4(M1)[0,1]; TMat4(M2)[0,0] := TMat4(M1)[1,1]/Det; TMat4(M2)[0,1] := -TMat4(M1)[0,1]/Det; TMat4(M2)[1,0] := -TMat4(M1)[1,0]/Det; TMat4(M2)[1,1] := TMat4(M1)[0,0]/Det; end; procedure Mul(var M1, M2, M3); begin TMat4(M3)[0,0] := TMat4(M1)[0,0]*TMat4(M2)[0,0] + TMat4(M1)[0,1]*TMat4(M2)[1,0]; TMat4(M3)[0,1] := TMat4(M1)[0,0]*TMat4(M2)[0,1] + TMat4(M1)[0,1]*TMat4(M2)[1,1]; TMat4(M3)[1,0] := TMat4(M1)[1,0]*TMat4(M2)[0,0] + TMat4(M1)[1,1]*TMat4(M2)[1,0]; TMat4(M3)[1,1] := TMat4(M1)[1,0]*TMat4(M2)[0,1] + TMat4(M1)[1,1]*TMat4(M2)[1,1]; end; procedure Neg(var M1); begin TMat4(M1)[0,0] := -TMat4(M1)[0,0]; TMat4(M1)[0,1] := -TMat4(M1)[0,1]; TMat4(M1)[1,0] := -TMat4(M1)[1,0]; TMat4(M1)[1,1] := -TMat4(M1)[1,1]; end; procedure SubNeg(var M1, M2); begin TMat4(M2)[0,0] := TMat4(M1)[0,0] - TMat4(M2)[0,0]; TMat4(M2)[0,1] := TMat4(M1)[0,1] - TMat4(M2)[0,1]; TMat4(M2)[1,0] := TMat4(M1)[1,0] - TMat4(M2)[1,0]; TMat4(M2)[1,1] := TMat4(M1)[1,1] - TMat4(M2)[1,1]; end; var T: TMat4; { Здесь бы очень пригодился сишный препроцессор тогда было бы возможно написать прозрачный и эффективный код практически без комментариев. #define A0 A[0,0] #define A1 A[0,2] #define A2 A[2,0] #define A3 A[2,2] #define B0 B[0,0] #define B1 B[0,2] #define B2 B[2,0] #define B3 B[2,2] #define InvA3 T[0,0] #define Temp T[0,2] #define InvA3A2 T[2,0] #define A1InvA3 T[2,2] } begin // Inv2(A3, InvA3); // Mul(InvA3, A2, InvA3A2); // Mul(A1, InvA3, A1InvA3); Inv2(A[2,2], T[0,0]); Mul(T[0,0], A[2,0], T[2,0]); Mul(A[0,2], T[0,0], T[2,2]); // B0 = Inv(A0 - A1 * (InvA3 * A2)) // Mul(A1, InvA3A2, Temp); // SubNeg(A0, Temp); // Inv2(Temp, B0); Mul(A[0,2], T[2,0], T[0,2]); SubNeg(A[0,0], T[0,2]); Inv2(T[0,2], B[0,0]); // B1 = - B0 * (A1 * InvA3) // Mul(B0, A1InvA3, B1); // Neg(B1); Mul(B[0,0], T[2,2], B[0,2]); Neg(B[0,2]); // B2 = - (InvA3 * A2) * B0 // Mul(InvA3A2, B0, B2); // Neg(B2); Mul(T[2,0], B[0,0], B[2,0]); Neg(B[2,0]); // B3 = InvA3 - B2 * (A1 * InvA3) // Mul(B2, A1InvA3, B3); // SubNeg(InvA3, B3); Mul(B[2,0], T[2,2], B[2,2]); SubNeg(T[0,0], B[2,2]); end; var AA, BB, CC: TMat4; i, j, k: Integer; begin for i := 0 to 3 do for j:= 0 to 3 do AA[i,j] := 10.0 * random - 5.0; Inv4(AA, BB); 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 © (2007-06-09 12:04) [28]

> palva © (08.06.07 13:21) [27]

Что интересно, реализация по ссылке [23] дает не намного худшие результаты, хотя количество умножений явно больше.

[23] - 3.0 kticks
[27] - 2.5 kticks

Вообще, при "ковырянии" в такой, казалось бы простой реализации сделал для себя несколько маленьких "открытий", применительно, по всей видимости, к особенностям выполнения программ на современных процессорах.

1.Вынос общих членов из тройных произведений порой приводит к ухудшению результатов, т.е. a*b*c + a*d*e = a*(b*c + d*e) - экономит одно умножение, но выполняется медленнее.

2. Использование Result: double в качестве временной переменной ухудшает время выполнения, по ср. с явно выделенной переменной x: double и в конце
Result := x;

3. Прибавку в скорости дает комбинирование обращений к членам массива
путем предварительного присвоения их временным переменным, а к некоторым
членам массива (по первым столбцам) стоит обращаться явно.

4. Улучшает дело и техника обращения к двумерным массивам, как к одномерным.

В качестве иллюстрации к 4 приведена реализация обращения матрицы 4х4
посредством иного разбиения на подматрицы 3-го порядка и работа в циклах.
Разумеется, скорострельность значительно ниже ( 14 kticks) за счет накладных
расходов на циклы и индексации, но в качестве иллюстрации сгодится.

{$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 GetCPUTicks_: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[dii * 3 + 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 );
TV4(mr)[i + j * 4] := mdet * ( m3_det( mtemp ) * sgn );
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 := GetCPUTicks_;
m4_inverse(AA, BB);
T2 := T1 - GetCPUTicks_;

Writeln(T2," 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.

Когда-то нас учили, что умножение выполняется на порядок медленнее сложения, а деление еще раз в 10 медленнее. Соответственно этому вся теория оценки сложности таких алгоритмов строилась на подсчете числа умножений. А уж несколько делений на одно и то же число всегда заменялось предварительным однократным вычислением обратной величины и дальнейшей заменой деления умножением. Нынешние сопроцессоры, похоже, устроены так, что время любой плавающей операции примерно одно и то же. На эту тему я, к сожалению, могу лишь задавать вопросы.

> palva © (09.06.07 14:45) [29]

Деление и сегодня медленней, чем умножение.
В Вашем примере [27] замена на обратную величину и умножение приводит к результату
было 2500 ticks
стало 2100 ticks

И видимо, деление заметно медленнее, раз такой эффект. Там ведь всего 8 делений!

Да, всего лишь замена на:

procedure Inv2(var M1, M2); var Det: Double; begin Det := 1.0/(TMat4(M1)[0,0]*TMat4(M1)[1,1] - TMat4(M1)[1,0]*TMat4(M1)[0,1]); TMat4(M2)[0,0] := TMat4(M1)[1,1]*Det; TMat4(M2)[0,1] := -TMat4(M1)[0,1]*Det; TMat4(M2)[1,0] := -TMat4(M1)[1,0]*Det; TMat4(M2)[1,1] := TMat4(M1)[0,0]*Det; 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; чем Result := Result + mr[n] * det * i;

← →
Sapersky (2007-06-09 16:03) [35]

Может, у тебя и матрица не общего вида, а, скажес, аффинного преобразования?

Ага:
http://delphimaster.net/view/9-1181032876/
И вариант из [20] именно для таких матриц (последний столбец 0,0,0,1).
В дополнение к нему:
Type
PMatrix = ^TMatrix;
TMatrix = record
Case integer of
0 : (_11, _12, _13, _14: Single;
_21, _22, _23, _24: Single;
_31, _32, _33, _34: Single;
_41, _42, _43, _44: Single);
1 : (m : array [0..3, 0..3] of Single);
end;
(Single для 3D-графики достаточно).

Ещё несколько ускоряет процесс установка соответствующей точности для FPU:

Type
TFPUPrecision = (fpDefault, fpSingle, fpDouble, fpExtended);

procedure SetPrecision(fp : TFPUPrecision = fpDefault);
Var cw : Word;
begin
Case fp of
fpDefault : cw := Default8087CW;
fpSingle : cw := Default8087CW and $FCFF;
fpDouble : cw := (Default8087CW and $FCFF) or $0200;
fpExtended : cw := Default8087CW or $0300;
// no exceptions: $133
end;
Set8087CW(cw);
end;

Другое дело, что в случае @!!ex_ особо оптимизировать обращение матрицы смысла не имеет, всё равно время построения теневого объёма будет на порядок больше (в матрице 16 чисел, а в модели несколько тысяч вершин).

Для меня было откровением, что это быстрее

FP-asm я почти не знаю, но с целочисленным алгоритмом была схожая ситуация, когда неочевидные изменения кода приводили к перемещению переменных из стека в регистры и довольно значительному ускорению.

Вот, кстати, один товарищ пишет, что Дельфи вообще толком не умеет оптимизировать floating-point вычисления:

http://dennishomepage.gugs-cats.dk/CodingForSpeedInDelphi.doc

Sapersky (09.06.07 16:03) [35]
> Вот, кстати, один товарищ пишет, что Дельфи вообще толком не умеет оптимизировать floating-point вычисления:

Я попробовал сишный вариант этой же программы (миллион обращений, поправка Jeer © (09.06.07 15:25) [33] учтена).
Borland С++ 5.5 дал 312 тиков, Delphi 7 - 281 тик.
Все компиляции производились с ключами по умолчанию.

> Jeer © (09.06.07 16:48) [37]
> Интересно, почему такая разница

Да странно. Вот сишный модуль. Переписал один в один.
void Inv4(double A[4][4], double B[4][4]); static void Inv2(double M1[4][4], double M2[4][4]) { double Det; Det = 1.0/(M1[0][0]*M1[1][1] - M1[1][0]*M1[0][1]); M2[0][0] = M1[1][1]*Det; M2[0][1] = -M1[0][1]*Det; M2[1][0] = -M1[1][0]*Det; M2[1][1] = M1[0][0]*Det; }; static void Mul(double M1[4][4], double M2[4][4], double M3[4][4]) { M3[0][0] = M1[0][0]*M2[0][0] + M1[0][1]*M2[1][0]; M3[0][1] = M1[0][0]*M2[0][1] + M1[0][1]*M2[1][1]; M3[1][0] = M1[1][0]*M2[0][0] + M1[1][1]*M2[1][0]; M3[1][1] = M1[1][0]*M2[0][1] + M1[1][1]*M2[1][1]; }; static void Neg(double M1[4][4]) { M1[0][0] = -M1[0][0]; M1[0][1] = -M1[0][1]; M1[1][0] = -M1[1][0]; M1[1][1] = -M1[1][1]; }; static void SubNeg(double M1[4][4], double M2[4][4]) { M2[0][0] = M1[0][0] - M2[0][0]; M2[0][1] = M1[0][1] - M2[0][1]; M2[1][0] = M1[1][0] - M2[1][0]; M2[1][1] = M1[1][1] - M2[1][1]; }; void Inv4(double A[4][4], double B[4][4]) { double T[4][4]; #define A0 (double (*)[4])&A[0][0] #define A1 (double (*)[4])&A[0][2] #define A2 (double (*)[4])&A[2][0] #define A3 (double (*)[4])&A[2][2] #define B0 (double (*)[4])&B[0][0] #define B1 (double (*)[4])&B[0][2] #define B2 (double (*)[4])&B[2][0] #define B3 (double (*)[4])&B[2][2] #define InvA3 (double (*)[4])&T[0][0] #define Temp (double (*)[4])&T[0][2] #define InvA3A2 (double (*)[4])&T[2][0] #define A1InvA3 (double (*)[4])&T[2][2] Inv2(A3, InvA3); Mul(InvA3, A2, InvA3A2); Mul(A1, InvA3, A1InvA3); // B0 = Inv(A0 - A1 * (InvA3 * A2)) Mul(A1, InvA3A2, Temp); SubNeg(A0, Temp); Inv2(Temp, B0); // B1 = - B0 * (A1 * InvA3) Mul(B0, A1InvA3, B1); Neg(B1); // B2 = - (InvA3 * A2) * B0 Mul(InvA3A2, B0, B2); Neg(B2); // B3 = InvA3 - B2 * (A1 * InvA3) Mul(B2, A1InvA3, B3); SubNeg(InvA3, B3); }

procedure Inv2(var M1, M2); передаётся указатель

static void Inv2(double M1[4][4], double M2[4][4]) передаётся весь массив вроде как. + static вроде не очень хорошо.

> static void Inv2(double M1[4][4], double M2[4][4]) передаётся весь массив вроде как
Вроде как нет:
; Inv2(A3, InvA3); ; ?live16385@16: ; EBX = B, ESI = A, EDI = &T @10: push edi lea eax,dword ptr [esi+80] push eax call @Inv2$qpa4$dt1 add esp,8

Обратная матрица. Метод Гаусса. Найти похожие ветки