PDF de programación - Librería secuencial de Álgebra Lineal Densa LAPACK

Librería secuencial de

Álgebra Lineal Densa LAPACK

Domingo Giménez

Javier Cuenca

Facultad de Informática
Universidad de Murcia

1

LAPACK

Independiente de la
plataforma
Linear
Algebra
Package

ScaLAPACK

Paso de mensajes

Direccionamiento
global

PBLAS

LAPACK

BLACS

Secuencial

BLAS

Comunicaciones: PVM, MPI
2

Dependiente de la
plataforma

Direccionamiento
local

LAPACK
 Conjunto de rutinas para resolver problemas
de los más frecuentes en álgebra lineal
densa: sistemas de ecuaciones y problemas
de valores propios
 Documentos:
Implementation Guide for LAPACK

UT, CS-90-101, April 1990. E. Anderson and J. Dongarra

LAPACK: A Portable Linear Algebra Library for High-Performance

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

3

LAPACK
 Algoritmos orientados a bloques

Basados en BLAS
Eficiencia
Portabilidad

4


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.

5

LAPACK
 Tipos de matrices:

 Densas.
 Banda.
 Reales y complejas.
 Tipos de sistemas:

… no escasas (dispersas, vacías...)

 Secuenciales.
 Memoria compartida (versiones
multithreaded).
6

LAPACK
 Tipos de rutinas:

 “Driver routines” – Rutinas conductoras.

 Resuelve un problema.

 “Computational routines” – Rutinas
computacionales.
 Realizan una tarea computacional

 “Auxiliary routines” – Rutinas auxiliares.
 Realizan una subtarea o trabajo de menor
nivel.

7

LAPACK. Tipos de rutinas
 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.

8

LAPACK. Tipos de rutinas
 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.

9

LAPACK. Tipos de rutinas
 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.

10

LAPACK

computacionales: XYYZZZ

 Formato de rutinas conductoras y
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 ...

11

LAPACK

 Tipos de matrices YY (1/2):

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;

12

LAPACK

 Tipos de matrices YY (2/2):

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

13

LAPACK

 Rutinas conductoras de resolución de

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.

14

LAPACK. Ejemplo dgesv
 Ejemplo dgesv

 Resuelve un sistema de ecuaciones
 Llamada en Fortran:

call dgesv( )

 En C:

dgesv_( )
y se pasan las referencias a los
parámetros

15

LAPACK. Ejemplo dgesv
 dgesv

 Rutina conductora de LAPACK
 Resolución de un sistema de ecuaciones AX=B
 Llamadas:

 dgetrf

 dgetrs

 Rutina computacional de LAPACK
 Factorización LU: Transforma A  LU

 Rutina computacional de LAPACK
 Resuelve el doble sistema triangular LU X = B

16

LAPACK. Ejemplo dgesv
 dgetrf

 Rutina computacional de LAPACK
 Factorización LU: Transforma A  LU
 Llamadas en cada pasada de bucle:

 dgetf2

 Rutina auxiliar de LAPACK
 Factorización LU sin bloques aplicada a determinados

bloques de A

 dtrsm (2 veces por pasada)
 Rutina del nivel 3 de BLAS
 Resuelve un sistema triangular de ecuaciones

 dgemm

 Rutina del nivel 3 de BLAS
 Multiplicación de matrices

17

 dgetrs

LAPACK. Ejemplo dgesv
 Rutina computacional de LAPACK
 Resuelve el doble sistema triangular LU X =B
 Llamadas en cada pasada de bucle:

 Rutina auxiliar de LAPACK
 Aplica a B los intercambios de filas realizados previamente a

las matrices L y U

 Rutina del nivel 3 de BLAS
 Resuelve un sistema triangular de ecuaciones LY=B

 dlaswp

 dtrsm

 dtrsm

 Rutina del nivel 3 de BLAS
 Resuelve un sistema triangular de ecuaciones UX=Y

18

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.



19

Factorización LU



Cada Aij, Lij, Uij de tamaño bb
Paso 1: L00 U00=A00  Factorización sin bloques
Paso 2: L00 U01=A01  Sistema múltiple triangular inferior (¿bloques?)
Paso 3: L10 U00=A10  Sistema múltiple triangular superior (¿bloques?)
Paso 4: A11 =L10 U01 + L11 U11  A’ 11 =A11 - L10 U01 , por bloques
y seguir trabajando con el nuevo valor de A11

20

Factorización LU

void lu_bloques(double *a,int fa,int ca,int lda,int tb)
{int i,j,k,f,c;
for(i=0;i<fa;i=i+tb) {
f=(tb<fa-i ? tb : fa-i); c=(tb<ca-i ? tb : ca-i);
lu(&a[i*lda+i],f,c,lda);
if(i+tb<fa) {
sistema_triangular_inferior(&a[i*lda+i],f,c,lda,&a[i*lda+i+c],f,ca-i-c,lda);
sistema_triangular_superior(&a[i*lda+i],f,c,lda,&a[(i+f)*lda+i],fa-i-f,c,lda);

multiplicar_restar_matrices(&a[(i+f)*lda+i],fa-i-f,c,lda,
&a[i*lda+i+f],f,ca-i-c,lda,&a[(i+f)*lda+i+c],fa-i-f,ca-i-c,lda);
} } }

21

Factorización LU

800
2.10
1.42
1.29
1.24
1.20
1.22
1.47
2.29
2.17
1.73

1000
4.01
2.78
2.27
2.37
2.00
2.32
2.24
3.47
3.67
3.43

22

 en mi portátil:
tamaño bloque\matriz
1
12
25
37
44
50
100
200
400

sin bloques

Practicar con los ejemplos

 Comprobar la compilación con MKL en el fichero

mkl.pbs.

 Probar lanzando a otros nodos.
 Comparar prestaciones de la multiplicación de

matrices con las rutinas de BLAS, dependiendo del
nivel de BLAS utilizado y variando el tamaño del
problema y el número de threads

23

Trabajo a entregar

Puntuación sobre 2.
 Hasta 9 de diciembre (1 punto), los experimentos se harán en el cluster y también puede ser en

sistemas propios:

 Comparar prestaciones de multiplicación secuencial por bloques y sin bloques con la

multiplicación de MKL



Identificar en la rutina de factorización LU funciones que se puedan realizar con

llamadas a BLAS, sustituir las llamadas y comparar el tiempo de ejecución.
Comparar los resultados con los de la factorización usando directamente MKL.

 Hasta fecha de examen de la convocatoria (no hay examen), en enero el 24 (1 punto). Hacer dos de

cinco posibles trabajos:

 LU con paralelismo OpenMP+MKL.
 LU con paralelismo MPI+MKL.
 LU combinando MKL en CPU con GPU.
 LU combinando MKL en CPU con Xeon Phi.



Instalación y práctica en multicore con otra librería de álgebra lineal distinta de

MKL.

24