Zugriffsgeschwindigkeit auf mehrdimensionale Arrays

[xcs]

Hallo!

Ich führe numerische Berechnungen durch und bin auf folgendes Phänomen gestoßen, das ich mir auf Anhieb nicht gleich erklären kann:

ich definiere ein statisches Array:

Delphi-Quellcode:

var mein_Array : array[0..x,0..y,0..z,0..x,0..y,0..z] of double;

Ich kann in mein_Array also (x+1)(x+1)(y+1)(y+1)(z+1)(z+1) Gleitkomma-Werte abspeichern.
So weit, so gut.
Es wird im Programm folgende Schleife durchlaufen:

Delphi-Quellcode:

for i:=0 to 50 do for j:=0 to 2 do for k:=0 to 2 do
begin
for l:=0 to 50 do for m:=0 to 2 do for n:=0 to 2 do
begin
//mache irgendwas, bei dem auf mein_Array[i,j,k,l,m,n] zugegriffen wird
end;
end;

Was mich jetzt wundert, ist folgendes:
Wenn x=50, y=14, z=14, dann wird die Schleife um ein mehrfaches (ca. Faktor 5-10) langsamer durchlaufen, als wenn x=1274, y=2, z=2.
Die Anzahl der Zugriffe ist exakt die selbe und auch die Arraygröße ist exakt die gleiche.

Woran liegt das?

[jfheins]

Es könnte sein, dass diese 1 Gigabyte-Array ausgelagert wird (= auf Festplatte) und dadurch das lesen von weit auseinanderliegenden Daten langsamer ist.

Ist aber nur ne Vermutung - btw.: Wieviel Arbeitsspeicher gibts denn?

[hadschi92]

Ein 6 dimensionales Array? Ist das nicht etwas übertrieben? Das macht mich jetzt echt neugierig was du programmierst. Vielleicht könnten wir dann auch noch eine andere Lösung finden.

[xcs]

Hmmm, das mit dem Arbeitsspeicher wäre ne Möglichkeit.
Zwar wird nix auf die Platte ausgelagert..., aber was da nun genau bei der Adressierung gemacht wird und wie sich das Programm durch das Array "hangelt", weiß man halt nicht so genau (ich jedenfalls nicht ;-))

Also, glaubt mir, ich habe mir schon überlegt, ob ich das nicht umgehen könnte.
Ich benötige es zur Berechnung des elektrischen Potentials im 3dimensionalen Raum, ausgehend von einer Ladungsverteilung.
An jeder Stelle i,j,k ist eine Ladungsdichte innerhalb eines diskreten Volumens dx*dy*dz gegeben.

Wenn das elektrische Potential am Ort i,j,k berechnet werden soll, müssen alle anderen Ladungen an den Stellen l,m,n mit jeweils einem Faktor multipliziert und aufsummiert werden (der Faktor ist eine Funktion des Abstandes zwischen i,j,k und l,m,n sowie der dielektrischen Permittivität des Mediums). Weil der Abstand ja zeitlich konstant ist (und die Permittivität auch), berechne ich die Faktoren f(i,j,k,l,m,n) einmal zu Beginn und speichere sie in dem Array ab.

Eine etwas elegantere Lösung wäre vielleicht, die Poisson-Gleichung zu lösen, aber das ist aufwändig zu programmieren und ich habe ein paar Probleme beim Übergang zwischen verschiedenen Werten (wenn sich die Permittivität ändert).

Aber wenn hier zufällig jemand Erfahrung mit der Numerik elektrostatischer Problemstellungen hat, würde ich mich über Rat freuen
:-)

Schönen Abend!

[xcs]

Aber falls jemand von euch Lust und Zeit hat, probiert es doch selbst mal aus.
Deklariert das Array wie ich und lasst die Schleife durchlaufen wie ich, mit irgend einem Zugriff auf das Array,
z.B.

Array[i,j,k,l,m,n] := i*j*k*l*m*n

o.ä.

die Zeit des Durchlaufens lässt sich mit QueryPerformanceCounter fast auf die ms genau bestimmen:

Start, Ende, f : int64;
Zeit : double;

QueryPerformanceFrequency(f)

QueryPerformanceCounter(Start);

//Schleife ausführen

QueryPerformanceCounter(Ende)

Zeit := (Ende-Start)/f;


Denn vielleicht liegt ja nur an meinem alten Delphi 7 Compiler ...

[jfheins]

Ich weis jetzt nicht ob ich das rtichtig verstanden habe ...

Du hast einen dreidimensionalen Raum und benötigst das array[x1, y1, z1, x2, y2, z2] um "Abstandsdaten" zwischen jedem Punktepaar zu speichern?

Wie kompliziert ist eigentlich die Berechnung der Werte? Wenn sie einfach ist, könnte man vll. auf das Array verzcihten :stupid:

Du könntest ja mal versuchen, einen 1 Gigabyte großen Speicherblock zu reservieren - wenn das fehlschlägt ist das Array wahrscheinlich in mehrere Teile aufgeteilt und der Zugriff entsprechend langsamer (normalerweise geschieht der Zugriff in O(1) also in konstanter Zeit)

[Medium]

Uff... Ich sehe das richtig, dass du ein Vektorfeld in einem Raum mit gewisser Auflösung komplett, also über das gesamte Volumen am Stück speichern möchtest? Warum brauchst du zu jedem Zeitpunkt alle diese Daten, und schnell verfügbar? Du hast zwar geschrieben was du da berechnest, aber dem "Normalohr" wie meinem erschließt sich noch immer nicht zu was das ganze letztlich führen soll.
Je nach dem was das werden soll, ließe sich evtl. über Näherungsverfahren, oder HD-gestützte Datenhaltung (evtl. mit prefetching wahrscheinlich demnächst gebrauchter Segmente), oder gar auf den Verzicht des ganzen nachdenken.

Volumendaten komplett am Stück lassen oft (nicht immer) auf konzeptionelle Schwächen schließen, da seltenst wirklich alles immer in voller Genauigkeit und super schnell gebraucht wird - so gut wie nie eigentlich. Ganz davon abgesehen ist es halt auch meist ziemlich unpraktikabel wie du merkst :)

Fazit: Erklär doch bitte nochmal in möglichst einfachen Worten was du berechnest, und vor allem wie du diese Werte anschließend benötigst (Simulation von einem/mehreren Teilchen in dem Feld, Echtzeit, Punktproben wählen, partiel summieren - das wären z.B. so Dinge wo ich mir vorstellen kann dass man sie in einem Vektorfeld tun wollen könnte, mal so als Beispiel).

[himitsu]

Zitat:

Zitat von jfheins

Du könntest ja mal versuchen, einen 1 Gigabyte großen Speicherblock zu reservieren - wenn das fehlschlägt ist das Array wahrscheinlich in mehrere Teile aufgeteilt und der Zugriff entsprechend langsamer (normalerweise geschieht der Zugriff in O(1) also in konstanter Zeit)

also Delphi selber sollte das Array nicht aufsplitten (im virtuellen Speicher)

aber was Windows macht, ist damit natürlich nicht sicher (dieses kann den Block ja aufteilen und auch auslagern ... z.B. in 8 bzw. 64 KB-Schritten)

[Blup]

Die Faktoren in (l,m,n) sind ausgehend von jedem Punkt im Raum (i,j,k) immer gleich.
Verschoben werden eigentlich nur die Punkte, auf die diese Faktoren anzuwenden sind.

Man kann also mit zwei getrennten Arrays arbeiten.

[xcs]

Danke für die bisherigen Ideen.

Was ist mit "Reservierung eines Speicherblocks" gemeint? Mit der Deklaration des Arrays wird doch der Speicher schon reserviert...

Zur Problemstellung:

- es geht um Elektrodiffusion von Ionen zwischen zwei Elektrolyten, welche durch eine teildurchlässige Membran getrennt sind
- die Membran hat eine andere relative dielektrische Permittivität ("Dielektrizitätszahl epsilon_r") als Wasser
- um Regionen, in denen "wenig" passiert (also weit weg von der Membran), gröber aufzulösen, und somit Rechenzeit zu sparen, will ich die Möglichkeit eines nicht-äquidistanten Gitters offen lassen

Daraus folgt:
- die diskreten Volumina sind nicht immer gleich groß --> f(i,j,k,l,m,n) ungleich f(l,m,n,i,j,k)
- wenn zwischen zwei Punkten gleichen Abstandes die Membran liegt, ändert sich der "Abstandsfaktor" auch: f(i,j,k,i+a,j+b,k+c) ist nicht gleich für alle i,j,k

Somit lässt sich das Problem nur sehr schlecht auf ein Array geringerer Dimension reduzieren.
(im homogenen, äquidistanten Fall würde ein 3D-Array reichen, weil ja jeder Punkt mittels eines Vektors vom Ursprung aus beschrieben werden könnte)

Eine Berechnung des "Abstandsfaktors" (ohne Array) würde viel länger dauern als es abzurufen. Das würde dann so aussehen:

Delphi-Quellcode:

faktor := sqrt(sqr(pos_x[i,j,k]-pos_x[l,m,n]) + sqr(pos_y[i,j,k]-pos_y[l,m,n]) + sqr(pos_z[i,j,k]-pos_z[l,m,n]))/(4*PI*epsilon_0*epsilon_r)

wobei epsilon_r streng genommen wieder eine Funktion von (i,j,k) ist.

Achso, die Schleife muss für jeden Zeitschritt delta_t erneut durchlaufen werden, weil sich dann ja die Ionenkonzentrationen (ergo die Ladungsdichten) durch die Diffusion verändert haben. Deswegen ist ja die Geschwindigkeit so entscheidend.

[jfheins]

Du könntest einfach mal versuchen single statt double zu benutzen sofern dir dessen Genauigkeit reicht - dann ist das ganze nurnoch halb so groß ;)

Ansonsten kann ich nur nochmal fragen: Wieviel Arbeitsspeicher hat der Rechner denn?

Oder ne ganz verrückte idee: Du könntest das ganze Zeug auf die GPU auslagern - gerade wenn die Probleme stark parallelisierbar sind, sind moderne GPUs um eine vielfaches schneller als x86-CPUs ;)

[xcs]

Ich hab 4 GB RAM und Windows XP (also 32 Bit).
Würde es was bringen, eine 64Bit-Version zu benutzen?

[xcs]

Aber um noch einmal auf das ursprüngliche Problem zurückzukommen:

Weiß jemand, wie genau bei der Programmausführung der Array-Zugriff funktionert (Pointer / Adressierung etc.)

Das Phänomen tritt schon bei sehr kleinen Arraygrößen auf:

Delphi-Quellcode:

implementation

{$R *.dfm}

var
  array1 : array[0..50,  0..4,  0..4, 0..50,  0..4,  0..4] of double;
  array2 : array[0..50, 0..14, 0..14, 0..50, 0..14, 0..14] of double;

procedure TForm1.Button1Click(Sender: TObject);
var
  i,j,k,l,m,n : integer;
  start1, ende1, start2, ende2 : int64;
  performance_factor : double;
begin

  queryperformancecounter(start1);
  for i := 0 to 50 do for j := 0 to 2 do for k := 0 to 2 do
  for l := 0 to 50 do for m := 0 to 2 do for n := 0 to 2 do
  begin
    array1[i,j,k,l,m,n] := array1[i,j,k,l,m,n] + i*j*k*l*m*n;
  end;
  queryperformancecounter(ende1);

  queryperformancecounter(start2);
  for i := 0 to 50 do for j := 0 to 2 do for k := 0 to 2 do
  for l := 0 to 50 do for m := 0 to 2 do for n := 0 to 2 do
  begin
    array2[i,j,k,l,m,n] := array2[i,j,k,l,m,n] + i*j*k*l*m*n;
  end;
  queryperformancecounter(ende2);

  performance_factor := (ende2-start2)/(ende1-start1);

  showMessage(floattostr(performance_factor));
end;

Probiert es mal aus. Schon hier ist der Unterschied von einem Faktor 2-5 zu erkennen (das kann entscheiden, ob ich 1 Tag oder 5 Tage warten muss, bis mein Rechner fertig ist). Bei größeren Arrays bis zu Faktor 10...

[jfheins]

Normalerweise geht das ungefähr so:

array[12] wird zu (@array + 12)^

und array[1, 4] wird zu (@array + 1 * <länge 2.Dimension> + 4)^

Oder so ähnlich. Auf jeden Fall sollte der Zugriff auf ein beliebiges Element in O(1) erfolgen also in konstanter Zeit (= unabhängig von der Menge der Daten)

[Medium]

Zitat:

Zitat von xcs

Das Phänomen tritt schon bei sehr kleinen Arraygrößen auf:

So "sehr klein" sind die gar nicht.
array1: 13.005.000 Bytes (~12,4 MB)
array2: 1.053.405.000 Bytes (~0,89 GB)

In deiner 6(!)-fach Schleife kommst du letztlich auf immerhin 210.681 Durchläufe. Das ist zunächst mal nicht so übermäßig viel, aber wenn dann auch noch im Speicher geswapped werden muss (was bei ca. 1GB auch passieren kann wenn man 4GB RAM hat) kann das durchaus mal etwas dauern.

Über wie viele Sekunden für je eine solche Schleife reden wir hier eigentlich? Evtl. lässt sich daraus ja zumindest ableiten ob da grundsätzlich nicht etwas schief läuft.


Edit: Und ich habe noch immer nicht verstanden was du vor hast. Die meisten sind hier vermutlich maximal Hobbyphysiker (unterstelle ich einfach mal ganz frech), und können mit den Fachtermini bestenfalls erahnen was da abgeht. Zumindest geht's mir so :?. Auch ging nicht so wirklich daraus hervor, was du am Ende nun haben/machen willst. Das ist hier u.U. wichtig, da man je nach Zweck noch immer ganz anders optimieren kann.

[Blup]

Um auf die ursprüngliche Fragestellung zurück zu kommen, hier mein Test:

Delphi-Quellcode:

const
 x=50; y=14; z=14; // Variante 1
// x=1274; y=2; z=2; // Variante 2
var
 mein_Array : array[0..x,0..y,0..z,0..x,0..y,0..z] of double;

procedure TForm1.Button1Click(Sender: TObject);
var
  i, j, k, l, m, n, idx: integer;
  x1, x2: LongWord;
begin
  Button1.Caption := '';
  x1 := GetTickCount;
  for i := 0 to x do
    for j := 0 to y do
      for k := 0 to z do
        for l := 0 to x do
          for m := 0 to y do
            for n := 0 to z do
              mein_Array[i,j,k,l,m,n] := i * j * k * l * m * n;
  x2 := GetTickCount;
  Button1.Caption := IntToStr(x2 - x1);
end;

Vor dem Test einer Variante habe ich Windows jeweils neu gestartet.
Die Funktion wurde nach jedem Programmstart 8 mal wiederholt aufgerufen.

Code:

Durchlauf          1  |   2  |   3  |   4  |   5  |   6  |   7  |   8
-----------------------------------------------------------------------
Variante 1
1.Programmstart  15,0 | 82,3 | 17,8 |  2,1 |  1,5 |  1,7 |  1,5 |  1,3
2.Programmstart  12,5 | 53,5 |  4,6 |  1,7 |  1,3 |  1,3 |  1,3 |  1,2

Variante 2
1.Programmstart  15,7 | 85,3 | 24,0 |  2,5 |  1,8 |  1,7 |  1,7 |  1,7
2.Programmstart  11,1 | 37,6 |  3,9 |  1,9 |  1,7 |  1,6 |  1,7 |  1,7

Beide Varianten von x,y,z sind in etwa gleich schnell.
Entscheident ist hier das Auslagern des Speichers durch Windows.

In der Variante 2 müssen allerdings deutlich mehr innere Schleifen initialisiert werden.
Das könnte die geringen Zeitabweichungen erklären.

[himitsu]

jupp, vom Zugriff auf die Elemente sollte bei der For-Schleifchen-Kombination eh alle Elemente in der selben Reinfolge und auch noch sequentiell erfolgen ... bei beiden Varianten.

für dieses hier

Delphi-Quellcode:

for i := 0 to x do
  for j := 0 to y do
    for k := 0 to z do
      for l := 0 to x do
        for m := 0 to y do
          for n := 0 to z do
            mein_Array[i,j,k,l,m,n] := i * j * k * l * m * n;

und falls ich mich jetzt nicht vertan hab, 'ne Variante mit nur einer Schleife

Delphi-Quellcode:

type PDoubleArray = array[0..0] of Double;

for f := 0 to (x+1) * (y+1) * (z+1) * (x+1) * (y+1) * (z+1) do
begin

//g := f;
//i := g mod (z+1);  g := g div (z+1);
//j := g mod (y+1);  g := g div (y+1);
//k := g mod (x+1);  g := g div (x+1);
//l := g mod (z+1);  g := g div (z+1);
//m := g mod (y+1);  g := g div (y+1);
//n := g mod (x+1);  g := g div (x+1);

  i := f mod (z+1);  g := f div (z+1);
  j := g mod (y+1);  g := g div (y+1);
  k := g mod (x+1);  g := g div (x+1);
  l := g mod (z+1);  g := g div (z+1);
  m := g mod (y+1);  n := g div (y+1);

  PDoubleArray(@mein_Array)[f] := i * j * k * l * m * n;
end;

[Blup]

Wenn die Membran eine sehr dünne Ebene ist und das restliche Medium weitgehend homogen, dürfte epsilon_r ja nur davon abhängen, ob die Membran zwischen den beiden Punkten liegt?
Diesen Term könnte man dann relativ leicht bestimmen:
(4*PI*epsilon_0*epsilon_r)

Der andere Term beschreibt nur den Abstand zwischen dem Punkt [i,j,k] und dem Punkt [l,m,n]:
sqrt(sqr(pos_x[i,j,k]-pos_x[l,m,n]) + sqr(pos_y[i,j,k]-pos_y[l,m,n]) + sqr(pos_z[i,j,k]-pos_z[l,m,n])

Um diese Berechnung zu sparen, würde ich(wie schon gesagt) eine extra Array erstellen, das genau diesen relativen Abstand enthält.
Abstand[i - l, j - m, k - n]