Política de uso de dados

Ao navegar por este site, você concorda com a política de uso de dados.

Política de uso de dados

Ao navegar por este site, você concorda com a política de uso de dados.

| | | |

Matemática Numérica Paralela

2 Multiprocessamento (MP)2.8 Resolução de Sistema Linear Triangular 2.10 Métodos iterativos para Sistemas Lineares

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

2.9 Decomposição LU

[[badge:remoção]]

Nesta seção, vamos discutir sobre a paralelização da decomposição LU para matrizes. A decomposição $L U$ de uma matriz $A$ com dimensões $n\times n$ é

A=LU

(2.24)

onde $L$ é uma matriz triangular inferior e $U$ é uma matriz triangular superior, ambas com dimensões $n\times n$ .

Para fixar as ideais, vamos denotar $A=[a_{i,j}]_{i,j=0}^{n-1}$ , $L=[l_{i,j}]_{i,j=0}^{n-1}$ sendo $l_{i,i}=1$ e $l_{i,j}=0$ para $i>j$ , e $U=[u_{i,j}]_{i,j=0}^{n}$ sendo $u_{i,j}=0$ para $i<j$ . O pseudoalgoritmo serial para computar a decomposição $L U$ é

1.

U = A, L = I
2.
Para $k=0,\dotsc,n-2$
1. (a)
  Para $i=k+1,\dotsc,n-1$
  1. i.
    
    $l_{i,k}=u_{i,k}/u_{k,k}$
  2. ii.
    
    Para $j=k,\dotsc,n-1$
    
    A.
    
    $u_{i,j}=u_{i,j}-l_{i,k}u_{k,j}$

A forma mais fácil de paralelizar este algoritmo em uma arquitetura MP é paralelizando um de seus laços (itens 2., 2.(a) ou 2.(a)ii.). O laço do item 2. não é paralelizável, pois a iteração seguinte depende do resultado da iteração imediatamente anterior. Agora, os dois laços seguintes são paralelizáveis. Desta forma, o mais eficiente é paralelizarmos o segundo laço 2.(a).

O seguinte código é uma versão paralela da decomposição LU. A matriz $A$ é inicializada como uma matriz simétrica de elementos randômicos (linhas 19-41), sendo que a decomposição é computada nas linhas 43-61.

Código: parallelLU.cc

⬇

1#include <omp.h>

2#include <stdio.h>

3#include <ctime>

4#include <algorithm>

6#include <gsl/gsl_matrix.h>

7#include <gsl/gsl_vector.h>

8#include <gsl/gsl_rng.h>

9#include <gsl/gsl_blas.h>

11int main(int argc, char *argv[]) {

13 int n = 5;

15 gsl_matrix *a = gsl_matrix_alloc(n,n);

16 gsl_matrix *u = gsl_matrix_alloc(n,n);

17 gsl_matrix *l = gsl_matrix_alloc(n,n);

19 // gerador randomico

20 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

21 gsl_rng_set(rng, time(NULL));

23 // inicializacao

24 printf('Inicializando ... \n');

25 for (int i=0; i<n; i++) {

26 for (int j=0; j<i; j++) {

27 int sig = 1;

28 if (gsl_rng_uniform(rng) >= 0.5)

29 sig = -1;

30 gsl_matrix_set(a, i, j,

31 sig*gsl_rng_uniform(rng));

32 gsl_matrix_set(a, j, i,

33 gsl_matrix_get(a, i, j));

34 }

35 int sig = 1;

36 if (gsl_rng_uniform(rng) >= 0.5)

37 sig = -1;

38 gsl_matrix_set(a, i, i,

39 sig*gsl_rng_uniform_pos(rng));

40 }

41 printf('feito.\n');

43 // U = A

44 gsl_matrix_memcpy(u,a);

45 // L = I

46 gsl_matrix_set_identity(l);

48 for (int k=0; k<n-1; k++) {

49 #pragma omp parallel for

50 for (int i=k+1; i<n; i++) {

51 gsl_matrix_set(l, i, k,

52 gsl_matrix_get(u, i, k)/

53 gsl_matrix_get(u, k, k));

54 for (int j=k; j<n; j++) {

55 gsl_matrix_set(u, i, j,

56 gsl_matrix_get(u, i, j) -

57 gsl_matrix_get(l, i, k) *

58 gsl_matrix_get(u, k, j));

59 }

60 }

61 }

63 gsl_matrix_free(a);

64 gsl_matrix_free(u);

65 gsl_matrix_free(l);

66 gsl_rng_free(rng);

68 return 0;

69}

2.9.1 Exercícios Resolvidos

Em remoção

ER 2.9.1.

Faça um código MP para computar a solução de um sistema linear $Ax=b$ usando a decomposição LU. Assuma $A$ uma matriz simétrica $n\times n$ de elementos randômicos, assim como os elementos do vetor $b$ .

Solução 0.

A decomposição LU da matriz $A$ nos fornece as matrizes $L$ (matriz triangular inferior) e $U$ (matriz triangular superior), com

A=LU

(2.25)

Logo, temos

	$\displaystyle Ax=b$		(2.26)
	$\displaystyle\Rightarrow(LU)x=b$		(2.27)
	$\displaystyle\Rightarrow L(Ux)=b$		(2.28)

Denotando $Ux=y$ , temos que $y$ é solução do sistema triangular inferior

Ly=b

(2.29)

e, por conseguinte, $x$ é solução do sistema triangular superior

Ux=y.

(2.30)

Em síntese, o sistema $Ax=b$ pode ser resolvido com o seguinte pseudocódigo:

1.

Computar a decomposição LU, A=LU.
2.

Resolver $Ly=b$ .
3.

Resolver $Ux=b$ .

Cada passo acima pode ser paralelizado. O código MP fica de exercício, veja E 2.9.1.

ER 2.9.2.

Considere a decomposição LU de uma matriz $A$ $n\times n$ . Em muitas aplicações, não há necessidade de guardar a matriz $A$ em memória após a decomposição. Além disso, fixando-se que a diagonal da matriz $L$ tem todos os elementos iguais a 1, podemos alocar seus elementos não nulos na parte triangular inferior (abaixo da diagonal) da própria matriz $A$ . E, a matriz $U$ pode ser alocada na parte triangular superior da matriz $A$ . Faça um código MP para computar a decomposição LU de uma matriz $A$ , alocando o resultado na própria matriz $A$ .

Solução 0.

O seguinte código faz a implementação pedida. Neste código, é necessário alocar apenas a matriz $A$ , sem necessidade de locar as matrizes $L$ e $U$ . Da linha 17 à 39, apenas é gerada a matriz randômica $A$ . A decomposição é computada da linha 41 a 54.

Código: parallelLU2.cc

⬇

1#include <omp.h>

2#include <stdio.h>

3#include <ctime>

4#include <algorithm>

6#include <gsl/gsl_matrix.h>

7#include <gsl/gsl_vector.h>

8#include <gsl/gsl_rng.h>

9#include <gsl/gsl_blas.h>

11int main(int argc, char *argv[]) {

13 int n = 5;

15 gsl_matrix *a = gsl_matrix_alloc(n,n);

17 // gerador randomico

18 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

19 gsl_rng_set(rng, time(NULL));

21 // inicializacao

22 printf('Inicializando ... \n');

23 for (int i=0; i<n; i++) {

24 for (int j=0; j<i; j++) {

25 int sig = 1;

26 if (gsl_rng_uniform(rng) >= 0.5)

27 sig = -1;

28 gsl_matrix_set(a, i, j,

29 sig*gsl_rng_uniform(rng));

30 gsl_matrix_set(a, j, i,

31 gsl_matrix_get(a, i, j));

32 }

33 int sig = 1;

34 if (gsl_rng_uniform(rng) >= 0.5)

35 sig = -1;

36 gsl_matrix_set(a, i, i,

37 sig*gsl_rng_uniform_pos(rng));

38 }

39 printf('feito.\n');

41 for (int k=0; k<n-1; k++) {

42 #pragma omp parallel for

43 for (int i=k+1; i<n; i++) {

44 gsl_matrix_set(a, i, k,

45 gsl_matrix_get(a, i, k)/

46 gsl_matrix_get(a, k, k));

47 for (int j=k+1; j<n; j++) {

48 gsl_matrix_set(a, i, j,

49 gsl_matrix_get(a, i, j) -

50 gsl_matrix_get(a, i, k) *

51 gsl_matrix_get(a, k, j));

52 }

53 }

54 }

55 gsl_matrix_free(a);

56 gsl_rng_free(rng);

58 return 0;

59}

Este algoritmo demanda substancialmente menos memória computacional que o código parallelLU.cc visto acima. Por outro lado, ele é substancialmente mais lento, podendo demandar até o dobro de tempo. Verifique!

O aumento no tempo computacional se deve ao mau uso da memória cache dos processadores. A leitura de um elemento da matriz, aloca no cache uma sequência de elementos próximos na mesma linha. Ao escrever em um destes elementos, a alocação do cache é desperdiçada, forçando o cache a ser atualizado. Note que o código parallelLU.cc requer menos atualizações do cache que o código parallelLU2.cc.

2.9.2 Exercícios

Em remoção

E. 2.9.1.

Implemente o código MP discutido no ER 2.9.1.

E. 2.9.2.

Implemente um código MP para computar a inversa de uma matriz simétrica de elementos randômicos usando decomposição LU.

E. 2.9.3.

Considere o pseudoalgoritmo serial da composição LU apresentado acima. Por que é melhor paralelizar o laço 2.(a) do que o laço o 2.(a)ii.?

E. 2.9.4.

Use o código MP discutido no ER 2.9.2 para resolver um sistema $Ax=b$ de $n$ equações e $n$ incógnitas. Assuma que a matriz $A$ seja simétrica.

E. 2.9.5.

Um algoritmo paralelo mais eficiente para computar a decomposição LU pode ser obtido usando-se a decomposição LU por blocos. Veja o vídeo https://youtu.be/E8aBJsC0bY8 e implemente um código MP para computar a decomposição LU por blocos.

Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Matemática Numérica Paralela

2 Multiprocessamento (MP)2.8 Resolução de Sistema Linear Triangular 2.10 Métodos iterativos para Sistemas Lineares

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

2.9 Decomposição LU

[[badge:remoção]]

Nesta seção, vamos discutir sobre a paralelização da decomposição LU para matrizes. A decomposição $L U$ de uma matriz $A$ com dimensões $n\times n$ é

A=LU

(2.24)

onde $L$ é uma matriz triangular inferior e $U$ é uma matriz triangular superior, ambas com dimensões $n\times n$ .

1.

U = A, L = I
2.
Para $k=0,\dotsc,n-2$
1. (a)
  Para $i=k+1,\dotsc,n-1$
  1. i.
    
    $l_{i,k}=u_{i,k}/u_{k,k}$
  2. ii.
    
    Para $j=k,\dotsc,n-1$
    
    A.
    
    $u_{i,j}=u_{i,j}-l_{i,k}u_{k,j}$

Código: parallelLU.cc

⬇

1#include <omp.h>

2#include <stdio.h>

3#include <ctime>

4#include <algorithm>

6#include <gsl/gsl_matrix.h>

7#include <gsl/gsl_vector.h>

8#include <gsl/gsl_rng.h>

9#include <gsl/gsl_blas.h>

11int main(int argc, char *argv[]) {

13 int n = 5;

15 gsl_matrix *a = gsl_matrix_alloc(n,n);

16 gsl_matrix *u = gsl_matrix_alloc(n,n);

17 gsl_matrix *l = gsl_matrix_alloc(n,n);

19 // gerador randomico

20 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

21 gsl_rng_set(rng, time(NULL));

23 // inicializacao

24 printf('Inicializando ... \n');

25 for (int i=0; i<n; i++) {

26 for (int j=0; j<i; j++) {

27 int sig = 1;

28 if (gsl_rng_uniform(rng) >= 0.5)

29 sig = -1;

30 gsl_matrix_set(a, i, j,

31 sig*gsl_rng_uniform(rng));

32 gsl_matrix_set(a, j, i,

33 gsl_matrix_get(a, i, j));

34 }

35 int sig = 1;

36 if (gsl_rng_uniform(rng) >= 0.5)

37 sig = -1;

38 gsl_matrix_set(a, i, i,

39 sig*gsl_rng_uniform_pos(rng));

40 }

41 printf('feito.\n');

43 // U = A

44 gsl_matrix_memcpy(u,a);

45 // L = I

46 gsl_matrix_set_identity(l);

48 for (int k=0; k<n-1; k++) {

49 #pragma omp parallel for

50 for (int i=k+1; i<n; i++) {

51 gsl_matrix_set(l, i, k,

52 gsl_matrix_get(u, i, k)/

53 gsl_matrix_get(u, k, k));

54 for (int j=k; j<n; j++) {

55 gsl_matrix_set(u, i, j,

56 gsl_matrix_get(u, i, j) -

57 gsl_matrix_get(l, i, k) *

58 gsl_matrix_get(u, k, j));

59 }

60 }

61 }

63 gsl_matrix_free(a);

64 gsl_matrix_free(u);

65 gsl_matrix_free(l);

66 gsl_rng_free(rng);

68 return 0;

69}

2.9.1 Exercícios Resolvidos

Em remoção

ER 2.9.1.

Solução 0.

A decomposição LU da matriz $A$ nos fornece as matrizes $L$ (matriz triangular inferior) e $U$ (matriz triangular superior), com

A=LU

(2.25)

Logo, temos

	$\displaystyle Ax=b$		(2.26)
	$\displaystyle\Rightarrow(LU)x=b$		(2.27)
	$\displaystyle\Rightarrow L(Ux)=b$		(2.28)

Denotando $Ux=y$ , temos que $y$ é solução do sistema triangular inferior

Ly=b

(2.29)

e, por conseguinte, $x$ é solução do sistema triangular superior

Ux=y.

(2.30)

Em síntese, o sistema $Ax=b$ pode ser resolvido com o seguinte pseudocódigo:

1.

Computar a decomposição LU, A=LU.
2.

Resolver $Ly=b$ .
3.

Resolver $Ux=b$ .

Cada passo acima pode ser paralelizado. O código MP fica de exercício, veja E 2.9.1.

ER 2.9.2.

Solução 0.

Código: parallelLU2.cc

⬇

1#include <omp.h>

2#include <stdio.h>

3#include <ctime>

4#include <algorithm>

6#include <gsl/gsl_matrix.h>

7#include <gsl/gsl_vector.h>

8#include <gsl/gsl_rng.h>

9#include <gsl/gsl_blas.h>

11int main(int argc, char *argv[]) {

13 int n = 5;

15 gsl_matrix *a = gsl_matrix_alloc(n,n);

17 // gerador randomico

18 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

19 gsl_rng_set(rng, time(NULL));

21 // inicializacao

22 printf('Inicializando ... \n');

23 for (int i=0; i<n; i++) {

24 for (int j=0; j<i; j++) {

25 int sig = 1;

26 if (gsl_rng_uniform(rng) >= 0.5)

27 sig = -1;

28 gsl_matrix_set(a, i, j,

29 sig*gsl_rng_uniform(rng));

30 gsl_matrix_set(a, j, i,

31 gsl_matrix_get(a, i, j));

32 }

33 int sig = 1;

34 if (gsl_rng_uniform(rng) >= 0.5)

35 sig = -1;

36 gsl_matrix_set(a, i, i,

37 sig*gsl_rng_uniform_pos(rng));

38 }

39 printf('feito.\n');

41 for (int k=0; k<n-1; k++) {

42 #pragma omp parallel for

43 for (int i=k+1; i<n; i++) {

44 gsl_matrix_set(a, i, k,

45 gsl_matrix_get(a, i, k)/

46 gsl_matrix_get(a, k, k));

47 for (int j=k+1; j<n; j++) {

48 gsl_matrix_set(a, i, j,

49 gsl_matrix_get(a, i, j) -

50 gsl_matrix_get(a, i, k) *

51 gsl_matrix_get(a, k, j));

52 }

53 }

54 }

55 gsl_matrix_free(a);

56 gsl_rng_free(rng);

58 return 0;

59}

2.9.2 Exercícios

Em remoção

E. 2.9.1.

Implemente o código MP discutido no ER 2.9.1.

E. 2.9.2.

Implemente um código MP para computar a inversa de uma matriz simétrica de elementos randômicos usando decomposição LU.

E. 2.9.3.

Considere o pseudoalgoritmo serial da composição LU apresentado acima. Por que é melhor paralelizar o laço 2.(a) do que o laço o 2.(a)ii.?

E. 2.9.4.

Use o código MP discutido no ER 2.9.2 para resolver um sistema $Ax=b$ de $n$ equações e $n$ incógnitas. Assuma que a matriz $A$ seja simétrica.

E. 2.9.5.

Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Pedro H A Konzen

| | | |