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

Dependiente de la
plataforma

Secuencial

Direccionamiento
local

BLAS

Comunicaciones: PVM, MPI

2

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.

no escasas

 Tipos de sistemas:

 Secuenciales.
 Memoria compartida.

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

 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 ...

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;
12
PP symmetric or Hermitian positive definite, packed storage;

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

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

 dgetrs

 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

LAPACK. Ejemplo dgesv
 dgetrs

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

 dlaswp

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

a las matrices L y U

 dtrsm

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

 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