// Programm zum Testen der Gleitkomma-Rechenleistung
// Benutzung:
// 1. mit Compiler nach Wahl übersetzen, möglichst mit allen Optimierungen
// 2. Programm starten und gewünschte Zahl der Threads eingeben
// 3. ein paar Sekunden laufen lassen, bis sich die Ausgabe
// der Rechenleistung stabilisiert hat
// 4. ermittelte Leistung im Forum mitteilen
// Es kann nötig sein, unportable Programmteile zu ändern, vor allem den
// Multithreading-Teil. Die Alternative ist, den Multithreading-Teil
// ganz rauszuwerfen und das Programm mehrmals zu starten (mehrere Prozesse).
// Dazu einfach die Präprozessor-Anweisung unten ändern (#if 0 statt #if 1).
// Neben der ermittelten Leistung bitte mindestens folgende Angaben machen:
// - Zahl der Threads/Prozesse
// - Prozessortyp; ggfs. mit Angabe von Takt, TDP-Einstellung, Kühlung usw.
// - verwendeter Compiler und dessen Parameter (vor allem Optimierung)
// Idealerweise testet man sowohl mit einem Thread als auch mit so vielen,
// wie Kerne vorhanden sind, als auch, wenn HT/SMT vorhanden ist,
// mit doppelt so vielen.
// Das Programm berechnet die Brechung von Lichtstrahlen an kugelförmigen
// Linsenoberflächen. Es gibt die Leistung je Thread aus,
// in millionen Strahlen je Sekunde. Ein Strahl wird der Reihe nach an 8
// Flächen gebrochen, was jeweils rund 45 Gleitkommaoperationen erfordert
// (siehe die Methode System::strahl, die einen Strahl durchrechnet).
// 1 mio Strahlen je Sekunde bedeutet also 360'000'000 Gleitkommaoperationen
// je Sekunde (und je Thread).
#include <iostream>
#include <time.h>
#include <math.h>
using std::cout;
using std::endl;
typedef unsigned long long U64;
class Rand { // fancy PRNG
U64 a, b, c;
static U64 ro(U64 a, U64 b) { return a << b | a >> (64-b); }
void _rand() {
a += c++;
b ^= ro(a * 0x3874a66b48a6ae55, 32);
a ^= ro(b * 0x406a27534e6ef1db, 31);
b ^= ro(a * 0x194ba598311e3193, 23);
a ^= ro(b * 0x65f61398d08ac887, 11);
}
public:
Rand(U64 s) : a(s), b(time(0)), c(0) {}
double rand(bool sign) { // PRN in [0,1) or (-1,1)
const double f = 0.5 / ((U64)1 << 63);
_rand();
return a * (sign && (b & 1) ? -f : f);
}
};
//----------- der Optik-Teil -------------------------------------------------
const int N = 8; // Zahl der Flächen
const double T = 0.05; // Hauptstrahlsteigung
const double O = 2.0; // Öffnungsradius
const double F = 50.0; // Brennweite
struct Strahl {
const double cx0, cy0;
double y0, z0, y, z;
Strahl() : cx0(1. / sqrt(T*T+1.)), cy0(T * cx0) {}
};
class System {
double _k[N]; // Krümmungen der Flächen
double _d[N]; // Abstände der Flächen (_d[N-1] ist Schnittweite)
double _v[N], _v2[N]; // Verhältnis der Brechzahlen an den Flächen
double _ep; // Position der Eintrittspupille
public:
System();
double ep() const { return _ep; }
void strahl(Strahl &) const;
};
System::System() {
const double p[] = { // Daten des Linsensystems
3.951754785, 5.725546581, 3.010066534, 8.537160394,
7.609756865, 1.946574948, 4.66783622, 32.9165898,
0.05455588075, 0.01527275136, 0.07943846653, 0.1098046131,
-0.07929071194, -0.0381540381, -0.00233718517,
1.520665213, 1.6171631, 1.663401941, 1.986201015, 25.71689858
};
int j = 0;
for (int i=0 ; i<N ; ++i) _d[i] = p[j++];
for (int i=0 ; i<N-1 ; ++i) _k[i] = p[j++];
for (int i=1 ; i<N ; i+=2) _v[i] = p[j++];
_ep = p[j];
for (int i=0 ; i<N ; i+=2) _v[i] = 1. / _v[i+1];
double a = 1, c = 0;
for (int i=0 ; i<N-1 ; ++i) {
c = c * _v[i] + a * _k[i] * (_v[i] - 1);
a += _d[i] * c;
}
_k[N-1] = (1/F + c * _v[N-1]) / (a - a * _v[N-1]); // Krümmung letzte Fl.
for (int i=0 ; i<N ; ++i) _v2[i] = _v[i] * _v[i];
}
void System::strahl(Strahl & str) const {
// berechnet den Weg eines Strahls durch das System
double x = 0;
double y = str.y0;
double z = str.z0;
double A = str.cx0;
double B = str.cy0;
double C = 0;
for (int i=0 ; i<N ; ++i) { // Schleife über die Linsenoberflächen
double e = A*x + B*y + C*z;
double h = _k[i] * (x*x + y*y + z*z - e*e) + 2 * (e*A - x);
double R = A * A - h * _k[i];
double E = sqrt(R);
double d = h / (A + E) - e; // Abstand zum vorherigen Schnittpunkt
x += d * A;
y += d * B;
z += d * C; // (x,y,z) == Schnittpunkt Strahl mit Fläche
double g = sqrt(1 + _v2[i] * (R - 1)) - E * _v[i];
double a = g * _k[i];
A = A * _v[i] - x * a + g;
B = B * _v[i] - y * a;
C = C * _v[i] - z * a; // (A,B,C) == normierter Strahlvektor
x -= _d[i];
}
x /= A;
// Schnittpunkt mit Bildebene:
str.y = y - x * B;
str.z = z - x * C;
}
//----------------------------------------------------------------------------
// Jeder Thread bekommt ein Objekt der Klasse Thread_t (mit einem eigenen
// Objekt der Klasse System und einem eigenen Zufallsgenerator)
// und ruft darauf die Methode test() auf.
class Thread_t {
const System sys;
Rand rn;
const time_t t0;
public:
const int threads; // Zahl der Threads
const int num; // dieser Thread (0 bis threads-1)
Thread_t(int thr, int n, time_t t) : rn(n), t0(t), threads(thr), num(n) {}
void test();
};
void Thread_t::test() {
const time_t interval = threads * (threads > 4 ? 1 : 2);
time_t tn = t0 + num*interval/threads;
const unsigned anz = 100;
Strahl s;
U64 c = 0;
for (unsigned n=0 ;;)
{
// Berechnung:
for (unsigned i=1 ; i <= anz ; ++i) {
// erzeuge zufälligen Strahl
s.z0 = O * rn.rand(false);
s.y0 = O * rn.rand(true) - sys.ep() * T;
// berechne Strahlweg:
sys.strahl(s);
double y = s.y - F * T;
// prüfe auf Fehler:
const double ylo = -983e-6;
const double yhi = -963e-6;
const double zlo = -5e-6;
const double zhi = 13e-6;
if (!(y > ylo && y < yhi && s.z > zlo && s.z < zhi)) {
// Negation nicht entfernen, damit auch NaN erkannt wird
cout << "error thread " << num << " after " << c+i;
cout << " rays: " << y << ", " << s.z << endl;
return;
}
}
c += anz; // gesamte Zahl der Strahlen
n += anz; // Zahl der Strahlen seit letzter Ausgabe
// Ausgabe:
const time_t t = time(0);
if (t >= tn) {
tn += interval; // Zeit für nächste Ausgabe
if (threads > 1) {
if (threads > 10 && num < 10) cout << ' ';
if (threads > 100 && num < 100) cout << ' ';
cout << "thread " << num << ": ";
}
cout.width(8);
cout << std::left << n/(interval*1e6) << " Mrays/sec" << endl;
//anz = ((unsigned)sqrt((double)n) + anz) / 2;
//if (anz < 1) anz = 1;
n = 0;
}
}
}
#if 1
#include <windows.h>
DWORD WINAPI thf(LPVOID p) {
Thread_t *t = (Thread_t *)p;
//cout << "thread " << t->num << " started" << endl;
t->test();
cout << " thread " << t->num << " ended" << endl;
return 0;
}
void main() {
// assert(sizeof(U64) == 8);
Thread_t *tv[256];
int th;
cout << "enter number of threads (1 ... 256):" << endl;
std::cin >> th;
if (th < 1 || th > 256) return;
time_t t0 = time(0);
for (int i=0 ; i<th ; ++i)
tv[i] = new Thread_t(th, i, t0);
for (int i=1 ; i<th ; ++i) {
DWORD tID;
if (CreateThread(0, 0, (LPTHREAD_START_ROUTINE)thf, (LPVOID)tv[i], 0, &tID) == NULL) {
cout << "CreateThread error: " << GetLastError() << endl;
abort();
}
}
thf((LPVOID)tv[0]);
}
#else
// Alternative ohne Multithreading:
void main() {
Thread_t th(1, 0, time(0));
th.test();
}
#endif
-- mit aggressiveren Compiler-Flags geht fast der gesamte Stack-Traffic weg, und weil 
