Entwickler-Ecke

Algorithmen, Optimierung und Assembler - Performance bei Mathematischen Gleichungen


Nova - Mi 07.12.05 02:32
Titel: Performance bei Mathematischen Gleichungen
Hallo zusammen,

ich bin gerade dabei einen Wellensimulator (für Wasserwellen) zu schreiben.
Auf die Idee bin ich gekommen als wir im Physikunterricht das Thema Wellen, Wellengleichung usw. hatten.

Eigentlich funktioniert schon alles soweit wie ich mir das wünsche, jedoch
möchte ich die Performance noch weiter erhöhen wenn es geht. Ich weiß nicht wie
stark die Rechner an unserer Schule sind, aber ich gehe 100%ig davon aus das meiner
zuhause schneller ist, und da der schon für einen Rechenzyklus 1/4 Sekunde braucht...
ohje ich will gar nicht daran denken *g*

Also die ganze Geschichte sieht so aus:

user defined image

Oben Links ist die Graphische Darstellung der Welle aus der Draufsicht.
Die unterschiedlichen Farben stehen für unterschiedliche Höhen der Welle (siehe Legende)

Berechnet wird das ganze wie folgt.
Jeder Pixel wir durch 2 for schleifen aufgerufen und einzelnt bemalt.


Delphi-Quelltext 
1:
2:
3:
4:
5:
6:
7:
8:
9:
 
  ImW.Picture.Bitmap.PixelFormat:= pf24Bit;
  for i:=0 to ImW.Height-3 do
  begin
    p := ImW.Picture.Bitmap.ScanLine[(i+1)];
    for a:=0 to ImW.Width-3 do
    begin
      inc(p);
      p^ := ColorToBColor(ElonToColor(Surface3d[i,a].Elongation));    end;
  end;


Jetzt gibt es zwei Funktionen die IMO optimiert werden können.
Das ist einemal die ElonToColor funktion, welche die jeweilige Elongation eines Pixels (also die Auslenkung nach oben oder nach unten) in eine Farbe umwandelt


Delphi-Quelltext 
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:
 
function TWave.ElonToColor(Elongation: Extended) : TColor;
var
  Elon: Extended;
  Maxi: Integer;
  Red, Green, Temp: Extended;
  C: String;
begin
  if TBMinElon.Position <> TBMaxElon.Position then
  begin
    Maxi := Max(abs(TBMaxElon.Position),abs(TBMinElon.Position));
    Elon := abs(Elongation);
    Temp := Elon/Maxi * 255;

    if Elongation > 0 then
    begin
      Green := 255 - Temp;
      Red := 255;
    end
    else
    begin
      Green := 255;
      Red := 255 - Temp;
    end;

    C := '$0000';
    C := C + inttohex(round(Green), 2);
    C := C + inttohex(round(Red), 2);

    Result := StringToColor(C);
  end;
end;


Und dann gibt es noch die getElon funktion, welche die Elongation ausrechnet


Delphi-Quelltext 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
 
function TWave.getElon(X,Y,EmitterId: Integer) : Extended;
var
 E: TEmitter;
 A,d,l,t: Extended;
 tX,tY: Integer;
begin
  if EmitterId < NumEmitterList then
  begin
    E := EmitterList[EmitterId];
    A := E.A;
    tX := abs(X-E.X);
    tY := abs(Y-E.Y);
    d := sqrt(tX*tX + tY*tY);
    t := SimTime;
    l := E.l;

    if t*l*E.f>d then
    begin
      Result := A * sin( 2*PI*( (t*E.f) - (d/l) ) );
    end
    else Result := 0;
  end
  else Result := 0;
end;


Weiterhin gibt es die Routine wo die jeweilige Elongation eines jeden Pixels berechnet wird. Vielleicht kann man hier auch etwas an der performance drehen?


Delphi-Quelltext 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
 
  for i:=0 to 278 do
  begin
    for j:=0 to 278 do
    begin
      Surface3d[i,j].Elongation := 0;
      for k:=0 to NumEmitterList do
      begin
        e := getElon(j,i,k);
        Surface3d[i,j].Elongation := Surface3d[i,j].Elongation + e;
      end;
    end;
  end;


Jetzt laufe ich halt vor das Problem, dass ich mir solche Dinge zwar ausdenken kann, ich aber keine Ahnung von Code Optimierung habe, also habe ich gehofft ihr könntet mir helfen die mathematischen Rechnungen zu optimieren?

Danke schonmal im Vorraus

gruß, Nova


Moderiert von user profile iconraziel: Topic aus Delphi Language (Object-Pascal) / CLX verschoben am Mi 07.12.2005 um 09:20


digi_c - Mi 07.12.05 09:17
Vielleicht verschieb das mal nach Algorithmen/Optimierung

Respekt, das sieht ja schon verdammt gut aus :-)

Ich würde evtl. ein Array anstatt der Pixel benutzen, auch wenn das ordentlich Speicher frist(relativ ;-)) ist es doch um einiges schneller anstatt für jedes Pixel ganze Routinen zu durchlaufen.
Und das wird dann gerendert.

Aber ich bin (noch nicht) so der Optimierungsfreak :-)

Achso für sowas bietet sich natürlich MMX und SSE und somit Inline Assembler an :nut:


DaRkFiRe - Mi 07.12.05 09:34
Einzelne Funktionen abschalten und Laufzeit testen (Anzeige der Laufzeit haste ja schon gut hinbekommen).

Dann entsprechend optimieren.


Kannst natürlich auch probieren, diverse kurze Funktionen (wie abs) auch inline zu implementieren, dann sparst Du Dir Prozeduraufruf (Programm wird schneller).

Ansonsten rate ich Dir - wie digi_c - dazu, die Berechnung erst ALLE im Array durchzuführen und DANN erst zu zeichnen.


Ansonsten scheint das ja ein qualitativ sehr hochwertiges Programm zu sein.


Horst_H - Mi 07.12.05 12:37
Hallo,

Das sieht ja ganz super aus.
Anbei, was ich bisher habe aendern konnte.
Aber die Hauptzeit ist das Zeichnen.
Fuer 100 Durchlaeufe kompletter Berechnung (tStart= 10 Sekunden) mit 4 Emittern:
Ohne zeichnen 2.5 Sekunden(egal ob double oder single)
mit zeichnen 8.5 Sekunden.


Delphi-Quelltext 
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:
 
unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, ExtCtrls;

type
  TForm1 = class(TForm)
    Button1: TButton;
    procedure Button1Click(Sender: TObject);
  private
    { Private-Deklarationen }
  public
    { Public-Deklarationen }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}
type
  tKoor = 0..238;
  tEmitter = record
               X0,Y0 : Tkoor;
               Amplitude,
               Frequenz,
               Wellenlaenge,
               Phase : double;
               istStehend : boolean;//Wellenfront ausserhalb des Bildes
             end;
 tElongation = array [tKoor,tKoor] of double;
 //globale Variable
var
  Elongation : tElongation;
  EmitterList: array of tEmitter;
  MAXAmp : double;

procedure BerechneElongation(var Elon:tElongation;var Emitter:tEmitter;Zeitpunkt:double);
var
  x,y : tKoor;
  dx,dy,dy2,Abstand,
  Grenzabstand,
  Auslenkung,
  WellenFakt,
  sinFakt : double;
  Berechnungen : integer;
begin
  With Emitter do
    begin
    sinfakt := (2*pi*Zeitpunkt*Frequenz-Phase);
    Wellenfakt := -2*pi/WellenLaenge;
   //schon fertig ausgebreitete Welle
    If IstStehend then
      begin
      for y:= low(tKoor) to High(tKoor) do
        begin
        dy := Y-Y0;
        dy2 := sqr(dy);
        for x:= low(tKoor) to High(tKoor) do
          begin
          dx := x-X0;
          Abstand := sqrt(sqr(dX)+dy2);
          Auslenkung := Amplitude*sin(sinfakt+Abstand*Wellenfakt);                             Elon[x,y] := Elon[x,y]+Auslenkung;
          end;//for x
        end;//for y
      end// if istStehend
    else
      begin
      Berechnungen := 0;
      Grenzabstand := Zeitpunkt*Wellenlaenge*Frequenz;
      for y:= low(tKoor) to High(tKoor) do
        begin
        dy := Y-Y0;
        dy2 := sqr(dy);
        for x:= low(tKoor) to High(tKoor) do
          begin
          dx := x-X0;
          Abstand := sqrt(sqr(dX)+dy2);
          //Welle schon ausgebreitet?
          If Abstand <= GrenzAbstand then
            begin
            Auslenkung := Amplitude*sin(sinfakt+Abstand*Wellenfakt);
            Elon[x,y]:=Elon[x,y]+Auslenkung;
            inc(Berechnungen);
            end;//If Abstand
          end;//For x
        end;//for y
      //Sind alle Punkte berechnet worden
      if Berechnungen >= sqr(High(tKoor)-low(tKoor)+1) then
        Iststehend := True;
      end;
    end;//with
end;

procedure ElonToColor(Elon:tElongation);
var

  x,y : tKoor;
  Red, Green, Temp: byte;
  Fakt,skalefaktor : double;
  BitMap : Tbitmap;

begin
  BitMap := Tbitmap.Create;
  Bitmap.Height := High(tKoor)+1;
  Bitmap.width := High(tKoor)+1;

  skalefaktor := 255/MAXAmp;//Maxi := Max(abs(TBMaxElon.Position),abs(TBMinElon.Position));
  for y := low(tKoor) to High(tKoor) do
    for x:= low(tKoor) to High(tKoor) do
      begin
      Fakt := Elon[x,y]*skaleFaktor;
      Temp := ABS(round(Fakt));
      If fakt < 0 then
        begin
        Green := 255;
        Red := 255 - Temp;
        end
      else
        begin
        Green := 255-Temp;
        Red := 255;
        end;
      //Form1.Canvas.Pixels[x,y] := tcolor(Green*256+Red);
      Bitmap.Canvas.Pixels[x,y] := tcolor(Green*256+Red);
      end;
   BitBlt(Form1.Canvas.Handle,0,0,High(tKoor)+1,High(tKoor)+1,Bitmap.Canvas.Handle,0,0,SRCCOPY);
  BitMap.free;
end;

procedure TForm1.Button1Click(Sender: TObject);
const
  //denn lambda*f = const = Wellenausbreitungsgeschwindigkeit
  //Licht 3e8 ,luft ~330 , wasser ~1500 [m/s]
  C = 300;//Pixel pro Sekunde
var
  k,i : integer;
  Zeit: double;
  t0 : TdateTime;
begin
  // Vier Emitter vorgeben
  setlength(Emitterlist,4);
  With EmitterList[0] do
    begin
    X0 := Low(tKoor)+20;
    Y0 := Low(tKoor)+20;
    Amplitude := 5.0;
    Frequenz :=  10.0;//+0.1;
    Wellenlaenge := c/Frequenz;//Punkte,
    Phase := 0.0;
    istStehend := false;
    end;
  With EmitterList[1] do
    begin
    X0 := High(tKoor)-20;
    Y0 := High(tKoor)-20;
    Amplitude := 5.0;
    Frequenz := 10.0;//-0.1;
    Wellenlaenge := c/Frequenz;//Punkte,
    Phase := 0.0;
    istStehend := false;
    end;
  With EmitterList[2] do
    begin
    X0 := Low(tKoor)+20;
    Y0 := High(tKoor)-20;
    Amplitude    := 5.0;
    Frequenz     := 5.0;//+0.1;
    Wellenlaenge := c/Frequenz;//Punkte,
    Phase        := 0.0;
    istStehend   := false;
    end;
  With EmitterList[3] do
    begin
    X0 := High(tKoor)-20;
    Y0 := Low(tKoor)+20;
    Amplitude := 5.0;
    Frequenz := 5.0;//-0.1;
    Wellenlaenge := c/Frequenz;//Punkte,
    Phase := 0.0;
    istStehend := false;
    end;
  MaxAmp := 0.0;
  for k:=0 to length(EmitterList)-1 do
    MaxAmp := MaxAMp+EmitterList[k].Amplitude;

  t0 :=time;
  For i := 1 to 100 do
    begin
    Zeit := i/200+10.0;
    for k:=0 to length(EmitterList)-1 do
      begin
      BerechneElongation(Elongation,EmitterList[k],Zeit);
      Button1.Caption := Inttostr(i);
      end;
    ElonToColor(Elongation);
    fillchar(Elongation,sizeof(Elongation),#0);
    end;
  t0 := time-t0;
  Button1.Caption := FormatDatetime('hh:mm:ss.zzz',t0);
end;

end.


Die Berechnung der Auslenkung ist wohl nicht richtig, aber das sieht auf dem Bild nicht falsch aus.Die Rechenzeiten werden dadurch wohl kaum veraendert.

Gruss Horst
[EDIT]
Auslenkung ist jetzt richtig gerechnet :-) Hoffe ich doch, oder ?


digi_c - Mi 07.12.05 13:26
Hilft dir da evtl. Scanline weiter? Schau mal nach Implementierungen von Antialiasing z.B. im Easy Delphi Helper http://www.dsdt.info/eh/

Wie wäre es mit Multithreading? Ich kann schlecht abschätzen, ob das so parallelisierbar ist sprich, ob es auf die Reihenfolge ankommt.

Ich behaupte GDI+ hilft dir da auch nicht groß weiter, da es ja nur höhere Funktionen anbietet, gell?

Bau doch ansonsten nen linienweisen Aufbau rein, wie bei PovRay/3DMax/... dann sieht das hochwissenschaftlich aus :D


Horst_H - Mi 07.12.05 15:11
Hallo,

durch Einsetzen von

Delphi-Quelltext 
1:
 
Bitmap.PixelFormat :=pf32bit;

,passend zu meiner Monitoraufloesung, konnte bei mir die Geschwindigkeit der Ausgabe extrem beschleunigt werden.
Bei zwei Emittern und wieder 100 Durchlaeufen betraegt die reine Rechenzeit 1,6 s (auch besser geworden) und mit Ausgabe 4.8 s.
Also 0,048 s pro Bild aus 20 Bilder pro Sekunde.
Die Ausgabegeschwindigkeit des Bildes ist ja konstant, also bremmsen mehr Emitter das Verfahren nicht so dramatisch.
Mist ich habe Elongation komplett uebergeben und nicht als Referenz const ...
hier als Grauwerte.

Delphi-Quelltext 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
 
procedure ElonToColor(const Elon:tElongation);
var
  x,y : tKoor;
  Gray, Temp: integer;
  Fakt,skalefaktor : double;
begin
  skalefaktor := 127/MAXAMp;//Maxi := Max(abs(TBMaxElon.Position),abs(TBMinElon.Position));
  for y := low(tKoor) to High(tKoor) do
    for x:= low(tKoor) to High(tKoor) do
      begin
      Fakt := Elon[x,y]*skaleFaktor;
      Temp := ABS(round(Fakt));
      If fakt < 0 then
        Gray := 128-temp
      else
        Gray := 128+Temp;
//      Gray := 128+Round(Fakt);
      Bitmap.Canvas.Pixels[x,y] := tcolor(gray*(256*256+256+1));
      end;
   BitBlt(Form1.Canvas.Handle,0,0,High(tKoor)+1,High(tKoor)+1,Bitmap.Canvas.Handle,0,0,SRCCOPY);
end;

Also sind es noch 4,4 s fuer 100 Bilder.
Bitmap wird in Form.create erzeugt und in Form.destroy freigegeben.

Gruss Horst
[Edit]
Nochmals geaendert, aber mwerkwuerdige Laufzeiten, nehmen immer mehr zu????
Bei Erststart beginnt es mit 5.0 Sekunden und nach mehrfachem Durchlauf 6.7 Sekunden ohne das der Speicherbedarf steigt.


Horst_H - Do 08.12.05 16:30
Hallo,

scanline ist der Bringer.
Das Zeichnen per draw spielt nur noch eine kleine Rolle.
Bei mir dauert(1.8 Ghz Athlon64) dauert es mit Sinustabelle 25s mit D7 bzw. 33s bei D2005 fuer 1000 Bilder mit Sinus gerechnet um 30s bzw 46s.
Also statt 0.25 s pro Bild nur noch 0,025 s also ein Faktor 10.
Mit 40 Bildern die Sekunde laesst es sich gut anschauen.
Die Ueberlaegerung mit einem anderen Bild hat merkwuerdige interessante Effekte.
Anbei der Quelltext.

@Nova: Sag mal etwas dazu.

Gruss Horst
P.S.:
Nach ein bisschen Zeigerarithmetik beim rechnen.
Die Schleife mit den 1000 Durchlaeufen (Button3Klick) braucht jetzt pro Bildpunkt ohne Anzeige 110 CPU-Takte pro BildPunkt Delphi 7 (125 mit Anzeige und SRCINVERT) mit sin-Berechnung sind es 154.5 Takte
Die 1.8e9 sind durch die eigene CPU-Takrate zu ersetzen.

Delphi-Quelltext 
1:
 
Label4.Caption := Format('%10.7f Takte',[t0*86400/1000/sqr(dist)/length(EmitterList)*1.8e9]);


Delphi 2005 ruft Int und sqrt und anderes viel umstaendlicher auf, das kostet zig- Takte macht 169 Takte pro Bildpunkt.
Das schon ganz schoen enttaeuschend, aber vielleicht genauer.
Mit sin-berechnug sind es extreme 253 Takte pro Bildpunkt

P.S.S 9.12.2005
Jetzt funktioniert es mit 91 Bildern pro Sekunde bei 240x240 Bildpunkten und 4 Emittern mit Sin Tabelle
Jetzt 80 Takte pro Pixel fuer die Berechnung , round(int(Abstand)) durch round(Abstand-conHalb) mit
const conHalb : double= 0.499999999, ersetzt.
Ca. 50 Bilder pro Sekunde werden mit sin Berechnung erreicht. Diese ist nur bei kurzen Wellenlaengen (In Pixeln auf dem Bildschirm) noetig.Eine Wellenlaenge kleiner als 2 Pixel zeigt sonst ein topfebenes Gewaesser.
Man kann es noch etwas beschleunigen wenn man Assembler einsetzt(immer dieses Wegspeichern und wiederholen der gleichen Daten), SSE(4 Wurzeln parallel in 38 Takten) oder 3dnow( 1.Naehreung in 5 Takten(rel. 2.5e-5 genau) bieten sich wegen der Wurzelberechnug an.

P.S.S....
Fuer D2005 muss man sqrt moeglichst vermeiden und mittel der Heron-Formel geht es ja auch ganz gut.
Statt 50 fps jetzt 89,7 fps disselbe Vesion unter D7 nur 87,6 fps, also fuer beide etwa gleich gut.
Die absolute Abweichung ist meist <1e-2 bei eine Amplitude von 4 ist das zu vernachlaessigen.
Das ist aber stark von der Anzahl der Wellenlaenge abhaengig, da die Tabelle nur Werte aufeinanderfolgender Bildpunkte beinhaltet.
Eine Wellenlaenge von 10 Punkten bedeutet, das der sin mit 10 Werten linear angenaehert wird und das ist nicht mehr allzu genau. Maximale Abweichung absolut : 0.17905 , relativ ~4,5%.
Bei einer Wellenlaenge von 25 ist nur noch absolut 0.03 , oder relativ 0.75%


Entwickler-Ecke.de based on phpBB
Copyright 2002 - 2011 by Tino Teuber, Copyright 2011 - 2025 by Christian Stelzmann Alle Rechte vorbehalten.
Alle Beiträge stammen von dritten Personen und dürfen geltendes Recht nicht verletzen.
Entwickler-Ecke und die zugehörigen Webseiten distanzieren sich ausdrücklich von Fremdinhalten jeglicher Art!