Lineární algebra¶
Typický numerický program je nabitý algebraickými operacemi s poli. Tehdy jsou vítané syntaktické možnosti Fortranu k manipulaci s poli bez cyklů (array language) – tím začneme. Pokud však bude záležet na rychlosti běhu programu, nemusí být interní implementace těchto operací dost rychlé a je pak záhodno použít optimalizované procedury z numerických knihoven. Ze soutěže numerických knihoven se již před časem vyvinulo jednotné rozhraní pro operace lineární algebry. Několik knihoven respektujících toto rozhraní, pro CPU i GPU, pro Linux i Windows, s několika fortranskými překladači, si ukážeme v chodu na vlajkovém příkladu – maticovém násobení větších plných matic.
Maticové násobení¶
Operace lineární algebry na vektorech a maticích patří k frekventovaným stavebním kamenům numerických algoritmů. Jejich rychlost je klíčová, zároveň je žádoucí jejich jednoduchý zápis ve zdrojovém kódu. Operace zahrnující vektory jsou obvykle dostatečně rychlé i bez další snahy, a lze tak těžit z jednoduché syntaxe Fortranu pro práci s poli:
! Fortran array language (no indexed loops)
real,allocatable :: a(:),b(:),c(:)
n=1000
allocate (a(n),b(n)) ! allocation of two vectors
a=1. ! assignment of a scalar to all array elements
b=a ! elemental assignment (element by element)
c=a ! implicit allocation and elemental assignment
a=a+b; a=a-b; a=a*b; a=a/c; a=a**a ! elemental arithmetic operations
a=sqrt(a); a=log(exp(a)); a=atan(tan(a)) ! elemental functions
a(1:n)=0.; a(2:n:2)=1. ! array sections (min:max[:stride])
print *,sum(a),count(a>0) ! reduction functions: real(n/2) int(n/2)
print *,product(a(2::2)),product(a,mask=a>0.) ! 1. 1.
print *,minval(a),minloc(a),maxval(a),maxloc(a) ! 0. 1 1. 2
print *,dot_product(b,b),sum(b*b),sum(hypot(b,0.)**2) ! 3x dot product
if (allocated(b)) deallocate (b) ! conditional deallocation
deallocate (c,stat=ierr); if (ierr/=0) print *,'ierr =',ierr ! secure deallocation
a=[(i,i=1,n)] ! array constructor with an implied-do list
print *,shape(a),size(a),lbound(a),ubound(a) ! inquiry functions: n n 1 n
allocate (b,source=a); allocate (c,mold=a) ! same-shape allocations with/without copying values
print *,a(1),b(1),c(1) ! 1. 1. undefined
call random_number(a) ! random numbers from <0,1)
print '(i0,a,i0,a,i0)',size(a),'=',count(a>0.5),'+',count(a<=.5) ! check of random-numbers distribution
b=merge(a,0.,a>0.5) ! merging arrays or scalars according to a mask
print '(*(f0.7,a))',sum(a),'=',sum(b),'+',sum(merge(mask=a<=.5,tsource=a,fsource=0.)) ! distribution of values
end program
Je-li počet aritmetických operací srovnatelný s počtem paměťových míst, jak bývá typické při operacích s vektory, je celková rychlost daná spíše rychlostí přístupu do paměti a algoritmus je memory-bound. Operace mezi maticemi mohou potřebovat více aritmetických operací, než kolik je zúčastněno paměťových míst, algoritmus je compute-bound, a záleží-li na rychlosti, může být žádoucí použít sofistikovanější postupy. Ukázkou je algoritmus pro maticové násobení,
Ten vyžaduje kubické množství aritmetických operací na kvadratickém množství prvků, vztaženo k počtu řádků či sloupců matic. Přímočará implementace maticového násobení, a často ani interní fortranská funkce matmul, neřeší kubický počet přístupů ke kvadratickému počtu prvků polí ekonomicky a jejich doba běhu (wallclock time) je zbytečně velká.
! matrix multiplication via loops
real,allocatable :: a(:,:),b(:,:),c(:,:)
mmax=1000; nmax=mmax; kmax=mmax
allocate (a(mmax,kmax),b(kmax,nmax),c(mmax,nmax))
a=1.; b=1.
do n=1,nmax; do m=1,mmax
t=0.; do k=1,kmax; t=t+a(m,k)*b(k,n); end do; c(m,n)=t
end do; end do
print *,c(1,1),sum(c)
! matrix multiplication via the matmul function
c=matmul(a,b)
print *,c(1,1),sum(c)
end program
Podstatnou roli, jako ostatně často, hraje schopnost fortranských překladačů překládaný kód optimalizovat, v případě funkce matmul dokonce nahradit jednoduchou interní implementaci sofistikovanou verzí z rychlých numerických knihoven. Procedury z těchto knihoven je samozřejmě možné volat i přímo.
Dostáváme se tak k dnešnímu tématu, kterým jsou numerické knihovny pro lineární algebru, jejich použití s Fortranem a ukázky porovnávající jejich rychlost na maticovém násobení matic. Vše pak můžeme porovnat s černými skříňkami Matlabu, Octavu a pythonovského NumPy, které přirozeně čerpají z týchž knihoven.
Numerické knihovny¶
Aritmetické operace 1. mezi vektory, 2. mezi maticemi a vektory, 3. mezi maticemi řeší knihovna BLAS Level 1, 2 a 3 (Basic Linear Algebra Subprograms, BLAS na wiki). Jde vlastně o definované rozhraní funkcí (Application Programming Interface, API), které pak různí tvůrci naplňují svými implementacemi. Tzv. referenční implementace obsahují elementární implementace bez nároků na rychlost (tedy jen s nároky na numerickou správnost) a nejsou žádoucí jinak než pro testování. Nás tak budou zajímat rychlé implementace BLASu (např. OpenBLAS, ATLAS) a také v numerické knihovně Intel MKL (Math Kernel Library) – s těmito knihovnami si ověříme explicitní propojení pro maticové násobení. Dále sáhneme po implementaci BLASu pro grafické karty (Graphics Processing Units, GPU), které se díky své rychlosti a zpřístupnění schopnosti počítat obecné numerické úlohy staly vítanou alternativou k běžným procesorům (Central Processing Units, CPU).
Budeme potřebovat fortranské rozhraní BLASu (definované je i podobné rozhraní pro jazyk C) k podprogramu dgemm pro maticové násobení:
! interface of the BLAS function for matrix multiplication of general double-precision matrices
subroutine dgemm(transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc)
Písmeno d v názvu podprogramu znamená verzi pro double precision neboli real(8) matice a naznačuje, že k dispozici jsou i verze sgemm pro single precision/real(4), cgemm pro complex(4) a zgemm pro complex(8). To se týká celého BLASu. Písmena ge stojí pro general neboli plné matice (v kontrastu k pásovým nebo třeba symetrickým pozitivně definitním maticím) a mm je matrix multiplication (nebo matrix × matrix). Podprogram poskytuje operaci c=alpha*op(a)*op(b)+beta*c, kde alpha a beta jsou skaláry, op() nabízí volitelnou transpozici nebo hermitovské sdružení a hvězdička mezi maticemi má význam maticového násobení, nikoliv (ve Fortranu realizovaného operátorem *) prvkového násobení (tedy prvek po prvku, elemental). K zjištění významu argumentů vede přímočará cesta přes zadání jména podprogramu v Googlu.
BLAS je nezbytnou součástí nadřazené knihovny LAPACK (Linear Algebra PACKage, LAPACK na wiki), přinášející řešiče (solvers) soustav lineárních algebraických rovnic a třeba také řešiče vlastních čísel (eigensolvers).
Připojování knihoven¶
Budeme pracovat s fortranskými překladači GCC Fortran (GNU Compiler Collection) alias gfortran, Intel Fortran alias ifort a ifx (Classic a LLVM-based) a Nvidia Fortran alias nvfortran (dřívější Portland neboli PGI Fortran, pgfortran). Všechny jsou dostupné v Linuxu, GCC a Intel i ve Windows.
V Linuxu budeme zdrojový kód z přednášky, a podobně i jiné fortranské kódy používající BLAS nebo LAPACK, překládat následujícími příkazy:
# compiler calls in Linux
# set variables LIBRARY_PATH and LD_LIBRARY_PATH
gfortran -O testMM.f90 -lblas
ifort -heap-arrays testMM.f90 -qmkl
nvfortran testMM.f90 -lblas
BLAS a LAPACK je přímočaré získat z linuxových repozitářů, obvykle v implementaci OpenBLAS. (Alternativou je ATLAS.) Intel poskytuje BLAS a LAPACK jako součást své knihovny Intel MKL, kterou lze jednoduše připojovat neboli linkovat volbou -qmkl nebo složitěji s pomocí stránky MKL Link Line Advisor, to pak i s jinými překladači. Default gfortranu je překlad bez optimalizace, je proto žádoucí uvádět patřičnou volbu explicitně, např. -O2, -O3 nebo agresivnější -Ofast; ostatní překladače v defaultu přiměřeně optimalizají, přidat lze např. volbami ifort -Ofast a nvfortran -fast. Volba Intelu -heap-arrays snižuje zatížení zásobníku, které se často projevuje chybou za běhu (Segmentation fault, segfault a podobně); obdobná voba gfortranu zní -fmax-stack-var-size=0, zde nezbytná při -Ofast. Přetékání zásobníku obecně lze v Linuxu čelit příkazem shellu ulimit -s unlimited. Poznamenejme zde, že pro ladění nehotových programů je vhodné místo optimalizace zapnout raději rozšířenou kontrolu chyb, a to volbami např.
# compiler flags for debugging
gfortran -fcheck=all -std=f2018
ifort -check all -std18 (or -C instead of -check all)
nvfortran -Mbounds -Mstandard (or -C instead of -Mbounds)
Ve Windows lze použít BLAS a LAPACK v implementacích OpenBLAS a Intel MKL, obojí jak s gfortranem, tak s Intel překladači:
:: compiler calls in Windows
:: set LIBRARY_PATH and PATH for OpenBLAS
:: call setvars.bat for Intel compilers or Intel MKL
gfortran -O testMM.f90 -lopenblas
gfortran -O testMM.f90 "%MKLROOT%"\lib\intel64\mkl_rt.lib
ifort -heap-arrays -o a testMM.f90 -Qmkl
ifort -heap-arrays -o a testMM.f90 C:\openblas\lib\libopenblas.lib
Přeložený OpenBLAS pro Windows je ke stažení na github.com jako OpenBLAS-*-x64.zip. Rozbalíme-li obsah archivu např. do adresáře C:\openblas, můžeme pomocí systempropertiesadvanced přidat do proměnné prostředí PATH adresář C:\openblas\bin a do proměnné LIBRARY_PATH uložit C:\openblas\lib, nebo lze všechny soubory z adresářů bin a lib přikopírovat ke zdrojovému souboru a překládat s volbami -L. -lopenblas. Intel MKL je součástí rozměrného Intel oneAPI Base Toolkitu. ifort a ifx ve Windows překládají do exe souboru téhož jména jako první zdrojový soubor; dosáhnout vytvoření souboru a.exe, jak je defaultem u jiných překladačů, lze volbou -o a.
Výkon fortranské funkce matmul je vyšší než rychlost maticového násobení zapsaného pomocí cyklů, a její výkon může dále zvýšit externí nahrazení z rychlého BLASu. Toho lze dosáhnout volbami překladačů (gfortran, resp. ifort) -fexternal-blas, resp. -qopt-matmul, zde včetně voleb -march=native, resp. -xHost pro maximální výkon na překládajícím hardwaru (více pro např. gfortran pod hesly GCC command options a x86 options):
# linking to external matmul
gfortran -fexternal-blas -Ofast -march=native -fmax-stack-var-size=0 testMM.f90 -lblas
ifort -qopt-matmul -Ofast -xHost -heap-arrays testMM.f90 -qmkl
Při honbě po maximálním výkonu hraje roli i možnost paralelizace, kterou knihovny nabízejí. BLAS bývá vláknově paralelizovaný pomocí knihovny OpenMP. Pokud není paralelizace aktivní v defaultu, lze počet vláken nastavit před spuštěním programu pomocí proměnných prostředí, např.
# variables for OpenMP parallelization
export OMP_NUM_THREADS=8 (for BLAS in Linux)
export MKL_NUM_THREADS=8 (for MKL in Linux)
set OMP_NUM_THREADS=8 (for BLAS in Windows)
set MKL_NUM_THREADS=8 (for MKL in Windows)
Drobnost: paralelizaci matmul linkovaného k MKL řídí proměnná OMP_NUM_THREADS.
Rychlost maticového násobení pomocí tří vnořených cyklů lze ovlivnit nastavením úrovně optimalizace. Intel Fortran nabízí navíc autoparalelizaci, s jejíž pomocí se práce vnějšího cyklu rozdělí mezi softwarová vlákna, o které se postarají hardwarová jádra. Lze vyzkoušet:
# auto-parallelization
ifort -parallel -Ofast -xHost -heap-arrays testMM.f90 -qmkl (check OMP_NUM_THREADS)
Zprovoznili jsme i volání knihovny cuBLAS (CUDA BLAS) z balíku CUDA Toolkit, nabízející BLAS pro GPU. Po nezbytných inicializacích prostředí (na našem klastru: source load_nv, podrobněji v záhlaví zdrojového kódu) se překládá, spolu s knihovnou CUDA Runtime, takto:
# cuBLAS for GPU
gfortran testcuBLAS.f90 -lcudart -lcublas
ifort -heap-arrays testcuBLAS.f90 -lcudart -lcublas (Linux)
ifort -heap-arrays testcuBLAS.f90 cudart.lib cublas.lib (Windows)
nvfortran testcuBLAS.f90 -cudalib
CUDA Toolkit nabízí také knihovnu NVBLAS, s jejíž pomocí lze některá volání BLASu (podprogramy *gemm k nim patří) převést na volání cuBLASu bez úpravy zdrojového kódu, jen linkováním; gfortran s volbou -fexternal-blas to dokáže i pro matmul. S povinným souborem nvblas.conf obsahujícím řádek NVBLAS_CPU_BLAS_LIB /opt/nvidia-20.11/hpc_sdk/Linux_x86_64/20.11/compilers/lib/libblas.so (adresa CPU BLASu se nastaví podle potřeby) můžeme překládat takto:
# Routing ordinary BLAS calls to GPU
gfortran testMM.f90 -lnvblas (routing sgemm/dgemm to GPU)
gfortran -fexternal-blas testMM.f90 -lnvblas (routing matmul to GPU)
ifort -heap-arrays testMM.f90 -lnvblas
nvfortran testMM.f90 -cudalib
Benchmark¶
Jsme u domácího úkolu. Přeložte zdrojové kódy z přednášky tolika překladači na tolika operačních systémech, na kolika dokážete. Dobu běhu pak porovnejte s toutéž akcí v Matlabu nebo Octavu,
% MATLAB and GNU Octave
tic; c=a*b; toc
a pythonovským ekvivalentem z NumPy,
# NumPy
t=time.time(); c=numpy.matmul(a,b); print(time.time()-t)
Fortranské zdrojové kódy z přednášky:
maticové násobení testMM.f90 pomocí BLASu, funkce
matmula tří vnořených cyklů,maticové násobení testcuBLAS.f90 pomocí BLASem inspirované funkce z knihovny cuBLAS a
modul mTicToc.f90 s měřičem času v duchu matlabovského
tic; toc;ještě jednou vše v zip souboru.
Pro jednotlivé varianty maticového násobení jsme připravili samostatné moduly, které lze z ukázek vyloupnout a vložit do vlastních programů. Moduly definují a do akce nabízejí fortranské operátory pracující se single/double precision operandy:
! using defined operators for matrix multiplication
use mBLAS ! get access to the modules
use mMatmul
use mLoops
use mcuBLAS
integer,parameter :: sp=4,dp=8
real(dp),allocatable :: a(:,:),b(:,:),c(:,:) ! double-precision allocatable arrays
allocate (a(mmax,kmax),b(kmax,nmax),c(mmax,nmax))
c=a.BLASgemm.b ! via BLAS subroutines sgemm/dgemm
c=real(a,sp).BLASgemm.real(b,sp) ! single-precision operands
c=a.matmul.b ! via matmul function
c=a.mmloops.b ! via indexed do loops
c=a.cuBLASgemm.b ! via BLAS implementation for GPUs
K měření času jsme na přednášce použili měřič založený na fortranském podprogramu system_clock s rozhraním inspirovaným matlabovským tic-toc:
! using wallclock timer
use mTicToc
call tic(); call task(); print *,toc() ! printing the elapsed time
call tic(); call task(); dt=toc() ! saving the time