Colabore! Saiba mais screen_rotation

Ao navegar por este site , você concorda com a Política de Uso de Dados.

| | | | |

2 Multiprocessamento (MP)2.7 Aplicação: Equação de Fisher 2.9 Decomposição LU

2.8 Resolução de Sistema Linear Triangular

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

Em remoção

Nesta seção, vamos discutir sobre a uma implementação em paralelo do método da substituição para a resolução de sistemas triangulares. Primeiramente, vamos considerar $A$ uma matriz triangular inferior quadrada de dimensões $n\times n$ , i.e. $A=[a_{i,j}]_{i,j=0}^{n-1}$ com $a_{i,j}=0$ para $i<j$ . Ainda, vamos coniderar que $A$ é invertível.

Neste caso, um sistema linear $Ax=b$ pode ser escrito na seguinte forma algébrica

	$\displaystyle a_{1,1}x_{1}=b_{1}$		(2.10)
	$\displaystyle\vdots$		(2.11)
	$\displaystyle a_{i,1}x_{1}+a_{i,2}x_{2}+\cdots+a_{i,i-1}x_{i-1}+a_{i,i}x_{i}=b% _{i}$		(2.12)
	$\displaystyle\vdots$		(2.13)
	$\displaystyle a_{n,1}x_{1}+a_{n,2}x_{2}+\cdots+a_{i,i}x_{i}+\cdots+a_{n,n}x_{n% }=b_{n}$		(2.14)

O algoritmo serial do método da substituição (para frente) resolve o sistema começando pelo cálculo de $x_{1}$ na primeira equação, então o cálculo de $x_{2}$ pela segunda equação e assim por diante até o cálculo de $x_{n}$ pela última equação. Segue o pseudocódigo serial.

1.
Para $i=0,\dotsc,n-1$ :
1. (a)
  Para $j=0,\dotsc,i-1$ :
  1. i.
    
    $b_{i}=b_{i}-A_{i,j}x_{j}$
2. (b)
  
  $\displaystyle x_{i}=\frac{b_{i}}{A_{i,i}}$

Implemente!

Para o algoritmo paralelo, vamos considerar uma arquitetura MP com $n_{p}\geq 1$ instâncias de processamento. Para cada instância de processamento $1\leq p_{id}<n_{p}-1$ vamos alocar as seguintes colunas da matriz $A$

	$\displaystyle t_{ini}=p_{id}\left\lfloor\frac{n}{n_{p}}\right\rfloor$		(2.15)
	$\displaystyle t_{fim}=(p_{id}+1)\left\lfloor\frac{n}{n_{p}}\right\rfloor-1$		(2.16)

e, para $p_{id}=n_{p}-1$ vamos alocar as últimas colunas, i.e.

	$\displaystyle t_{ini}=p_{id}\left\lfloor\frac{n}{n_{p}}\right\rfloor$		(2.17)
	$\displaystyle t_{fim}=n-1$		(2.18)

Segue o pseudocódigo em paralelo.

1.
Para $i=0,\dotsc,n-1$
1. (a)
  
  $s=0$
2. (b)
  Região paralela
  1. i.
    
    Para $j\in\{t_{ini},\dotsc,t_{fim}\}\land\{0,\dotsc,i-1\}$
    
    A.
    
    $s=s+a_{i,j}x_{j}$
3. (c)
  
  $\displaystyle x_{i}=\frac{b_{i}-s}{a_{i,i}}$

O código MP C/C++ que apresentaremos a seguir, faz uso do construtor threadprivate

#pragma omp threadprivate(list)

Este construtor permite que a lista de variáveis (estáticas) list seja privada para cada thread e seja compartilhada entre as regiões paralelas. Por exemplo:

x = 0
#pragma omp parallel private(x)
  x = 1
#pragma omp parallel private(x)
  x vale 0

Agora, com o construtor threadprivate:

static x = 0
#pragma omp threadprivate(x)
#pragma omp parallel
  x = 1
#pragma omp parallel private(x)
  x vale 1

Ainda, apenas para efeito de exemplo, vamos considerar que $a_{i,j}=(-1)^{i+j}(i+j)/(ij+1)$ para $i<j$ , $a_{i,i}=2[(i-n/2)^{2}+1]/n$ e $b_{i}=(-1)^{i}/(i+1)$ para $i=0,\dotsc,n-1$ .

Segue o código paralelo para a resolução direta do sistema triangular inferior. Verifique!

Código: sistria1dcol.cc

⬇

1#include <omp.h>

2#include <stdio.h>

3#include <ctime>

4#include <algorithm>

6#include <gsl/gsl_spmatrix.h>

7#include <gsl/gsl_vector.h>

8#include <gsl/gsl_rng.h>

10int np, pid;

11int ini, fim;

12#pragma omp threadprivate(np,pid,ini,fim)

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

16 int n = 9999;

18 // vetores

19 gsl_spmatrix *a = gsl_spmatrix_alloc(n,n);

20 gsl_vector *b = gsl_vector_alloc(n);

21 gsl_vector *x = gsl_vector_alloc(n);

23 // inicializacao

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

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

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

28 gsl_spmatrix_set(a, i, j,

29 pow(-1.0,i+j)*(i+j)/(i*j+1));

30 }

31 gsl_spmatrix_set(a, i, i,

32 (pow(i-n/2,2)+1)*2/n);

33 gsl_vector_set(b, i,

34 pow(-1.0,i)/(i+1));

35 }

37 printf('feito.\n');

39 printf('Executando em paralelo ... \n');

41 time_t t = time(NULL);

42 #pragma omp parallel

43 {

44 np = omp_get_num_threads();

45 pid = omp_get_thread_num();

47 ini = pid*n/np;

48 fim = (pid+1)*n/np;

49 if (pid == np-1)

50 fim = n;

51 }

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

54 double s = 0;

55 #pragma omp parallel reduction(+: s)

56 {

57 for (int j=std::max(0,ini); j<i and j<fim; j++)

58 s += gsl_spmatrix_get(a,i,j) *

59 gsl_vector_get(x,j);

60 }

61 gsl_vector_set(x, i,

62 (gsl_vector_get(b,i) - s) /

63 gsl_spmatrix_get(a,i,i));

64 }

66 t = time(NULL)-t;

68 printf('feito. %ld s\n', t);

71 gsl_spmatrix_free(a);

72 gsl_vector_free(b);

73 gsl_vector_free(x);

75 return 0;

76}

2.8.1 Exercícios resolvidos

Em remoção

ER 2.8.1.

Seja $Ax=b$ um sistema triangular inferior de dimensões $n\times n$ . O seguinte pseudocódigo paralelo é uma alternativa ao apresentado acima. Por que este pseudocódigo é mais lento que o anterior?

1.
Região paralela
1. (a)
  Para $j=0,\dotsc,n-1$
  1. i.
    
    Se $j\in\{t_{ini},\dotsc,t_{fim}\}$
    
    A.
    
    $\displaystyle x_{j}=\frac{b_{j}}{a_{j,j}}$
  2. ii.
    
    Para $i\in\{t_{ini},\dotsc,t_{fim}\}\land\{j+1,\dotsc,n-1\}$
    
    A.
    
    $b_{i}=b_{i}-a_{i,j}x_{j}$

Solução.

Este código tem um número de operações semelhante ao anterior, seu desempenho é afetado pelo chamado compartilhamento falso (false sharing). Este é um fenômeno relacionado ao uso ineficiente da memória cache de cada thread. O último laço deste pseudocódigo faz sucessivas atualizações do vetor $b$ , o que causa sucessivos recarregamentos de partes do vetor $b$ da memória RAM para a memória cache de cada um dos threads. Verifique!

ER 2.8.2.

Seja $A$ uma matriz triangular inferior e invertível de dimensões $n\times n$ . Escreva um pseudocódigo MP para calcular a matriz inversa $A^{-1}$ usando o método de substituição direta.

Solução.

Vamos denotar $A=[a_{i,j}]_{i,j=1}^{n-1}$ e $A^{-1}=[x_{i,j}]_{i,j=1}^{n-1}$ . Note que $x^{\prime}s$ são as incógnitas. Por definição, $AA^{-1}=I$ , logo

	$\displaystyle a_{1,1}x_{1,k}=\delta_{1,k}$		(2.19)
	$\displaystyle\cdots$		(2.20)
	$\displaystyle a_{i,1}x_{1,k}+\cdots+a_{i,i-1}x_{i-1,k}+a_{i,i}x_{i,k}=\delta_{% i,k}$		(2.21)
	$\displaystyle\cdots$		(2.22)
	$\displaystyle a_{n-1,1}x_{1,k}+\cdots+a_{n-1,n-1}x_{n-1,k}=\delta_{n-1,k}$		(2.23)

onde, $k=0,\dotsc,n-1$ e $\delta_{i,j}$ denota o Delta de Kronecker. Ou seja, o cálculo de $A^{-1}$ pode ser feito pela resolução de $n$ sistemas triangulares inferiores tendo $A$ como matriz de seus coeficientes.

Para construirmos um pseudocódigo MP, podemos distribuir os sistemas lineares a entre os threads disponíveis. Então, cada thread resolve em serial seus sistemas. Segue o pseudocódigo, sendo $x_{k}=(x_{1,k},\dotsc,x_{n-1,k})$ e $b_{k}=(\delta_{1,k},\dotsc,\delta_{n-1,k})$ .

1.
Região paralela
1. (a)
  Para $k\in\{t_{ini},\dotsc,t_{fim}\}$
  1. i.
    
    resolve $Ax_{k}=b_{k}$

2.8.2 Exercícios

Em remoção

E. 2.8.1.

Implemente um código MP do pseudocódigo discutido no ER 2.8.1. Compare o tempo computacional com o do código sistria1dcol.cc.

E. 2.8.2.

Implemente um código MP para computar a inversa de uma matriz triangular inferior de dimensões $n\times n$ .

E. 2.8.3.

Implemente um código MP para computar a solução de um sistema linear triangular superior de dimensões $n\times n$ .

E. 2.8.4.

Implemente um código MP para computar a inversa de uma matriz triangular superior de dimensões $n\times n$ .

Envie seu comentário

As informações preenchidas são enviadas por e-mail para o desenvolvedor do site e tratadas de forma privada. Consulte a Política de Uso de Dados para mais informações. Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

| | | |