Optimierung: Determinanten-Berechnung

BenBE · Verfasst: Sa 09.10.04 01:06

Wollte mal fragen, ob ihr irgendeine Möglichkeit seht, folgenden Source NOCH schneller zu bekommen.

Zur Erklärung stehen an entsprechenden Stellen die Pascal-Befehle dran.

markieren

Bei Omorphia unter maths\OMathMatrix.pas

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
203:
204:
205:
206:
207:
208:
209:
210:
211:
212:
213:
214:
215:
216:
217:
218:
219:
220:

Type
    //Indizes are handled Row, Col
    TMatrix = Array Of Array Of Extended;

    //Dimensions of a matrix
    TMatrixSize = Packed Record
        X: Integer;
        Y: Integer;
    End;

Function GetMatrixDim(Const M: TMatrix): TMatrixSize; Register;
//Begin
//  Result.X := 0;
//    Result.Y := Length(M);
//    If Result.Y <> 0 Then
//        Result.X := Length(M[0]);
//End;
//
//This routine is intended to avoid the internal size checks due to redundant code
//
//EAX  Ptr to the Matrix
Asm
  MOV    ECX, EAX

  //Test if the array is assigned
  TEST  EAX, EAX
  JZ    @@Undef

  //Read the number of rows
  MOV    ECX, DWORD PTR [EAX-$00000004]

@@Undef:
  //Store the number of rows in the result
  MOV    Result.TMatrixSize.Y, ECX

  //If no rows in the array, skip testing the number of cols
  TEST  ECX, ECX
  JZ    @@Undef2

  //Dereference the first row
  MOV   EAX, DWORD PTR [EAX]

  //Read the number of columns
  MOV    ECX, DWORD PTR [EAX-$00000004]
  
@@Undef2:
  //Store the number of columns
  MOV    Result.TMatrixSize.X, ECX

  //Copy the result into the right register
  MOV   EAX, Result
End;

Function MatrixDeterminant(Const M: TMatrix): Extended;
Var
    S: TMatrixSize;

    Procedure LocalRaiseExcept;
    Begin
        Raise EMathError.Create(FmtInvalidDim);
    End;

Asm
//  =>  EAX     Pointer to Matrix Structure
//  <=  ST(0)   Determinant of this Matrix

    PUSH    EAX
//  S := GetMatrixDim(M);
//  MOV     EAX, DWORD PTR [M]                  //EAX is already the right value
  LEA     EDX, DWORD PTR [S]
  CALL    GetMatrixDim

    //If (S.X = 0) Or (S.X <> S.Y) Then
    MOV     EAX, [S.TMatrixSize.X]
    TEST    EAX, EAX
    JLE     @@MatrixDeterminant_DimError
    CMP     EAX, [S.TMatrixSize.Y]
    JZ      @@MatrixDeterminant_CorrectFormat

@@MatrixDeterminant_DimError:
    //Raises an exception due to invalid input matrix dimensions
    //Raise EMathError.Create(FmtInvalidDim);
    CALL    LocalRaiseExcept            //Call the local function for raising the exception

@@MatrixDeterminant_CorrectFormat:
    POP     EAX

  PUSH    ESI
  PUSH    EDI

    CMP     DWORD PTR [S.TMatrixSize.X], 2
    JNZ     @@AnySizeMatrix

    //Result := M[0, 0] * M[1, 1] - M[0, 1] * M[1, 0];
  MOV     ESI, EAX                            //DWORD PTR [M]
    MOV     EDI, DWORD PTR [ESI+$00000004]
    MOV     ESI, DWORD PTR [ESI+$00000000]
    FLD     TBYTE PTR [ESI]
    FLD     TBYTE PTR [EDI+$0A]
    FMULP
    FLD     TBYTE PTR [ESI+$0A]
    FLD     TBYTE PTR [EDI]
    FMULP
    FSUBP

    JMP     @@DeterminantFinish

@@AnySizeMatrix:
//  Result := 0;
  FLDZ

  MOV     ESI, EAX                            //DWORD PTR [M]
  
  MOV     ECX, DWORD PTR [S.TMatrixSize.X]
  TEST    ECX, ECX
  JZ      @@DeterminantFinish
//  For X := 0 To S.X-1 do
//  Begin
@@LoopDiagShift:
  DEC     ECX

//    Z := X;
  MOV     EDX, ECX
//    Tmp := 1;
  FLD1
  
//    For Y := S.X-1 DownTo 0 do
//    Begin
  MOV     EAX, DWORD PTR [S.TMatrixSize.X]
@@LoopDiagSum:
  DEC     EAX
  
//      Tmp := Tmp * M[Y,Z];
  PUSH    EAX
  MOV     EDI, DWORD PTR [ESI+4*EAX]
  LEA     EAX, DWORD PTR [EDX+4*EDX]
  FLD     TBYTE PTR [EDI+2*EAX]
  FMULP
  
//      If Tmp = 0 Then
//        Break;
  FTST
  FSTSW   AX
  SAHF
  POP     EAX
  JZ      @@LoopDiagSum_Break
  
//      If Z <= 0 Then
  TEST    EDX, EDX
  JNZ     @@LoopDiagSum_End
//        Z := Z + S.X;
//  ADD     EDX, DWORD PTR [S.TMatrixSize.X]  //Liefert sowieso S.X, da EDX = 0
  MOV     EDX, DWORD PTR [S.TMatrixSize.X]

@@LoopDiagSum_End:
//      Z := Z - 1;
  DEC     EDX

//    end;
  TEST    EAX, EAX
  JNZ     @@LoopDiagSum
  
@@LoopDiagSum_Break:
  
//  Result := Result + Tmp;
  FADDP

//    Z := S.X - X - 1;
  MOV     EDX, DWORD PTR [S.TMatrixSize.X]
  SUB     EDX, ECX
  DEC     EDX
//    Tmp := 1;
  FLD1
  
//    For Y := S.X-1 DownTo 0 do
//    Begin
  MOV     EAX, DWORD PTR [S.TMatrixSize.X]
@@LoopDiagDiff:
  DEC     EAX
  
//      Tmp := Tmp * M[Y,Z];
  PUSH    EAX
  MOV     EDI, DWORD PTR [ESI+4*EAX]
  LEA     EAX, DWORD PTR [EDX+4*EDX]
  FLD     TBYTE [EDI+2*EAX]
  FMULP
  
//      If Tmp = 0 Then
//        Break;
  FTST
  FSTSW   AX
  SAHF
  POP     EAX
  JZ      @@LoopDiagDiff_Break
  
//      Z := Z + 1;
  INC     EDX
//      If Z >= S.X Then
  CMP     EDX, DWORD PTR [S.TMatrixSize.X]
  JB      @@LoopDiagDiff_End
//        Z := Z - S.X;
//  SUB     EDX, DWORD PTR [S.TMatrixSize.X]  //Liefert sowieso 0, da EDX = S.X
  XOR     EDX, EDX
//    end;
@@LoopDiagDiff_End:
  TEST    EAX, EAX
  JNZ     @@LoopDiagDiff
  
@@LoopDiagDiff_Break:

//    Result := Result - Tmp;
  FSUBP
//  end;
  TEST    ECX, ECX
  JNZ     @@LoopDiagShift

@@DeterminantFinish:
  POP     EDI
  POP     ESI
End;

Zeitwertung:
Durchschnittlich 0.0009 Sec (~2.8 Mio. Taktzyklen) auf einer AMD CPU bei 1.4 GHz. und einer zufällig generierten 1000*1000-Matrix.

Viel Spaß beim Unterbieten.

P.S.: Das Code-Fragment gehört zum Omorphia-Projekt und unterliegt daher den Bedingungen der LGPL. Zu finden ist dieses unter sf.net/projects/omorphia/

_________________

Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.

Alni · Verfasst: Mo 11.10.04 11:34

So, ich hab mir das jetzt mal in Ruhe angeschaut und keine möglichkeit gesehen das noch weiter zu optimieren. Ich denke wenn man es weiterbeschleunigen will, bleibt nur noch eine prozessorspezifische Optimierung übrig. Oder du kannst natürlich die Determinatenberechnung ohne Probleme auf mehrere Threads oder Prozesse aufteilen. Das würde aber nur bei P4s mit HT oder Multiprozessorsystemen etwas bringen. Bei deinem AMD, vermute ich, würden mehrere Threads die Berechnug eher verlangsamen.

_________________
MfG Alex

BenBE · Verfasst: Di 12.10.04 21:06

THX erstmal für die Info, aber die Optimierung für spezifische CPUs ist im konkreten Anwendungsfall glaube eher zweitrangig, da das allgemein gehalten werden sollte. Sonst sieht der Code noch grausamer aus, als sowieso schon der Fall ;-)

Hab ich mich selbst schon gewundert hab, ob da nicht wirklich noch irgendwo eine Zeile zu viel drin sein könnte :twisted:

Naja, solange es schneller als Delphi ist, ist's ja in Ordnung (das kriegt aber glaub ich eigentlich jeder hin, der ASM relativ gut kann) ;-)

Werd mich dann mal noch an andere Algos machen... Dazu meld ich mich aber noch in anderen Topics hier.

MfG,
BenBE.

_________________

Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.

BenBE · Verfasst: Fr 22.10.04 19:21

Ich hab den Quelltext der Routine mal aktualisiert, da sie einen Bug für S.X = S.Y = 2 enthielt: Ergebnis war immer 0

Für einzellige Matrizzen bleibt ein weitere Fehler, den ich demnächst auch noch entferne.

_________________

Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.

Alni · Verfasst: Sa 23.10.04 01:20

Kannst du bitte die Stellen hervorheben die du geändert hast? Ich frage nur deshalb da Spezialfälle oft die allgemeine Lösung verlangsamen, obwohl diese eigentlich nie wirklich gebraucht werden.

_________________
MfG Alex

uall@ogc · Verfasst: Sa 23.10.04 01:55

ohne genauer den source anzugucken, wenn du aber das vorher in delphi geschriben hast und dann gedebuugt hast und du in delphi

// For Y := S.X-1 DownTo 0 do

benutzt hast, dann gibbet auf jedenfall ne optimierung, da downto 3 mal so langsam ist wie to

BenBE · Verfasst: Sa 23.10.04 02:31

Also, zu den Änderungen

58-61(+),
64(+),
73-83(#),
86(+),
91-108(+)

+ = hinzugefügt
# = geändert

Zu 58-61: Geht nicht anders. Begründung findet sich in einem anderen Thread dieser Sparte.

@Allgemeine Lösung: Nutzt eigentlich die Rechenvorschrift für Determinanten, indem das Summe vom Produkt der lo-ru-Diagonalen minus der Summe vom Produkt ro-lu-Diagonalen gebildet wird.

@Rechengenauigkeit: Könnte jemand mal prüfen, wo die größten (Rechen-)Fehlerquellen liegen könnten? Hab diverse Unstimmigkeiten mit dem eigentlichen Ergebnis entdeckt und kann mir nur eine Ursache wirklich erklären, jedoch denk ich, kann diese für S.X < 8 vernachlässig werden.

@Oben erwähnter Bug: Süezialbehandlung S.X = 1 fehlt noch, folgt aber demnächst noch.

_________________

Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.

	Mitgliederliste
	Gruppen
	Das Team
	Richtlinien
	Synonyme