Форум: "Прочее";
Текущий архив: 2007.07.15;
Скачать: [xml.tar.bz2];
ВнизОбратная матрица. Метод Гаусса. Найти похожие ветки
← →
@!!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.
← →
palva © (2007-06-09 14:45) [29]Когда-то нас учили, что умножение выполняется на порядок медленнее сложения, а деление еще раз в 10 медленнее. Соответственно этому вся теория оценки сложности таких алгоритмов строилась на подсчете числа умножений. А уж несколько делений на одно и то же число всегда заменялось предварительным однократным вычислением обратной величины и дальнейшей заменой деления умножением. Нынешние сопроцессоры, похоже, устроены так, что время любой плавающей операции примерно одно и то же. На эту тему я, к сожалению, могу лишь задавать вопросы.
← →
Jeer © (2007-06-09 15:16) [30]
> palva © (09.06.07 14:45) [29]
Деление и сегодня медленней, чем умножение.
В Вашем примере [27] замена на обратную величину и умножение приводит к результату
было 2500 ticks
стало 2100 ticks
← →
palva © (2007-06-09 15:20) [31]Jeer © (09.06.07 15:16) [30]
Ух ты! Лоханулся. Забыл суровую школу работы с матрицами.
← →
palva © (2007-06-09 15:22) [32]И видимо, деление заметно медленнее, раз такой эффект. Там ведь всего 8 делений!
← →
Jeer © (2007-06-09 15:25) [33]Да, всего лишь замена на:
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;
и такой результат:))
← →
Jeer © (2007-06-09 15:27) [34]Для меня было откровением, что это быстрее
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
← →
palva © (2007-06-09 16:40) [36]
Sapersky (09.06.07 16:03) [35]
> Вот, кстати, один товарищ пишет, что Дельфи вообще толком не умеет оптимизировать floating-point вычисления:
Я попробовал сишный вариант этой же программы (миллион обращений, поправка Jeer © (09.06.07 15:25) [33] учтена).
Borland С++ 5.5 дал 312 тиков, Delphi 7 - 281 тик.
Все компиляции производились с ключами по умолчанию.
← →
Jeer © (2007-06-09 16:48) [37]palva © (09.06.07 16:40) [36]
Интересно, почему такая разница
D7 на 1E6 обращений - 880 ticks
← →
palva © (2007-06-09 16:53) [38]
> 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);
}
← →
ferr © (2007-06-09 17:09) [39]procedure Inv2(var M1, M2); передаётся указатель
static void Inv2(double M1[4][4], double M2[4][4]) передаётся весь массив вроде как. + static вроде не очень хорошо.
← →
palva © (2007-06-09 17:17) [40]> 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
← →
palva © (2007-06-09 17:39) [41]Нашел разницу. Delphi использует передачу параметров в регистрах, а C кладет в стек. Надо бы поставить в C __fastcall.
← →
palva © (2007-06-09 17:49) [42]> Надо бы поставить в C __fastcall.
Не помогло. Параметры стали передаваться через регистры, но работать стало даже медленнее - 328 тиков. Не пойму, в чем дело.
← →
ferr © (2007-06-09 18:01) [43]Может стоит включить режим RELEASE, если он не включен?
← →
Alx2 © (2007-06-09 21:45) [44]Так и не понимаю, чего добиваетесь с этой банальщиной. Скорострельность работы? Оно вроде бы давно и успешно на подобных задачках делается компиляторами.
Интересно погрешность было б минимизировать. Но лениво. Совсем не вижу мелодрамы "мир и матрица 4x4". :)
← →
@!!ex_ (2007-06-09 21:50) [45]> [44] Alx2 © (09.06.07 21:45)
Ну так и какой из приведенных вариантов на ваш взгляд лучше?
← →
Alx2 © (2007-06-09 22:00) [46]>@!!ex_ (09.06.07 21:50) [45]
Насчет погрешности? Думать надо. Не готов сказать сейчас.
Насчет скорости - я через Maple прогоняю формулы с минимизацией на промежуточные вычисления и потом на Intel C++ 8.0 с SSE2 - оптимизацией получаю шустрые вещи. Если сходимость есть - забиваю. Иначе - ад. :)
← →
palva © (2007-06-09 22:21) [47]> Так и не понимаю, чего добиваетесь с этой банальщиной.
Думаем, как уважить просьбу автора - сделать быстрее.
← →
@!!ex_ (2007-06-09 22:24) [48]Кстати, я хочу заметить, что я - олень, в лучшем смысле этого слова...
В своей процедуре выполняю преобразования...
И оставляю матрицу в старом виде...
ТОт у меня нифига не получаются правильные теневые объемы, если объект повернут... :))))))))
← →
Alx2 © (2007-06-10 11:16) [49]Тоже увлекся :)
Получилось инвертировать за 500 тиков на P4 HT
← →
Alx2 © (2007-06-10 11:30) [50]Ого! Inv4 идет за 376 тиков.
А код от Maple за 500 тиков.
Да, есть смысл думать :)
← →
Alx2 © (2007-06-10 11:36) [51]Прогнал код через встроенный в компилятор профилировщик, потом попросил его оптимизировать на этой базе.
Inv4 получился всего 220 тиков. Класс!
← →
palva © (2007-06-10 23:17) [52]Alx2 © (10.06.07 11:36) [51]
А какой это компилятор и какой был процессор?
← →
Alx2 © (2007-06-11 02:44) [53]>palva © (10.06.07 23:17)
процессор P4 HT 2800 Mhz
Компилятор Intel C++ 8.0 с опциями /G7 /QxW /O3 /Qprof_use с первоначальной подгонкой опцией /Qprof_gen
← →
Jeer © (2007-06-13 18:34) [54]
> Alx2 © (10.06.07 11:36) [51]
> Inv4 получился всего 220 тиков. Класс!
Действительно - "умненький" icl. Но это уже тяжелое орудие:))
На P4 2.8GHz без HT посредством IC v.5 получаем 450 тактов CPU.
D7 устойчиво дает 850 тактов.
Кстати, FreePascal v.2.0.0 хуже - 950 тактов
Страницы: 1 2 вся ветка
Форум: "Прочее";
Текущий архив: 2007.07.15;
Скачать: [xml.tar.bz2];
Память: 0.69 MB
Время: 0.05 c