Librerías secuenciales de
Álgebra Lineal Densa: BLAS y
LAPACK
Domingo Giménez
Departamento de Informática y Sistemas
Universidad de Murcia, Spain
dis.um.es/~domingo
07 de octubre de
2015
Universidad Politécnica de
Valencia
1
Contenido






Jerarquía de librerías
Obteniendo información
Algoritmos por bloques
BLAS
LAPACK
Otras librerías
07 de octubre de 2015
Universidad Politécnica de Valencia
2
Jerarquía de librerías
ScaLAPACK
Paso de mensajes
Direccionamiento
global
PBLAS
Independiente de la
plataforma
LAPACK
BLACS
Secuencial
Direccionamiento
local
BLAS
07 de octubre de 2015
Dependiente de la
plataforma
Comunicaciones: PVM, MPI
Universidad Politécnica de Valencia
3
Jerarquía de librerías: historia
ScaLAPACK 1994
PRISM
LAPACK2
Templates 1995
PBLAS 1994
Templates
LAPACK 1987
BLACS 1991
LINPACK 1979
EISPACK 1974
07 de octubre de 2015
BLAS1 1988
BLAS3 1990
Comunicaciones: PVM, MPI
1993
Universidad Politécnica de Valencia
Rutinas de
comunicación
específicas
4
Jerarquía de librerías:
extensión
PDE Solver
Se puede extender
la jerarquía resolviendo
problemas de alto coste
computacional. Necesarios
algoritmos eficientes
en sistemas de altas
prestaciones.
Least Square Problem
ScaLAPACK
PBLAS
LAPACK
BLAS
07 de octubre de 2015
Inverse Eigenvalue Problem
Universidad Politécnica de Valencia
BLACS
Comunicaciones
5
Obteniendo información


www.netlib.org/liblist.html
www.netlib.org/utk/people/JackDongarr
a/la-sw.html
07 de octubre de 2015
Universidad Politécnica de Valencia
6
Algoritmos por bloques


En vez de realizar operaciones elemento a
elemento realizarlas con bloques de
elementos: menos accesos a memoria para el
mismo volumen de computación  menor
tiempo de ejecución.
Técnica utilizada desde los años 80. Se utiliza
en LAPACK para obtener rutinas eficientes
independientemente del sistema donde se
ejecuten.
07 de octubre de 2015
Universidad Politécnica de Valencia
7
Algoritmos por bloques

Multiplicación de matrices (en SUN Ultra 1):
Método\tamaño
Normal
Traspuesta
Bloques
10
25
50
Bloq tras 10
25
50
Almac blo 10
25
50
Bl tr al bl 10
25
50
Bloq dob 20 5
20 10
50 10
50 25
07 de octubre de 2015
200
0.2179
0.2013
0.2880
0.2192
0.2161
0.2937
0.2195
0.2152
0.2949
0.2277
0.2296
0.2925
0.2244
0.2231
0.6105
0.3206
0.3039
0.2370
400
13.4601
3.3653
2.5901
1.8347
1.7709
2.5026
1.8009
1.7628
2.5122
1.8490
1.8429
2.4985
1.8082
1.7147
4.9363
2.6669
2.4542
1.9221
Universidad Politécnica de Valencia
800
217.5464
27.9945
21.9029
14.9642
14.2502
20.4405
14.6415
14.1806
20.3762
14.8625
14.7314
20.1975
14.5282
13.6553
39.9594
19.7044
19.7044
15.5190
8
Algoritmos por bloques

Las prestaciones de los distintos
algoritmos varían con:





Tamaño del problema
Parámetros del algoritmo
Sistema en que se ejecuta
Forma de la matriz
…
07 de octubre de 2015
Universidad Politécnica de Valencia
9
Algoritmos por bloques

Necesario:



Aprender a usar librerías eficientes
Aprender a desarrollar algoritmos por bloques
Desarrollar técnicas de autooptimización, que
seleccionen el algoritmo y los parámetros a usar,
para obtener buenas prestaciones independiente
de:



El sistema donde se ejecuta
Las condiciones actuales del sistema
Los conocimientos del usuario
07 de octubre de 2015
Universidad Politécnica de Valencia
10
Algoritmos por bloques
mult(a,fa,ca,lda,b,fb,cb,ldb,c,fc,cc,ldc)
double *a; int fa,ca,lda; double *b; int fb,cb,ldb; double *c; int fc,cc,ldc;
{
Algoritmo sin bloques (normal).
int i,j,k; double s;
Acceso elemento a elemento.
Problemas pequeños buenas
prestaciones pues caben en
memoria de niveles bajos de
la jerarquía.
Problemas grandes peores
prestaciones.
for(i=0;i<fa;i++)
for(j=0;j<cb;j++) {
s=0.;
for(k=0;k<ca;k++)
s+=a[i*lda+k]*b[k*ldb+j];
c[i*ldc+j]=s;
}
}
07 de octubre de 2015
Universidad Politécnica de Valencia
11
Algoritmos por bloques
multbloques(a,fa,ca,lda,b,fb,cb,ldb,c,fc,cc,ldc,tb)
double *a; int fa,ca,lda; double *b; int fb,cb,ldb; double *c; int fc,cc,ldc, tb;
{
Algoritmo por bloques.
int i,j,k; double *s;
s=(double *) malloc(sizeof(double)* tb * tb); Acceso y operaciones por
bloques .
Buenas prestaciones
independiente del tamaño.
El tamaño de bloque es
parámetro a determinar.
for(i=0;i<fa;i=i+ tb)
for(j=0;j<cb;j=j+ tb) {
ceros(s, tb, tb, tb);
for(k=0;k<ca;k=k+ tb)
multsumar(&a[i*lda+k], tb, tb,lda,&b[k*ldb+j], tb, tb,ldb,s, tb, tb, tb);
copiar(s, tb, tb, tb,&c[i*ldc+j], tb, tb,ldc);
}
free(s);
}
07 de octubre de 2015
Universidad Politécnica de Valencia
12
Algoritmos por bloques
A
i
tb
tb
B
k
k
j
C
tb
j
i
s
07 de octubre de 2015
Universidad Politécnica de Valencia
13
Algoritmos por bloques
multbloquesgrandes(a,fa,ca,lda,b,fb,cb,ldb,c,fc,cc,ldc,tb,tbp)
double *a; int fa,ca,lda; double *b; int fb,cb,ldb; double *c; int fc,cc,ldc, tb, tbp;
{
Algoritmo por bloques dobles.
int i,j,k; double *s;
La operación sobre bloques
s=(double *) malloc(sizeof(double)* tb * tb);
no es la multiplicación
directa, sino por bloques.
Tenemos dos tamaños
de bloque.
for(i=0;i<fa;i=i+ tb)
for(j=0;j<cb;j=j+ tb) {
ceros(s, tb, tb, tb);
for(k=0;k<ca;k=k+ tb)
multsumarbloques(&a[i*lda+k], tb, tb,lda,&b[k*ldb+j], tb, tb,ldb,s, tb, tb, tb, tbp);
copiar(s, tb, tb, tb,&c[i*ldc+j], tb, tb,ldc);
}
free(s);
}
07 de octubre de 2015
Universidad Politécnica de Valencia
14
Algoritmos por bloques
Almacenamiento por bloques:
almacenamiento
matriz
0
1
2
3
0
1
4
5
4
5
6
7
2
3
6
7
8
9 10
11
8
9 12
12 13 14 15
13
10 11 14 15
posible acceso más rápido a los datos dentro de las operaciones
por bloques
07 de octubre de 2015
Universidad Politécnica de Valencia
15
Algoritmos por bloques

Trasposición de la matriz B y multiplicación:
mult(a,fa,ca,lda,b,fb,cb,ldb,c,fc,cc,ldc)
double *a; int fa,ca,lda; double *b; int fb,cb,ldb; double *c; int fc,cc,ldc;
{
int i,j,k; double s;
for(i=0;i<fa;i++)
for(j=0;j<cb;j++) {
s=0.;
for(k=0;k<ca;k++)
s+=a[i*lda+k]*b[j*ldb+k];
c[i*ldc+j]=s;
}
}
07 de octubre de 2015
Universidad Politécnica de Valencia
16
BLAS
ScaLAPACK
Paso de mensajes
Direccionamiento
global
PBLAS
Independiente de la
plataforma
LAPACK
Basic Linear
Algebra
Subprograms
07 de octubre de 2015
BLACS
Direccionamiento
local
Secuencial
BLAS
Dependiente de la
plataforma
Comunicaciones: PVM, MPI
Universidad Politécnica de Valencia
17
BLAS

Conjunto de rutinas para realizar operaciones
básicas sobre vectores y matrices

1.
2.
3.
Publications/references for the BLAS?
C. L. Lawson, R. J. Hanson, D. Kincaid, and F. T. Krogh, Basic Linear
Algebra Subprograms for FORTRAN usage, ACM Trans. Math. Soft., 5
(1979), pp. 308--323.
J. J. Dongarra, J. Du Croz, S. Hammarling, and R. J. Hanson, An
extended set of FORTRAN Basic Linear Algebra Subprograms, ACM
Trans. Math. Soft., 14 (1988), pp. 1--17.
J. J. Dongarra, J. Du Croz, I. S. Duff, and S. Hammarling, A set of Level
3 Basic Linear Algebra Subprograms, ACM Trans. Math. Soft., 16
(1990), pp. 1--17.
07 de octubre de 2015
Universidad Politécnica de Valencia
18
BLAS

Hay tres niveles según el coste
computacional:
tipo
operaciones
coste
computacional
accesos
memoria
BLAS1 vector-vector
n
n
BLAS2 matriz-vector
n2
n2
BLAS3 matriz-matriz
n3
n2
07 de octubre de 2015
Universidad Politécnica de Valencia
19
BLAS 1
07 de octubre de 2015
Universidad Politécnica de Valencia
20
BLAS 1
Ejemplo ddot.f
Calcula el producto escalar de dos
vectores
Se puede usar en el bucle más interno de
la multiplicación de matrices, dando
lugar a una versión con BLAS 1
Se compila con
icc –O3 mb1.c –lgslcblas -lm
07 de octubre de 2015
Universidad Politécnica de Valencia
21
BLAS

X:
Formato de las funciones (niveles 2 y 3): XYYZZZ
Tipo de datos:
S : REAL
D : DOUBLE PRECISION
C : COMPLEX
Z : DOUBLE COMPLEX
YY:
Tipo de matriz:
GE, GB, HE, HP, HB, SY, SP, TR, TP, TB
ZZZ:
Operación:
MV: productor matriz vector
MM: producto matriz matriz
SV: sistema de ecuaciones ...
07 de octubre de 2015
Universidad Politécnica de Valencia
22
BLAS 2
07 de octubre de 2015
Universidad Politécnica de Valencia
23
BLAS 2
07 de octubre de 2015
Universidad Politécnica de Valencia
24
BLAS 2
Ejemplo dgemv.f
Calcula el producto de una matriz por un
vector
Se puede usar en el segundo bucle en la
multiplicación de matrices, dando lugar
a una versión con BLAS 2
Se compila con
icc –O3 mb2.c –lgslcblas -lm
07 de octubre de 2015
Universidad Politécnica de Valencia
25
BLAS 3
07 de octubre de 2015
Universidad Politécnica de Valencia
26
BLAS 3
Ejemplo dgemm.f
Calcula el producto de una matriz por un
vector
Se puede hacer la multiplicación de
matrices llamando directamente a la
rutina correspondiente de BLAS
Se compila con
icc –O3 mb3.c –lgslcblas -lm
07 de octubre de 2015
Universidad Politécnica de Valencia
27
BLAS

Multiplicación de matrices (en kefren,
pentium 4):
Método\tamaño
Normal
Traspuesta
200
0.0463
0.0231
400
0.7854
0.2875
800
7.9686
2.3190
Bloques 10
25
50
Bloq dob 20 5
20 10
50 10
50 25
Blas 1
Blas 2
Blas 3
0.0255
0.0265
0.0219
0.0393
0.0269
0.0316
0.0215
0.0536
0.0501
0.0429
0.2493
0.2033
0.1785
0.3669
0.3090
0.2232
0.1755
0.8190
0.5861
0.6115
2.0327
1.6928
1.6594
3.4955
2.4424
2.2768
1.4726
8.2311
5.9997
4.7252
07 de octubre de 2015
Universidad Politécnica de Valencia
28
BLAS - Práctica

Hacer llamada a dgemm desde
programa en C y en Fortran:


Ver las formas diferentes en la llamada
Comparar los resultados, teniendo en
cuenta que en C los elementos se
almacenan por filas y en Fortran por
columnas (en la multiplicación de matrices
estas no deben ser simétricas)
07 de octubre de 2015
Universidad Politécnica de Valencia
29
LAPACK
ScaLAPACK
Paso de mensajes
Direccionamiento
global
PBLAS
Independiente de la
plataforma
Linear
Algebra
Package
LAPACK
Dependiente de la
plataforma
Direccionamiento
local
Secuencial
BLAS
07 de octubre de 2015
BLACS
Comunicaciones: PVM, MPI
Universidad Politécnica de Valencia
30
LAPACK

Conjunto de rutinas para resolver problemas
de los más frecuentes en álgebra lineal
densa: sistemas de ecuaciones y problemas
de valores propios

1.
Documentos:
Implementation Guide for LAPACK
UT, CS-90-101, April 1990.
E. Anderson and J. Dongarra
2.
LAPACK: A Portable Linear Algebra Library for HighPerformance Computers
UT, CS-90-105, May 1990.
E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. DuCroz, A.
Greenbaum, S. Hammarling, A. McKenney, D. Sorensen
07 de octubre de 2015
Universidad Politécnica de Valencia
31
LAPACK
Algoritmos orientados a bloques


Basados en BLAS
Eficiencia
Portabilidad
07 de octubre de 2015
Universidad Politécnica de Valencia
32
LAPACK

Problemas que resuelve:





Sistemas de ecuaciones lineales
Problemas de mínimos cuadrados
Problemas de valores propios
Problemas de valores singulares
Otros: factorización de matrices,
estimación del número de condición, etc.
07 de octubre de 2015
Universidad Politécnica de Valencia
33
LAPACK

Tipos de rutinas:

“Driver routines” – Rutinas conductoras.


“Computational routines” – Rutinas
computacionales.


Resuelve un problema.
Realizan una tarea computacional
“Auxiliary routines” – Rutinas auxiliares.

Realizan una subtarea o trabajo de menor
nivel.
07 de octubre de 2015
Universidad Politécnica de Valencia
34
LAPACK

Tipos de matrices:



Densas.
Banda.
Reales y complejas.
… no escasas

Tipos de sistemas:


Secuenciales.
Memoria compartida.
07 de octubre de 2015
Universidad Politécnica de Valencia
35
LAPACK

Rutinas conductoras:

Para la resolución completa de problemas
estándar:



Sistemas de ecuaciones lineales.
Problemas de valores propios.
Siempre que sea posible es recomendable
usar estas rutinas para resolver un
problema.
07 de octubre de 2015
Universidad Politécnica de Valencia
36
LAPACK

Rutinas computacionales:

Realizan tareas computacionales:



Factorizaciones LU y QR, reducción de matriz
simétrica a tridiagonal, ...
Cada rutina conductora realiza una
secuencia de llamadas a las rutinas
computacionales.
El usuario también puede llamar en sus
programas a rutinas computacionales.
07 de octubre de 2015
Universidad Politécnica de Valencia
37
LAPACK

Rutinas auxiliares:

Son rutinas que hacen operaciones de bajo
nivel:



Versiones no orientadas a bloques de
algoritmos orientados a bloques.
Computaciones de bajo nivel (escalar una
matriz, generación de matriz de Householder).
Extensiones de BLAS.
07 de octubre de 2015
Universidad Politécnica de Valencia
38
LAPACK
Formato de rutinas conductoras y
computacionales: XYYZZZ
X: Tipo de datos:

S : REAL
D : DOUBLE PRECISION
C : COMPLEX
Z : DOUBLE COMPLEX
YY: Tipo de matriz
ZZZ: Operación:
SV: sistemas de ecuaciones
EV: valores propios ...
07 de octubre de 2015
Universidad Politécnica de Valencia
39
LAPACK

Tipos de matrices:
BD bidiagonal; GB general band; GE general (i.e.,
unsymmetric, in some cases rectangular); GG general matrices,
generalized problem (i.e., a pair of general matrices); GT
general tridiagonal; HB (complex) Hermitian band; HE
(complex) Hermitian; HG upper Hessenberg matrix, generalized
problem (i.e a Hessenberg and a triangular matrix); HP
(complex) Hermitian, packed storage; HS upper Hessenberg; OP
(real) orthogonal, packed storage; OR (real) orthogonal; PB
symmetric or Hermitian positive definite band; PO symmetric
or Hermitian positive definite; PP symmetric or Hermitian
positive definite, packed storage; PT symmetric or Hermitian
positive definite tridiagonal; SB (real) symmetric band; SP
symmetric, packed storage; ST (real) symmetric tridiagonal;
SY symmetric; TB triangular band; TG triangular matrices,
generalized problem (i.e., a pair of triangular matrices); TP
triangular, packed storage; TR triangular (or in some cases
quasi-triangular); TZ trapezoidal; UN (complex) unitary; UP
(complex) unitary, packed storage
07 de octubre de 2015
Universidad Politécnica de Valencia
40
LAPACK

Ecuaciones lineales: AX = B

Rutina simple: xyySV


Factoriza A y sobreescribe B con X
Rutina experta: xyySVX. Puede llevar a cabo
otras funciones:




ATX=B o AHX=B
Número de condición, singularidad, ...
Refina la solución y hace análisis de error.
Equilibrado del sistema.
07 de octubre de 2015
Universidad Politécnica de Valencia
41
LAPACK
Ejemplo dgesv.f
Resuelve un sistema de ecuaciones
Llamada en Fortran:

call dgesv( )
En C:
dgesv_( )
y se pasan las referencias a los parámetros
07 de octubre de 2015
Universidad Politécnica de Valencia
42
LAPACK - Prácticas
Resolver un sistema de ecuaciones
usando dgesv.f
Utilizando C y Fortran, y comparar las dos
formas de llamada

07 de octubre de 2015
Universidad Politécnica de Valencia
43
LAPACK


Problema de mínimos cuadrados

Problema:

Rutinas: xyyLS, xyyLSX, xyyLSS
Problema generalizado de mínimos
cuadrados


Problema:
Rutinas: xyyLSE, xyyGLM
07 de octubre de 2015
Universidad Politécnica de Valencia
44
LAPACK

Problema simétrico de valores propios

Problema:

Rutinas:



Conductora simple: xyyEV
Conductora divide y vencerás: xyyEVD
Conductora experta: xyyEVX
07 de octubre de 2015
Universidad Politécnica de Valencia
45
LAPACK

También:





Valores propios no simétrico.
Descomposición en valores singulares.
Valores propios simétrico generalizado.
Valores propios no simétrico generalizado.
Descomposición en valores singulares
generalizado.
07 de octubre de 2015
Universidad Politécnica de Valencia
46
Algoritmos por bloques Práctica

Factorización LU por bloques:
A11 A12 A13
A21 A22 A23
=
A31 A32 A33
Paso
Paso
Paso
Paso
1:
2:
3:
4:
U11 U12 U13
L11
L21 L22
L31 L32 L33
U22 U23
U33
A11  L11U11 (factorización LU sin bloques)
A1i  L11U1i
(sistema triangular inferior múltiple)
Ai1  Li1U11
(sistema triangular superior múltiple)
Aij  Aij  Li1U1 j (actualización de bloque sur-este)
07 de octubre de 2015
Universidad Politécnica de Valencia
47
Algoritmos por bloques Práctica



Programar la factorización LU por bloques
utilizando llamadas a rutinas de BLAS y
LAPACK. Usar la LU sin bloques como
referencia. Hacerla sin pivotamiento.
Comparar el tiempo de ejecución con el
obtenido con la rutina sin bloques.
Comparar el tiempo de ejecución
sustituyendo la multiplicación de matrices
por la versión por bloques programada.
07 de octubre de 2015
Universidad Politécnica de Valencia
48
ATLAS
ScaLAPACK
Paso de mensajes
Direccionamiento
global
PBLAS
Independiente de la
plataforma
Automatic
Tuned Linear
Algebra
Software
07 de octubre de 2015
LAPACK
BLACS
Direccionamiento
local
Secuencial
ATLAS
Dependiente de la
plataforma
Comunicaciones: PVM, MPI
Universidad Politécnica de Valencia
49
ATLAS

Rutinas básicas de Álgebra Lineal con
autooptimización.




En la instalación se realizan una serie de tests para decidir
los mejores parámetros a utilizar en la resolución de los
problemas, dependiendo de la rutina y del tamaño de los
vectores y matrices.
Tiempo de instalación grande.
Se tendría disponible software eficiente para nuevas
plataformas sin esperar el tedioso proceso de desarrollar
software optimizado.
La misma funcionalidad que BLAS.
Se puede hacer interface para que las llamadas a las
sean iguales.Universidad Politécnica de Valencia
07 derutinas
octubre de 2015

50
FLAME
ScaLAPACK
Paso de mensajes
Direccionamiento
global
PBLAS
Independiente de la
plataforma
LAPACK
Formal Linear
Algebra Methods
Environment
FLAME
07 de octubre de 2015
BLACS
Dependiente de la
plataforma
Direccionamiento
local
Secuencial
Comunicaciones: PVM, MPI
Universidad Politécnica de Valencia
51
FLAME




Desarrollo de rutinas de Álgebra Lineal
manera sistemática.
Metodología orientada a objeto.
Código más sencillo.
Para algunos tamaños (matrices
rectangulares) mejores prestaciones que
ATLAS y BLAS.
07 de octubre de 2015
Universidad Politécnica de Valencia
52
Descargar

Designing Polylibreries to Speed Up Linear Algebra