PDF de programación - Algoritmos matriciales por bloques - Computación Matricial y Paralela

Master en Nuevas Tecnologías en Informática

Programación Paralela y Computación de Altas Prestaciones



Algoritmos Matriciales por Bloques



Javier Cuenca & Domingo Giménez

Grupo de Investigación en Computación Científica y Programación Paralela



Universidad de Murcia

1

Contenido

 Trabajo por bloques

 Multiplicación de matrices

 Factorización LU

 Optimización automática: Tamaño de bloque óptimo

 Trabajo de los alumnos. Uso de colas en HeteroSolar

 Trabajo de los alumnos. Ejercicios



2

Trabajo por bloques

 En las operaciones anteriores los costes son:



 Vector-vector

 Matriz-vector

 Matriz-matriz



Coste Computacional



n

n2

n3



Memoria



n

n2

n2

 Desarrollando algoritmos con operaciones matriz-matriz, para el mismo número de

operaciones aritméticas menos accesos a memoria  menor tiempo de ejecución



 Usado en el desarrollo de librerías desde los 80 (BLAS-3, LAPACK)



 Posibilidad de migración a rutinas paralelas más escalables



3

Trabajo por bloques

 La reducción varía de un sistema a otro

 ¿Cómo se sabe el tamaño de bloque óptimo? Varía con:

 Sistema

 Problema a resolver

 Tamaño del problema



 Con lo que el método preferido también varía con el tamaño y el sistema

(polialgoritmos)



 Difícil determinar a priori mejor método y parámetros  métodos de

optimización automática



4

Contenido

 Trabajo por bloques

 Multiplicación de matrices

 Factorización LU

 Optimización automática: Tamaño de bloque óptimo

 Trabajo de los alumnos. Uso de colas en HeteroSolar

 Trabajo de los alumnos. Ejercicios



5

Multiplicación de matrices

void matriz_matriz_ld (double *a,int fa,int ca,int lda,

double *b,int fb,int cb,int ldb,double *c,int fc,int

cc,int ldc)

{ int i,j,k;



double s;

for(i=0;i<fa;i++)

for(j=0;j<cb;j++)



{



s=0.;

Algoritmo sin bloques (normal).
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(k=0;k<ca;k++)



s+=a[i*lda+k]*b[k*ldb+j];



c[i*ldc+j]=s;

}

6

Multiplicación de matrices

A

tb

i

k

B

tb

j

C

tb

j

k

i

s

7

Multiplicación de matrices

void matriz_matriz_bloques (double *a,int fa,int ca,int lda,double
*b,int fb,int cb,int ldb,double *c,int fc,int cc,int ldc,int tb)

{ int i,j,k; double *s;



s=(double *) malloc(sizeof(double)* tb * tb);

for(i=0;i<fa;i=i+ tb)



for(j=0;j<cb;j=j+ tb)

{



matriz_cero(s, tb, tb, tb);

for(k=0;k<ca;k=k+ tb)

Algoritmo por bloques.
Acceso y operaciones por
bloques .
Buenas prestaciones
independiente del tamaño.
El tamaño de bloque es
parámetro a determinar.



}

free(s);

}

multiplicar_acumular_matrices(&a[i*lda+k], tb,
tb,lda,&b[k*ldb+j],tb, tb,ldb,s, tb, tb, tb);



copiar_matriz(s,tb, tb, tb,&c[i*ldc+j], tb, tb,ldc);

8

void matriz_matriz_bloques (double *a,int fa,int ca,int lda,double
*b,int fb,int cb,int ldb,double *c,int fc,int cc,int ldc,int tb)

Multiplicación de matrices

void multiplicar_acumular_matrices(double *a,int fa,int ca,int lda,double *b,int
fb,int cb,int ldb,double *c,int fc,int cc,int ldc)
{
int i,j,k,kb;
double *da,*db,s;

for(i=0;i<fa;i++)
{
da=&a[i*lda];
for(j=0;j<cb;j++)
{
db=&b[j];
s=c[i*ldc+j];
for(k=0,kb=0;k<ca;k++,kb=kb+ldb)
{
s=s+da[k]*db[kb];
}

c[i*ldc+j]=s;
}
}
}

Algoritmo por bloques.
Acceso y operaciones por
bloques .
Buenas prestaciones
independiente del tamaño.
El tamaño de bloque es
parámetro a determinar.

s=(double *) malloc(sizeof(double)* tb * tb);

{ int i,j,k; double *s;

matriz_cero(s, tb, tb, tb);

for(i=0;i<fa;i=i+ tb)

for(j=0;j<cb;j=j+ tb)



{



for(k=0;k<ca;k=k+ tb)



}

free(s);

}

multiplicar_acumular_matrices(&a[i*lda+k], tb,
tb,lda,&b[k*ldb+j],tb, tb,ldb,s, tb, tb, tb);



copiar_matriz(s,tb, tb, tb,&c[i*ldc+j], tb, tb,ldc);

9

Multiplicación de matrices

void matriz_matriz_bloquesdobles (double *a,int fa,int ca,int

lda,double *b,int fb,int cb,int ldb,double *c,int fc,int cc,int
ldc,int tb,int tbp)

{ int i,j,k; double *s;



s=(double *) malloc(sizeof(double)* tb * tb);

for(i=0;i<fa;i=i+ tb)



for(j=0;j<cb;j=j+ tb)

{



matriz_cero(s, tb, tb, tb);

for(k=0;k<ca;k=k+ tb)

Algoritmo por bloques dobles.
La operación sobre bloques
no es la multiplicación
directa, sino por bloques.
Tenemos dos tamaños
de bloque.



multiplicar_acumular_bloques(&a[i*lda+k],tb,


tb,lda,&b[k*ldb+j],tb,tb,ldb,s,tb,tb, tb, tbp);

copiar_matriz(s, tb, tb, tb,&c[i*ldc+j], tb, tb,ldc);



}

free(s);

}

10

void matriz_matriz_bloquesdobles (double *a,int fa,int ca,int

Multiplicación de matrices


void multiplicar_acumular_bloques(double *a,int fa,int ca,int lda,double *b,int fb,int cb,int ldb,double
*c,int fc,int cc,int ldc,int tbp)
{
int i,j,k;
double *s;

s=(double *) malloc(sizeof(double)*tbp*tbp);

for(i=0;i<fa;i=i+tbp)
{
for(j=0;j<cb;j=j+tbp)
{
copiar_matriz(&c[i*ldc+j],tbp,tbp,ldc,s,tbp,tbp,tbp);
for(k=0;k<ca;k=k+tbp)
{
multiplicar_acumular_matrices(&a[i*lda+k],tbp,tbp,lda,&b[k*ldb+j],tbp,tbp,ldb,s,tbp,tbp,tbp);
}
copiar_matriz(s,tbp,tbp,tbp,&c[i*ldc+j],tbp,tbp,ldc);
}
}
free(s);


Algoritmo por bloques dobles.
La operación sobre bloques
no es la multiplicación
directa, sino por bloques.
Tenemos dos tamaños
de bloque.

lda,double *b,int fb,int cb,int ldb,double *c,int fc,int cc,int
ldc,int tb,int tbp)

s=(double *) malloc(sizeof(double)* tb * tb);

{ int i,j,k; double *s;

for(i=0;i<fa;i=i+ tb)



for(j=0;j<cb;j=j+ tb)



{



matriz_cero(s, tb, tb, tb);

for(k=0;k<ca;k=k+ tb)



multiplicar_acumular_bloques(&a[i*lda+k],tb,


tb,lda,&b[k*ldb+j],tb,tb,ldb,s,tb,tb, tb, tbp);

copiar_matriz(s, tb, tb, tb,&c[i*ldc+j], tb, tb,ldc);



}

free(s);

}

11

Multiplicación de matrices



Almacenamiento por bloques:

matriz

0 1 2 3

4 5 6 7

8 9 10 11

12 13 14 15


almacenamiento

0 1 4 5

2 3 6 7

8 9 12 13

10 11 14 15


posible acceso más rápido a los datos dentro de las operaciones
por bloques

12

Multiplicación de matrices

1000
12.70

 Multiplicación de matrices, en portátil antiguo:
1400
Método\tamaño
normal
36.41


bloques 25
50
100


3.69
3.56
4.25

6.30
5.90
6.33

9.25
8.71
8.95

1200
21.95



Reducción 76%



bloques 25 5
dobles
50 10
50 25


100 20
100 25
100 50



4.67
5.03
4.53
4.87
4.78
3.90



7.87
8.08
7.16
7.33
7.06
5.85



10.89
12.93
11.11
10.97
9.92
8.92

13

Multiplicación de matrices

 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



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



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

Reducción 93%

14

Multiplicación de matrices



0.7854

7.9686



Reducción 79%

400



800

0.0463

0.0231

 Multiplicación de matrices,en kefren, pentium 4:

Método\tamaño 200

Normal

Traspuesta

Bloques 10

25

50

Bloq dob 20 5
20 10


50 10

50 25



0.0393
0.0269
0.0316
0.0215

0.3669
0.3090
0.2232
0.1755

3.4955
2.4424
2.2768
1.4726

0.0255
0.0265
0.0219

0.2493
0.2033
0.1785

2.0327
1.6928
1.6594



0.2875

2.3190



15

Multiplicación de matrices

 Multiplicación de matrices,en PC actual (2016):

Método\tamaño 500
1.02
Normal
Traspuesta
0.60

Bloques
50

100



1500
54.85
17.08

1000
13.62
5.04

20.94
21.08

6.55
6.29

0.87
0.86



16

Contenido

 Trabajo por bloques

 Multiplicación de matrices

 Factorización LU

 Optimización automática: Tamaño de bloque óptimo

 Trabajo de los alumnos. Uso de colas en HeteroSolar

 Trabajo de los alumnos. Ejercicios



17

Factorización LU



Cada Aij, Lij, Uij de tamaño b×n :

Factorización sin bloques

Paso 1: L00 U00=A00
Paso 2: L00 U01=A01
Paso 3: L10 U00=A10
Paso 4: A11 =L10 U01 + L11 U11  A’11 =A11 - L10 U01 , por bloques
y seguir trabajando con el nuevo valor de A’11

18

 Sistema múltiple triangular inferior (¿bloques?)

 Sistema múltiple triangular superi