effiziente numerische algorithmen und softwareentwicklung ...€¦ · effiziente numerische...
TRANSCRIPT
Effiziente numerische Algorithmenund Softwareentwicklung fürhochparallele Systeme
PD Dr.-Ing. habil. Harald Köstler
17.4.2015, Tag der Informatik
Contents
Überblick Forschungsbereich
Von der abstrakten Problembeschreibung zur Implementierung
Beispiel: Entrauschen von Bilddaten
2
Überblick
Ziele
Entwicklung von effizienten Lösungsverfahren für Problemstellungen aus der Datenanalyse
Entwicklung von Tools zur Erhöhung der Portabilität auf verschiedene Rechnerarchitekturen und der Produktivität bei der Softwareerstellung
Exploration
Validierung
Erkenntnis
Design
Problem
Model
4
5
Computational Science and Engineering
Anwendungsfelder
• Naturwissenschaften• Ingenieurwissenschaften• Medizintechnik• …
Angewandte Mathematik
• Numerik• Statistik• Optimierung
Informatik
• High Performance Computing
• Parallele Algorithmen
• Software Design
USE_SweepSection( getLBMsweepUID() )USE_Sweep()
swUseFunction(„LBM",sweep::LBMsweep,FSUIDSet::all(),hsCPU,BSUIDSet::all()); USE_After()//Communication
6
High Performance Computing (HPC)
Bildverarbei-tung
Simulation
Cluster
Gegebene Ressourcen
Schwache Skalierung
Physik
GPU, FPGA
Gegebene Datengrößen
Starke Skalierung
Medizin
Modelle und numerische Methoden
Probleme im HPC
Hardware: Moderne HPC Plattformen sind massiv parallelParallelität auf allen Ebenen: Intra-core, intra-node, inter-node
Heterogene Architekturen
7
Probleme im HPC
Numerische Algorithmen: Die Anzahl von unterschiedlichen Verfahren wächst, viele davon sind problemspezifisch
Die Auswahl des besten (?) Verfahrens ist nicht trivial
Parametrisierung der Verfahren ist häufig komplex
Beschreibungssprache ist die Mathematik
Aufgrund der Anwendungsdomäne ist Datenparallelität (meist) wichtiger als Nebenläufigkeit
8
Probleme im HPC
Software: Die Anforderungen der Anwendungen an die Software sind extrem
Performance: Entweder die Laufzeit sollte möglichst 0 sein (Echtzeitbildverarbeitung) oder die Problemgrößen möglichst unendlich (Simulation der Entstehung des Universums)
Portabilität: Neue Architekturen sollen möglichst schnell unterstützt werden
Produktivität: Neue, komplexere Modelle sollen einfach integrierbar sein
Die Softwareentwicklung findet meist in interdisziplinären Teams statt
9
10
Anwendungsorientierte Projekte
Mathematiker
Hardwarespezialist
Softwarespezialist
Anwender Mathematische Beschreibung
Lösungsverfahren
Parallele Implementierung (Framework)
Effiziente Implementierung auf spezieller Architektur
Anwendungsbeispiele
12
waLBerla Framework
13
Strömung in komplexer Geometrie
Godenschwager, Schornbaum, Bauer, Köstler, Rüde, A framework for hybrid parallel flow simulations with a trillion
cells in complex geometries. Proc. of the Int. Conf. on HPC, Networking, Storage and Analysis. ACM, 2013.
14
Ternäre eutektische Erstarrung
Hötzer, Bauer, Köstler, Nestler et al, Large scale phase-field simulations of directional ternary eutectic solidification, Acta
Materialia, 2015. accepted.
15
Starrkörpersimulation
Zweites newtonschesGesetz
Discrete Element Method
OpenCL + DirectX
Microsoft Kinect Sensor
amF
JM
Kuckuk, Köstler, Preclik, Interactive particle dynamics using OpenCL and Kinect, Int.
J. of Parallel, Emergent and Distributed Systems, 2013.
Variationelle Bildverarbeitung
Problemstellung: Finde Abbildung , so dass das Energiefunktional
minimiert wird.
16
Segmentierung von Muskelfasern
Ziel: Segmentierung der 3D Faserstruktur ausSchichtbildern eines Mausmuskels
17
Data provided by O. Röhrle, Universität Stuttgart from Dane Gerneke, Auckland Bioengineering Institute at the
University of Auckland, New Zealand
18
2D Bildsegmentierung
Röhrle, Köstler et al, Segmentation of skeletal muscle fibers for applications in computational skeletal muscle mechanics,
Comp. Biomechanics for Medicine, 2011.
19
3D Bildsegmentierung
Osterwald, Using computational techniques to analyze structural components within skeletal muscles, 2014
Industrieprojekte I
Babycare: Signalanalyse
Femcare: Segmentierung, Klassifikation
Source: http://www.pg.com/de_DE/20
Industrieprojekte II
GPU Computing: Bildregistrierung, High dynamic rangecompression, Bildrekonstruktion
Source: http://www.siemens.com/entry/de/de/
21
Industrieprojekte III
Bildverarbeitung Nanoskopie: Entfaltung, Segmentierung
Source: http://www.abberior.com/
22
Von der mathematischen Beschreibung zur
effizienten Implementierung
24
• Sebastian Kuckuk
• Harald Köstler
• Ulrich Rüde
• Christian Schmitt
• Frank Hannig
• Jürgen Teich
• Hannah Rittich
• Matthias Bolten
• Alexander Grebhahn
• Sven Apel
• Stefan Kronawitter
• Armin Größlinger
• Christian Lengauer
Projekt ExaStencils: gefördert von der Deutschen Forschungs-
gemeinschaft (DFG) im Rahmen des Schwerpunktprogramms 1648
(Software for Exascale Computing) -> http://www.exastencils.org
Softwareentwicklungsansätze in CSE
Ein Problem / Architektur – eine Implementierung Eine neue Software wird entwickelt, die sich ggfs. existierender
Bibliotheken bedient
Einfach zu spezialisieren und damit zu optimieren
Problemstellung muss häufig auftreten
Eine Bibliothek – mehrere Probleme / Architekturen Höherer Wartungsaufwand und Nutzer müssen geschult und
unterstützt werden
Schwierig die Anforderungen aller Nutzer zu erfüllen und gleichzeitig optimale Performance zu erreichen
Eine ganze Community profitiert davon
Ein Generator – mehrere Probleme / Architekturen Die Entwicklung des Generators (domänen-spezifische Sprache)
erfordert sehr hohen Aufwand
Problem-spezifische Optimierungen sind möglich
25
26
Domänen-spezifische Projekte
Mathematiker
Hardware Spezialist
Software Spezialist
Verschiedene Anwender Beschreibung in domänen-spezifischer Sprache
Automatische Auswahl der Algorithmen
Codegenerierung für konkrete Anwendung
Automatische Anpassung an konkrete
Architektur
DomänenwissenVarianten Modell
Domänenexperte
PDE
Operators::Laplacian(Data::solution) = Data::rhs
Ebenenmodell
27
28
Kontinuierliche Domäne
Ω
∂Ω
In welchem Raum-Zeit Bereich
soll die Berechnung durchgeführt
werden?
Physik: Temperaturverteilung?
Datanalyse: Wo ist der Ofen?
Ω
∂Ω
Spezialisierung: zeitunabhängig
29
Kontinuierliches Modell
Ω
∂Ω
onTu
infu
Suche eine Funktion u, die
jedem Punkt im Raum die
passende Temperatur
im stationären Zustand zuordnet
30
Diskrete Domäne
Ω
∂Ω
∂Ω
Spezialisierung: regulär
An welchen diskreten Punkten
soll u berechnet werden?
Ω
In welchen Teilbereichen
soll u berechnet werden?
31
Diskretes Modell
Ω
∂ΩAn welchen diskreten Punkten
soll u berechnet werden?
hh fAu
Die Matrix A beschreibt die direkten
Abhängigkeiten der Punkte
voneinander, diese sind nur lokal- In welcher Reihenfolge und wo sind
die Punkte abgespeichert?
- Welche Einträge hat A und an welchen
Stellen?
32
Parallelisierung
Ω
∂Ω
∂Ω
Spezialisierung: regulär
Zerlegung in Teilgebiete
Ω
33
Hierarchische Zerlegung
Ω
∂Ω
Einführung einer hierarchischen
Zerlegung um globalen
Informationsaustausch zu
beschleunigen
Lösen des Gleichungssystems
34
Lokales Lösen
Vergröberung Verfeinerung
Korrektur
Lösen auf der gröbsten Stufe
Mehrgitterzyklus
Skalierung von parallelen Mehrgitterlöser
35Köstler et al, A Geometric Multigrid Solver on Tsubame 2.0, in Efficient Algorithms for Global Optimization Methods in
Computer Vision, 2014.
Zielsetzung
1. Das Webinterface öffnen
2. Das kontinuierliche Modell in einer latex-ähnlichen Sprache beschreiben
3. Anwenden von symbolischer Algebra um das diskrete Modell abzuleiten
4. Beschreibung des numerischen Algorithmus in einer Matlab-ähnlichen Sprache
5. Automatische Übersetzung in C-ähnliche Sprache und low-level Optimierung für die konkrete Architektur
6. Berechnung der numerischen Ergebnisse in der HPC Cloud und Visualisierung im Browser
36
Funktionen auf Ebene 4 – Beispiel
37
function VCycle@coarsest ( ) : Unit /* coarse grid solver */
function VCycle@((coarsest + 1) to finest) ( ) : Unit
repeat up 1
Smoother@current ( )
UpResidual@current ( )
Restriction@current ( )
SetSolution@coarser ( 0 )
VCycle@coarser ( )
Correction@current ( )
repeat up 2
Smoother@current ( )
Lokaler Jacobi Löser
38
function Smoother@((coarsest + 1) to finest) ( ) : Unit
communicate Solution[0]@(current)
loop over inner on Solution@(current)
Solution[1]@current =
Solution[0]@current
+ ( ( ( 1.0 / diag ( Lapl@(current) ) ) * 0.8 )
* ( RHS@current - Lapl@(current) * Solution[0]@current ) )
…
Resulting Code (w/o basic Opt)
39
#include "MultiGrid/MultiGrid.h"
void Smoother_4()
exchsolutionData_4(0);
#pragma omp parallel for schedule(static) num_threads(8)
for (int fragmentIdx = 0; fragmentIdx < 8; ++fragmentIdx)
if (isValidForSubdomain[fragmentIdx][0])
for (int y = iterationOffsetBegin[fragmentIdx][0][1];
y < (iterationOffsetEnd[fragmentIdx][0][1]+17); y +=1)
for (int x = iterationOffsetBegin[fragmentIdx][0][0];
x < (iterationOffsetEnd[fragmentIdx][0][0]+17); x +=1)
slottedFieldData_Solution[1][fragmentIdx][4][(((y*19)+19)+(x+1))] =
(slottedFieldData_Solution[0][fragmentIdx][4][(((y*19)+19)+(x+1))]
+(((1.0e+00/fieldData_LaplCoeff[fragmentIdx][4][((y*17)+x)])*8.0e-01)
*(fieldData_RHS[fragmentIdx][4][((y*17)+x)]
-(((((fieldData_LaplCoeff[fragmentIdx][4][((y*17)+x)]
*slottedFieldData_Solution[0][fragmentIdx][4][(((y*19)+19)+(x+1))])
+(fieldData_LaplCoeff[fragmentIdx][4][(((y*17)+289)+x)]
*slottedFieldData_Solution[0][fragmentIdx][4][(((y*19)+19)+(x+2))]))
+(fieldData_LaplCoeff[fragmentIdx][4][(((y*17)+578)+x)]
*slottedFieldData_Solution[0][fragmentIdx][4][(((y*19)+19)+x)]))
+(fieldData_LaplCoeff[fragmentIdx][4][(((y*17)+867)+x)]
*slottedFieldData_Solution[0][fragmentIdx][4][(((y*19)+38)+(x+1))]))
+(fieldData_LaplCoeff[fragmentIdx][4][(((y*17)+1156)+x)]
*slottedFieldData_Solution[0][fragmentIdx][4][((y*19)+(x+1))])))));
…
Resulting Code (w/ basic Opt)
40
#include "MultiGrid/MultiGrid.h"
void Smoother_4()
exchsolutionData_4(0);
#pragma omp parallel for schedule(static) num_threads(8)
for (int fragmentIdx = 0; fragmentIdx < 8; ++fragmentIdx)
if (isValidForSubdomain[fragmentIdx][0])
for (int c0 = iterationOffsetBegin[fragmentIdx][0][1];
(c0<=(iterationOffsetEnd[fragmentIdx][0][1]+16)); c0 = (c0+1))
double* slottedFieldData_Solution_1_fragmentIdx_4_p1 =
&(slottedFieldData_Solution[1][fragmentIdx][4][(19*c0)]);
double* fieldData_RHS_fragmentIdx_4_p1 = &(fieldData_RHS[fragmentIdx][4][(17*c0)]);
double* slottedFieldData_Solution_0_fragmentIdx_4_p1 = …
double* fieldData_LaplCoeff_fragmentIdx_4_p1 = …
for (int c1 = iterationOffsetBegin[fragmentIdx][0][0];
(c1<=(iterationOffsetEnd[fragmentIdx][0][0]+16)); c1 = (c1+1))
slottedFieldData_Solution_1_fragmentIdx_4_p1[(c1+20)] =
(slottedFieldData_Solution_0_fragmentIdx_4_p1[(c1+20)]
+(((1.0e+00/fieldData_LaplCoeff_fragmentIdx_4_p1[c1])*8.0e01)*(fieldData_RHS_fragmentIdx_4_p
1[c1] -
(((((fieldData_LaplCoeff_fragmentIdx_4_p1[c1]*slottedFieldData_Solution_0_fragmentIdx_4_p1[(c1
+20)])+(fieldData_LaplCoeff_fragmentIdx_4_p1[(c1+289)]*slottedFieldData_Solution_0_fragmentId
x_4_p1[(c1+21)]))+(fieldData_LaplCoeff_fragmentIdx_4_p1[(c1+578)]*slottedFieldData_Solution_0
_fragmentIdx_4_p1[(c1+19)]))+(fieldData_LaplCoeff_fragmentIdx_4_p1[(c1+867)]*slottedFieldData
_Solution_0_fragmentIdx_4_p1[(c1+39)]))+(fieldData_LaplCoeff_fragmentIdx_4_p1[(c1+1156)]*slott
edFieldData_Solution_0_fragmentIdx_4_p1[(c1+1)])))));
…
41
Stochastische variable Wärmeverteilung
Erwartungswert der Lösung (100
samples, Monte Carlo simulation)GRMF Einzellösung
Loh, Solving Stochastic PDEs with Approximate Gaussian Markov Random Fields using Different Programming
Environments, 2014.
Zusammenfassung
Für eine klar abgesteckte Problemklasse ist es möglichaus der mathematischen Beschreibung automatisch eineeffiziente Implementierung für eine bestimmte Architekturabzuleiten
Domänenexperten werden benötigt, um verschiedeneDisziplinen zu verbinden und automatisierte“Übersetzungstools” bereitzustellen
Variationelle Bildverarbeitung ist beeinflusst von physikalischen Modellen als a priori Wissen
42
Acknowledgements
Funded by
Bundesministerium für Bildung und Forschung
KONWIHR. Bavarian project
DFG SPP 1648/1 – Software for Exascale computing
Industry
Supercomputing centers
http://www.exastencils.org/
43
References (full list on http://www.exastencils.org/)
Stefan Kronawitter et al. Domain-Specific Optimization of two Jacobi Kernels
and their Evaluation in the ECM Performance Model. Parallel Processing
Letters, 24(3):Article 1441004, 2014. Special issue: HiStencils 2014.
Harald Köstler et al. A Scala Prototype to Generate Multigrid Solver
Implementations for Different Problems and Target Multi-Core
Platforms. Computing Research Repository (CoRR), 2014. arXiv:1406.5369.
Christian Schmitt et al. ExaSlang: A Domain-Specific Language for Highly
Scalable Multigrid Solvers. In Proceedings of the 4th International Workshop
on Domain-Specific Languages and High-Level Frameworks for High
Performance Computing (WOLFHPC), 2014.
Alexander Grebhahn et al. Experiments on Optimizing the Performance of
Stencil Codes with SPL Conqueror. Parallel Processing Letters, 24(3):Article
1441001, 2014.
Sebastian Kuckuk et al. A Generic Prototype to Benchmark Algorithms and
Data Structures for Hierarchical Hybrid Grids. Proceedings of the
International Conference on Parallel Computing (ParCo), 2013.44
Thank you for yourAttention!
Questions?