Pokro čilé programovací techniky v Matlabu
T. Kozubek
MI21 14.11. 2010
Katedra aplikované matematikyFakulta elektrotechniky a informatikyVŠB-Technická univerzita Ostrava
Obsah1. Matlab
1. Funkce, řízení toku2. Homogenní pole (číselné matice, řetězce)3. Struktury4. Pole buněk5. OOP
MI21 24.11. 2010
6. GUI7. Import, export, MEX
2. Fourierova řada (spektrální analýza)
3. Diskrétní Fourierova transformace (spektrální analýza)
4. Konvoluce (vyhlazování signálu)
5. Laplaceova transformace (řešení DR)
6. Úlohy mechaniky – MatSol, paralelní programování
Matlab – funkcefunction [mean_x,std_x,var_x]=statistics(x)% STATISTICS Spocte zakladni statistiky vstupniho v ektoru x% UZITI: [mean,std,var]=statistics(x)
% Test na pocet vstupnich argumentuif nargin~=1; error( 'Spatny pocet vstupnich argumentu!' ); end;
% Test na korektnost vstupnich argumentu (x je cise lny vektor)if ~isvector(x); error( 'Vstupni argument x neni vektor!' ); end;
MI21 34.11. 2010
if ~isvector(x); error( 'Vstupni argument x neni vektor!' ); end;if ~isnumeric(x); error( 'Vstupni argument x neni ciselny vektor!' ); end;
% Spocti zakladni statistikymean_x=mean(x); std_x=std(x); var_x=var(x);
end
% Podfunkce MEAN (upraveny vypocet stredni hodnoty)function f=mean(x)% Vypocet stredni hodnotyf=sum(x)/(length(x)+1);end
Matlab – homogenní poleVytvo ření matice
A=[1 2 32 3 45 6 7];
B=[0,0,0; 2,3,4; 5,6,7];
Přístup k prvk ům, rušení
A(1:2,3)=[1;3]
A(A<7)=0
A(2,:)=[]
MI21 44.11. 2010
B=[0,0,0; 2,3,4; 5,6,7];
C=[A,B; B,A]; S=[‘abc’;’def’]
A(2,:)=[]
Řídké maticee=ones(5,1); A = spdiags([-e 2*e -e], -1:1, 5, 5)B=sprand(5,4,0.1)I=[1 1 2 3]; J=[1 3 2 4]; V=[1 1.5 2 3.7]; S=sparse(I,J,V,5,6)
spy(S,'o',8); imagesc(S); colorbar;
Matlab - struktury• Struktura je heterogenní pole, jehož prvky jsou libovolné objekty –
instance tříd (řetězce, číselné pole, buněčné pole, struktury, atd). • Přístup k položkám je prostřednictvím názvů položek užitých v
definici struktury.>>student=struct('jmeno','Jan Madlo', 'predmet',‘LAM','znamky',[1 2 1]);
MI21 54.11. 2010
• Paměťové požadavky : není třeba souvislý blok paměti pro celou strukturu, ale pouze pro položky struktury.
• Využití : snížení počtu vstupních a výstupních argumentů funkce, zvýšení čitelnosti kódu.
• Přístup k položkám>>student.jmeno>>student.jmeno=‘Miro Vilik'
Matlab – pole buněk• Pole bun ěk je heterogenní pole, jehož prvky jsou libovolné objekty
– instance tříd (řetězce, číselné pole, buněčné pole, struktury, atd). • Přístup k položkám je prostřednictvím indexů buněk.>>c={‘Hello’,rand(3),student};
• Paměťové požadavky : není třeba souvislý blok paměti pro celé
MI21 64.11. 2010
• Paměťové požadavky : není třeba souvislý blok paměti pro celé pole, ale pouze pro jeho položky.
• Využití : uložení řetězců různé velikosti, uložení matic různého řádu u metody rozložení oblastí, nahrazení vstupních a výstupních argumentů funkce, zvýšení čitelnosti kódu.
• Přístup k položkám>>c{1}, c{1}(2:end), c{3}.jmeno>>c{1}=‘Ahoj’, c{3}.jmeno=‘Ivan Leden’
Matlab – OOP
MATLAB má zabudovánu kompletní podporu OO přístupu
• Jak systém roste, je čím dál složitější sledovat předávání dat mezi funkcemi.
MI21 74.11. 2010
• Opakování kódu copy&paste –redundance – nekonzistence.
• Deklarace třídních proměnných – předem se ví, jaké má proměnné.
• Srozumitelnost a přirozenost návrhu.
Matlab – OOPclassdef MyClass < MySuperClass
propertiesProperty1Property2 = sin(pi/12); % default value
endproperties (SetAccess = private, GetAccess = private)
StressStrain
end
MI21 84.11. 2010
endmethods
function value = get.Property1(obj)%optionally implement get value
endfunction obj = set.Property1(obj,value)
%optionally implement set valueend
endend
Dědičnost, přetěžování operátorů, konstruktory, proměnné s omezeným přístupem.
Matlab – GUI
1. Pomocí LAYOUT editoru (.fig, .m)
2. Pomocí Matlabovských příkazů
MI21 94.11. 2010
(figure, uicontrol, uibuttongroup, guidata, guihandles, uimenu, …)
Matlab – import, export, MEX
1. Podpora práce s textovými i binárními soubory2. Načtení obrázků, zvuků, videa, MS Excell, …
Načtení a zobrazení obrázkurgb = imread('hokej.jpg');
Načtení a p řehrání zvukového souboru[y,fs,nbits]=wavread('babycry.wav');
MI21 104.11. 2010
rgb = imread('hokej.jpg');image(rgb);
[y,fs,nbits]=wavread('babycry.wav');sound(y,fs); plot(y);
MEX – přilinkování k Matlabu samostatných funkcí i celých knihoven v C, C++, Fortranu, Javě
( ) CaCzzazf ∈∈= ,,
Lineární zobrazení
MI21 114.11. 2010
( ) ))arg()(arg()arg()arg( zaiziai ezaezeazazf +===
Fourierova řada (spektrální analýza){ }
( ) RtectfTLtf
TT
,TLTe
ntin
Znnti
∈=∈
=
∑∞+
∞−
∈
,),,0()(
perioda...rychlost,úhlová...2
)0(fcísystémortonorm..../
2
2
ω
ω
πω
( ) ,,1
ZndtetfcT
nti ∈∫= − ω
MI21 124.11. 2010
( ) ,,1
0
ZndtetfT
c ntin ∈∫= − ω
• Dvoustranné spektrum : {cn,ϕ
n}
Fourierova řada
( ) ,,1
0
ZndtetfT
cT
ntin ∈∫= − ω( ) ,, Rtectf nti
n ∈= ∑+∞
∞−
ω
( ) 1,...,0,11 1
0
21
0
2
−==∆∆
= ∑∑−
=
−−
=
∆∆
−Nnef
Ntetf
tNc
N
k
nkN
i
k
N
k
tnktN
i
kn
ππ
Diskrétní Fourierova transformace (spektr. analýza)
MI21 134.11. 2010
00∆ == NtN kk
1,...,0,)(1
0
21
0
2
−=== ∑∑−
=
−
=
∆∆ Nkecectf
N
n
nkN
i
n
N
n
tnktN
i
nk
ππ
1,...,0,1
,1,...,0,1
0
21
0
2
−==−== ∑∑−
=
−−
=
NnefN
cNkecfN
k
nkN
i
kn
N
n
nkN
i
nk
ππ
DFT
Konvoluce posloupností (vyhlazování signálu)
{ } { }∞=
∞= 00 , llkk ba
{ } { } { } { }lknMN
nn bacc *,1
0 =−+=
{ } { } { } { } ∑=
−∞
= =∗=n
kkknnnnlkn baccbac
00,,
{ } { } 1
0
1
0 , −=
−=
N
llM
kk ba
MI21 144.11. 2010
==
−
−
−
−−
−
1
2
2
1
0
1
21
01
01
0
2
1
0
000
00
00
000
0000
00000
0
0
0
N
N
M
MM
M
b
b
b
b
b
a
aa
aa
aa
a
a
a
a
Abc
⋮
⋮
⋯⋯
⋯⋯
⋮⋱⋮⋱⋮⋮
⋯⋯
⋯
⋯
⋯
⋮
( ) ( ) ( ) ,tfyayayaya nn
nn =+′+++ −
− 011
1 ⋯
( ) ( ) ( ) ( )10 0 0 1 0 1
nny t c , y t c , ,y t c .−−′= = =⋯
( ) ( )∞
Laplaceova transformace ( řešení DR)
MI21 154.11. 2010
L(DR)DR Y(s) y(t)
( ) ( )∫∞
−==0
,.}{ dtetfsFfLT st
ANSA
modelování geometrie
ImporterSTL, IGIS,
Specifikace problému a generování sít ě.
(typ úlohy, materiálové vlastnosti, počáteční a okrajové podmínky, …)
Exporter
MI21 164.11. 2010
CAD System Integration
•Autodesk Inventor
•CATIA
•Pro/ENGINEER
•Solid Edge
•SolidWorks
•Unigraphics
STL, IGIS, VRML, … CDB format
T. Kozubek, T. Brzobohatý, A. MarkopoulosZ. Dostál, V. Vondrák, R. Ku čera, M. Sadowská
knihovna se škálovatelnými algoritmy založena na metodách rozložení oblasti FETI a BETI
Fproblem
FETI
1.podoblasti jsou volné
nebo uchycené
(předem není znám
defekt matic K)
F F
Metody rozložení oblasti
MI21 174.11. 2010
FETI-
DP
TFETI
2.
3.
defekt matic K)
FETI-DP (částečné
rozdělení, regulární
matice tuhosti)
všechny podoblasti
defekt matic K je znám
F
F
F
F
Fproblem
FETI
1.subdomains are fixed or free but
with different defectsF F
BETIIntroduction to Total FETI/BETI
MI21 184.11. 2010
FETI-DP
TBETI
2.
3.
FETI-DP (partial splitting,
nonsingular)
všechny podoblasti
jsou volné se shodnou
dimenzí jádra
FF
F
fuKuu TT −21min I I
G G
B B
B u c
B u c
B u c
≤==
u pole posuvůK matice tuhostiB matice omezeníc vektor omezení
funkcionál energie
Total FETI – primární proměnné
MI21 194.11. 2010
FF
obstacle
c vektor omezení
FETITotalDirichletovy okrajové podmínky jsou také
vynuceny pomoci LM, B
Bu= c
B
Vše a mnohem více m ůžete najít Table of Contents - Preface.
Part I. Background1. Linear Algebra.- 2. Optimization.
Part II. Algorithms3. CG for Unconstrained Minimization4. Equality Constrained Minimization 5. Bound Constrained Minimization 6. Bound
MI21 204.11. 2010
Bound Constrained Minimization 6. Boundand Equality Constrained Minimization
Part III. Applications to VariationalInequalities7. Solution of a Coercive VariationalInequality by FETI-DP method8. Solution to a Semicoercive VariationalInequality by TFETI Method.- References.-Index.
Seminář výpočetní mechaniky, 1-3. prosince 2010, VŠB-TU & Ústav geoniky AV ČR, v.v.i.
Škálovatelnost: kostka nad tuhou překážkou
MI21 214.11. 2010
MI21 224.11. 2010
23
OOK ⋮
1
blokově diagonální struktura
matice tuhosti =>
vhodná pro paralelní
implementaci
1. ingredience: paralelní implementace
MI21 234.11. 2010
=
nKOO
OKO
OOK
K
⋮
⋯⋱⋯⋯
⋮
⋮
2
Mohou být řešeny koercivní
a semikoercivní úlohy.
• Paralelně smyčky for– mnoho cyklů
Paralelní programování
Typické případyCluster HP model BLc7000 (c-class).
Konfigurace výpo četních uzl ů:2x dual core CPU AMD Opteron 2210 HE 8GB ECC DDR2 667MHz RAM
MI21 244.11. 2010
– dlouhé iterační cykly
• Offloading Work• Objemné datové sady
IT4Innovations, Centrum excelence, Ostrava, Czech Republic, cluster with more than 30,000 processors.http://www.it4i.eu
Matlab Distributed Computing Server
ClientParallel
Scheduler
WorkerDistributed
Computing Server
Worker
MI21 254.11. 2010
ParallelComputing
Toolbox
Scheduleror
Job manager
WorkerDistributed
Computing Server
WorkerDistributed
Computing Server
• Local job manager
• MathWorks job manager
• Third-party schedulers
(Windows CCS, PBS Pro, LSF)
Job
Job
Job
Job
Job
Job
Scheduler
Pending
Queued RunningWorker
Worker
Life cycle of a job
MI21 264.11. 2010
Client
Job
Job
Job
Job
Job
Job
Job
Job
Job
JobgetAllOutputArguments
submitFinished
Worker
Worker
Worker
createJob
Distributed computations1. Find scheduler or
job manager2. Create a job3. Create tasks and
associate them withthe job
>>sched=findResource(‘scheduler‘,
’type’,’local’);
>>job=createJob(sched);
>>createTask(job,@rand,1,{3,3});
>>createTask (job,@eye,1,{4});
MI21 274.11. 2010
the job
4. Send job to the front
5. Wait until job finishes
6. Gather results7. Destroy job
>>createTask (job,@eye,1,{4});
>>createTask (job,@ones,1,{{4},{3}});
>>submit(job);
>>waitForState(job);
>>results=getAllOutputArguments(job )
>>destroy(job);
Parallel computations>>sched=findResource(‘scheduler‘,…’type’,’local’);
>>pjob=createParallelJob(sched);set(pjob,‘MaximumNumberOfWorkers,4);set(pjob,‘MinimumNumberOfWorkers,4);
1. Find scheduler or job manager
2. Create a parallel job and sets the numer of tasks
3. Create parallel tasks
MI21 284.11. 2010
>>createTask(pjob,@MyParTask,1,{});
>>submit(pjob);
>>waitForState(pjob);
>>results=getAllOutputArguments(pjob)
>>destroy(pjob);
MPT functions: numlabs, labindex, labBarrierlabBroadcast, labProbe, labReceive, labSend,
labSendReceive
tasks4. Send job to the
front5. Wait until job
finishes6. Gather results7. Destroy job
Parallel computation scenario
MasterSlave 1
(subdom. 1)Slave 2
(subdom. 2)Slave 3
(subdom. 3)
Send FEM
Assemble K,K+,R,f
Build
Send FEM models to workers and parallelize assembling of K, K+, R, fv, stresses, searching
contact pairs, multiplication procedures.
MI21 294.11. 2010
Build
contact
conditions Gather f
Multiply by K+
Gather resultIterate
Compute stress
Gather stress
X
Důlní pr ůmysl: sv ěrný spoj d ůlní výztuže
MI21 304.11. 2010
tlak vyvolaný okolním materiálem
ocelová výztužpřekryv
R1
R2
t2
t1
Coulombovské t ření
f = 0.1
E = 2.1e5 MPa,
µ = 0.3
107m
m
BB
C
F
AA
B
C
A
Důlní pr ůmysl: sv ěrný spoj d ůlní výztuže
MI21 314.11. 2010
rovina symetrie
tA = tD = 10 mm
tB = tC = 40 mm
107m
m
F
DD
počet:
binárních proměnných 65562
duálních proměnných 3112
D
Důlní pr ůmysl: sv ěrný spoj d ůlní výztuže
MI21 324.11. 2010
8 těles7 volných
redukované napětí HMHPosunutí
2D
Důlní pr ůmysl: sv ěrný spoj d ůlní výztuže
MI21 334.11. 2010
3D
1NODAL SOLUTION
bez korekcí kontaktních směru 3 korekce
σN = 985.4 MPa σN = 1020.6 MPa
MI21 344.11. 2010
MN
MX
PR03.db FKOP= -0.04 f=0.0 PLANE182
.012349126.416
252.821379.225
505.629632.033
758.437948.043
DEC 9 200811:06:26
PLOT NO. 8
STEP=2SUB =34TIME=.5CONTSTOT (AVG)DMX =.100117SMN =.012349SMX =948.043
ANSYS (paralel)dynamický výpočet: 27 hodin
ANSYS (paralel) statický výpočet: 1 hodina
MATSOL statická úloha, 3 korekce:
sekvenční verze 15 min.
paralelní 7 min.
TBETI TFETI1.548
1.611
0.774
0.387Celková posunutí [mm]
Srovnání TBETI a TFETI
MI21 354.11. 2010
TBETI TFETI
0.387
0.017
1.548
1.611
0.774
0.387
0.017
Celková posunutí [mm]
Čas výpočtu 2.41 h.
Celkový čas 2.55 h.
násobení matice
Čas výpočtu 2.54 h.
Celkový čas 3.33 h.
násobení matice
250 podoblastí f = 0.1
Srovnání TBETI a TFETI
MI21 364.11. 2010
násobení matice vektor
2 126
primární proměnné
1 592 853
duální proměnné 261 553
max(abs(u)) 1.303 mm
max(totalU) 1.548 mm
násobení matice vektor
1 882
primární proměnné
713 751
duální proměnné 261 553
max(abs(u)) 1.378 mm
max(totalU) 1.616 mm
FETI BETI
vnější kroužek
klec
vnitřní kroužek
hřídel
zatížení
Automobilový pr ůmysl: kuli čkové ložisko
MI21 374.11. 2010
uchycení
3
Válečkové ložisko v ětrné elektrárny
MI21 384.11. 2010
Statisticspočet těles 73
počet podoblastí 700
primární prom ěnné 2,73 M
duální prom ěnné 459,8 k
počet iterací 4270
Automobilový pr ůmysl: pneumatika 3.9 mil stupňů volnosti,4000 iterací
MI21 394.11. 2010
celková posunutí redukované napětí HMH
Světlomet – Visteon (Autopal)
12 těles
MI21 404.11. 2010 VŠB – TU Ostrava
Světlomet – Visteon
64 Oblastí
MI21 414.11. 2010 VŠB – TU Ostrava
Světlomet – Visteon (Test solution-Total Displacement)
MI21 424.11. 2010 VŠB – TU Ostrava
Sestava složená z 29 těles
Dveře auta
MI21 434.11. 2010 VŠB – TU Ostrava
Dveře auta
MI21 444.11. 2010 VŠB – TU Ostrava
Dveře auta
128 Oblastí
MI21 454.11. 2010 VŠB – TU Ostrava
Dveře auta (Test solution-Total Displacement)
MI21 464.11. 2010 VŠB – TU Ostrava
Auto – Sestava obsahuje okolo 700 těles
MI21 474.11. 2010 VŠB – TU Ostrava
3000 subdomains
MI21 484.11. 2010 VŠB – TU Ostrava
700 bodies
MI21 494.11. 2010 VŠB – TU Ostrava
3000 subdomains
Rozšíření na d ynamic ké úlohy
MI21 504.11. 2010
Děkuji za pozornost
MI21 514.11. 2010