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.9 Decomposição LU 2.11 Métodos iterativos para problemas não lineares

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

2.10 Métodos iterativos para Sistemas Lineares

Em remoção

Nesta seção, vamos discutir sobre a paralelização MP de alguns métodos iterativos para sistemas lineares

Ax=b

(2.31)

com $A=[a_{i,j}]_{i,j=0}^{n-1}$ , $x=(x_{i})_{i=0}^{n-1}$ e $b=(b_{i})_{i=0}^{n-1}$ .

2.10.1 Método de Jacobi

Em remoção

Nós podemos escrever a $i$ -ésima equação do sistema $Ax=b$ como

\sum_{j=1}^{n}a_{i,j}x_{j}=b_{i}.

(2.32)

Isolando $x_{i}$ , obtemos

x_{i}=-\frac{1}{a_{i,i}}\left[\sum_{j\neq i}a_{i,j}x_{j}-b_{i}\right].

(2.33)

Nesta última equação, temos que $x_{i}$ pode ser diretamente calculado se todos os elementos $x_{j}$ , $j\neq i$ , forem conhecidos. Isso motiva o chamado método de Jacobi que é dado pela seguinte iteração

	$\displaystyle x_{i}(0)=\text{aprox. inicial},$		(2.34)
	$\displaystyle x_{i}(t+1)=-\frac{1}{a_{i,i}}\left[\sum_{j\neq i}a_{i,j}x_{j}(t)% -b_{i}\right],$		(2.35)

para cada $i=0,1,\dotsc,n-1$ e $t=0,1,2,\ldots$ . O número máximo de iterações $t_{\text{max}}$ e o critério de parada podem ser escolhidos de forma adequada.

O pseudocódigo serial para o método de Jacobi pode ser escrito como segue

1.

Alocar a aproximação inicial $x^{0}$ .
2.
Para $t=0,1,2,\cdots,t_{\text{max}}$ :
1. (a)
  Para $i=0,1,2,\cdots,n$ :
  1. i.
    
    $\displaystyle x_{i}=-\frac{1}{a_{i,i}}\left[\sum_{j\neq i}a_{i,j}x^{0}_{j}-b_{% i}\right]$ .
2. (b)
  
  Verificar o critério de parada.
3. (c)
  
  $x0=x$ .

A paralelização MP no método de Jacobi pode ser feita de forma direta e eficaz pela distribuição do laço 2.(a) do pseudocódigo acima. O seguinte código é uma implementação MP do método de Jacobi. Os vetores $b$ e $x0$ são inicializados com elementos randômicos $(0,1)$ . A matriz $A$ é inicializada como uma matriz estritamente diagonal dominante¹¹endnote: ¹O método de Jacobi é convergente para matriz estritamente diagonal dominante. com elementos randômicos $(-1,1)$ . O critério de parada é

\|x-x^{0}\|_{2}<\text{tol},

(2.36)

onde tol é a tolerância.

Código: pJacobi.cc

⬇

1#include <omp.h>

2#include <stdio.h>

3#include <time.h>

5#include <gsl/gsl_matrix.h>

6#include <gsl/gsl_vector.h>

7#include <gsl/gsl_blas.h>

8#include <gsl/gsl_rng.h>

10// random +/- 1

11double randsig(gsl_rng *rng);

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

15 int n = 999;

16 int tmax = 50;

17 double tol = 1e-8;

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

20 gsl_vector *b = gsl_vector_alloc(n);

22 gsl_vector *x = gsl_vector_alloc(n);

23 gsl_vector *x0 = gsl_vector_alloc(n);

25 // gerador randomico

26 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

27 gsl_rng_set(rng, time(NULL));

29 // Inicializacao

30 // Matriz estritamente diagonal dominante

31 printf('Inicializacao ... \n');

32 double sig;

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

34 double s = 0;

35 for (int j=0; j<n; j++) {

36 double aux = gsl_rng_uniform(rng);

37 gsl_matrix_set(a, i, j,

38 randsig(rng)*aux);

39 s += aux;

40 }

41 gsl_matrix_set(a, i, i,

42 randsig(rng) * s);

43 gsl_vector_set(b, i,

44 randsig(rng) *

45 gsl_rng_uniform(rng));

46 gsl_vector_set(x0, i,

47 randsig(rng) *

48 gsl_rng_uniform(rng));

49 }

50 printf('feito.\n');

52 // Jacobi

53 for (int t=0; t<tmax; t++) {

54 #pragma omp parallel for

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

56 double s = 0;

57 for (int j=0; j<i; j++)

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

59 gsl_vector_get(x0, j);

60 for (int j=i+1; j<n; j++)

61 s += gsl_matrix_get(a, i, j) *

62 gsl_vector_get(x0, j);

63 gsl_vector_set(x, i,

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

65 gsl_matrix_get(a, i, i));

66 }

67 // criterio de parada

68 // ||x-x0||_2 < tol

69 gsl_blas_daxpy(-1.0, x, x0);

70 double e = gsl_blas_dnrm2(x0);

71 printf('Iter. %d: %1.0e\n', t, e);

72 if (e < tol)

73 break;

74 gsl_vector_memcpy(x0, x);

75 }

77 gsl_matrix_free(a);

78 gsl_vector_free(b);

79 gsl_vector_free(x);

80 gsl_vector_free(x0);

81 gsl_rng_free(rng);

83 return 0;

84}

86double randsig(gsl_rng *rng)

87{

88 double signal = 1.0;

89 if (gsl_rng_uniform(rng) >= 0.5)

90 signal = -1.0;

91 return signal;

92}

2.10.2 Método tipo Gauss-Seidel

Em remoção

No algoritmo serial, observamos que ao calcularmos $x_{i}$ pela iteração de Jacobi(2.33), as incógnitas $x_{j}$ , $j<i$ , já foram atualizadas. Isto motivo o método de Gauss-Seidel, cujo algoritmo é descrito no seguinte pseudocódigo:

1.

Alocar a aproximação inicial $x^{0}$ .
2.
Para $t=0,1,2,\cdots,t_{\text{max}}$ :
1. (a)
  Para $i=0,1,2,\cdots,n$ :
  1. i.
    
    $\displaystyle x_{i}=-\frac{1}{a_{i,i}}\left[\sum_{j<i}a_{i,j}x_{j}+\sum_{j>i}a% _{i,j}x^{0}_{j}-b_{i}\right]$ .
2. (b)
  
  Verificar o critério de parada.
3. (c)
  
  $x0=x$ .

Embora este método seja normalmente muito mais rápido que o método de Jacobi, ele não é paralelizável. Isto se deve ao fato de que o cálculo da incógnita $x_{i}$ depende dos cálculos precedentes das incógnitas $x_{j}$ , $j<i$ .

No entanto, a paralelização do método de Gauss-Seidel pode ser viável no caso de matrizes esparsas. Isto ocorre quando o acoplamento entre as equações não é total, podendo-se reagrupar as equações em blocos com base nos seus acoplamentos. Com isso, os blocos podem ser distribuídos entre as instâncias de processamento e, em cada uma, o método de Gauss-Seidel é aplicado de forma serial.

Uma alternativa baseada no Método de Gauss-Seidel, é utilizar o dado atualizado $x_{j}$ loco que possível, independentemente da ordem a cada iteração. A iteração do tipo Gauss-Seidel pode-se ser escrita da seguinte forma

x_{i}=-\frac{1}{a_{i,i}}\left[\sum_{\hat{j}\neq i}a_{i,\hat{j}}x_{\hat{j}}+% \sum_{j\neq i}a_{i,j}x^{0}_{j}-b_{i}\right],

(2.37)

onde arbitrariamente $\hat{j}$ correspondem aos índices para os quais $x_{\hat{j}}$ já tenham sido atualizados na iteração corrente e $j$ corresponde aos índices ainda não atualizados. O pseudocódigo MP deste método pode ser descrito como segue:

1.

Alocar a aproximação inicial $x$ .
2.
Para $t=0,1,2,\cdots,t_{\text{max}}$ :
1. (a)
  
  $x^{0}=x$ .
2. (b)
  distribuição de laço em paralelo:
  1. i.
    
    Para $i=0,1,2,\cdots,n$ :
    
    A.
    
    $\displaystyle x_{i}=-\frac{1}{a_{i,i}}\left[\sum_{j\neq i}a_{i,j}x_{j}-b_{i}\right]$ .
3. (c)
  
  Verificar o critério de parada.

Este método tipo Gauss-Seidel converge mais rápido que o método de Jacobi em muitos casos. Veja [Bertsekas2015a, p. 151–153], para alguns resultados sobre convergência.

A implementação MP do pseudocódigo acima é apresentada no código abaixo. Os elementos dos vetores $b$ , $x^{0}$ e da matriz $A$ são inicializados da mesma forma que no código pJacobi.cc acima.

Código: pGSL.cc

⬇

1#include <omp.h>

2#include <stdio.h>

3#include <time.h>

5#include <gsl/gsl_matrix.h>

6#include <gsl/gsl_vector.h>

7#include <gsl/gsl_blas.h>

8#include <gsl/gsl_rng.h>

10// random +/- 1

11double randsig(gsl_rng *rng);

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

15 int n = 999;

16 int tmax = 50;

17 double tol = 1e-8;

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

20 gsl_vector *b = gsl_vector_alloc(n);

22 gsl_vector *x = gsl_vector_alloc(n);

23 gsl_vector *x0 = gsl_vector_alloc(n);

25 // gerador randomico

26 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

27 gsl_rng_set(rng, time(NULL));

29 // Inicializacao

30 // Matriz estritamente diagonal dominante

31 printf('Inicializacao ... \n');

32 double sig;

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

34 double s = 0;

35 for (int j=0; j<n; j++) {

36 double aux = gsl_rng_uniform(rng);

37 gsl_matrix_set(a, i, j,

38 randsig(rng)*aux);

39 s += aux;

40 }

41 gsl_matrix_set(a, i, i,

42 randsig(rng) * s);

43 gsl_vector_set(b, i,

44 randsig(rng) *

45 gsl_rng_uniform(rng));

46 gsl_vector_set(x, i,

47 randsig(rng) *

48 gsl_rng_uniform(rng));

49 }

50 printf('feito.\n');

52 // Random Gauss-Seidel

53 for (int t=0; t<tmax; t++) {

54 gsl_vector_memcpy(x0, x);

55 #pragma omp parallel for

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

57 double s = 0;

58 for (int j=0; j<i; j++)

59 s += gsl_matrix_get(a, i, j) *

60 gsl_vector_get(x, j);

61 for (int j=i+1; j<n; j++)

62 s += gsl_matrix_get(a, i, j) *

63 gsl_vector_get(x, j);

64 gsl_vector_set(x, i,

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

66 gsl_matrix_get(a, i, i));

67 }

68 // criterio de parada

69 // ||x-x0||_2 < tol

70 gsl_blas_daxpy(-1.0, x, x0);

71 double e = gsl_blas_dnrm2(x0);

72 printf('Iter. %d: %1.0e\n', t, e);

73 if (e < tol)

74 break;

75 }

77 gsl_matrix_free(a);

78 gsl_vector_free(b);

79 gsl_vector_free(x);

80 gsl_vector_free(x0);

81 gsl_rng_free(rng);

83 return 0;

84}

86double randsig(gsl_rng *rng)

87{

88 double signal = 1.0;

89 if (gsl_rng_uniform(rng) >= 0.5)

90 signal = -1.0;

91 return signal;

92}

2.10.3 Método do Gradiente Conjugado

Em remoção

O Método do Gradiente Conjugado pode ser utilizado na resolução de sistemas lineares $Ax=b$ , onde $A$ é uma matriz simétrica e positiva definida. No caso de sistemas em gerais, o método pode ser utilizado para resolver o sistema equivalente $A^{\prime}Ax=A^{\prime}b$ , onde $A$ é uma matriz inversível, com $A^{\prime}$ denotando a transposta de $A$ .

O pseudocódigo deste método é apresentado como segue:

1.

Alocar a aproximação inicial $x$ .
2.

Calcular o resíduo $r=Ax-b$ .
3.

Alocar a direção $d=r$ .
4.
Para $t=0,1,\cdots,t_{\text{max}}$ :
1. (a)
  
  $\displaystyle\alpha=-\frac{r\cdot d}{d\cdot Ad}$ .
2. (b)
  
  $x=x+\alpha d$ .
3. (c)
  
  $r=Ax-b$ .
4. (d)
  
  $\displaystyle\beta=\frac{r\cdot Ad}{d\cdot Ad}$ .
5. (e)
  
  $d=-r+\beta d$

Uma versão MP deste método pode ser implementada pela distribuição em paralelo de cada uma das operações de produto escalar, multiplicação matriz-vetor e soma vetor-vetor. O seguinte código é uma implementação MP do Método do Gradiente Conjugado. Os elementos do vetor $b$ e da matriz $A$ são inicializados de forma randômica e é garantida que matriz é simétrica positiva definida.

Código: pGC.cc

⬇

1#include <omp.h>

2#include <stdio.h>

3#include <time.h>

4#include <math.h>

6#include <gsl/gsl_matrix.h>

7#include <gsl/gsl_vector.h>

8#include <gsl/gsl_blas.h>

9#include <gsl/gsl_rng.h>

11int n = 999;

12int tmax = 50;

13double tol = 1e-8;

15// inicializacao

16void init(gsl_matrix *a,

17 gsl_vector *b);

19// random +/- 1

20double randsig(gsl_rng *rng);

22// residuo

23void residuo(const gsl_matrix *a,

24 const gsl_vector *x,

25 const gsl_vector *b,

26 gsl_vector *r);

28// Metodo do Gradiente Conjugado

29void pGC(const gsl_matrix *a,

30 const gsl_vector *b,

31 gsl_vector *x);

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

35 // sistema

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

37 gsl_vector *b = gsl_vector_alloc(n);

39 // incognita

40 gsl_vector *x = gsl_vector_alloc(n);

42 // inicializacao

43 init(a, b);

45 // Metodo do Gradiente Conjugado

46 pGC(a, b, x);

48 gsl_matrix_free(a);

49 gsl_vector_free(b);

50 gsl_vector_free(x);

52 return 0;

53}

55/******************************************

56Inicializacao

57******************************************/

58void init(gsl_matrix *a,

59 gsl_vector *b)

60{

61 printf('Inicializacao ... \n');

62 // gerador randomico

63 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

64 gsl_rng_set(rng, time(NULL));

66 // C - Matriz estritamente diagonal positiva

67 double sig;

68 gsl_matrix *c = gsl_matrix_alloc(n,n);

69 #pragma omp parallel for

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

71 double aux;

72 double s = 0;

73 for (int j=0; j<n; j++) {

74 aux = gsl_rng_uniform(rng);

75 gsl_matrix_set(c, i, j,

76 randsig(rng) * aux);

77 s += aux;

78 }

79 gsl_matrix_set(c, i, i,

80 randsig(rng) * s);

81 gsl_vector_set(b, i,

82 randsig(rng) *

83 gsl_rng_uniform(rng));

84 }

85 // A = C'C: Simetrica positiva definida

86 #pragma omp parallel for

87 for (int i=0; i<n; i++)

88 for (int j=0; j<n; j++) {

89 double s;

90 gsl_vector_const_view ci =

91 gsl_matrix_const_column(c, i);

92 gsl_vector_const_view cj =

93 gsl_matrix_const_column(c, j);

94 gsl_blas_ddot(&ci.vector, &cj.vector, &s);

95 gsl_matrix_set(a, i, j, s);

96 }

98 gsl_rng_free(rng);

99 gsl_matrix_free(c);

100

101 printf('feito.\n');

102}

103/*****************************************/

104

105/*****************************************

106Sinal randomico

107*****************************************/

108double randsig(gsl_rng *rng)

109{

110 double signal = 1.0;

111 if (gsl_rng_uniform(rng) >= 0.5)

112 signal = -1.0;

113 return signal;

114}

115/***************************************/

116

117/***************************************

118residuo

119***************************************/

120void residuo(const gsl_matrix *a,

121 const gsl_vector *x,

122 const gsl_vector *b,

123 gsl_vector *r)

124{

125 #pragma omp parallel for

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

127 double s = 0;

128 for (int j=0; j<n; j++)

129 s += gsl_matrix_get(a, i, j) *

130 gsl_vector_get(x, j);

131 gsl_vector_set(r, i,

132 s - gsl_vector_get(b, i));

133 }

134}

135/*************************************/

136

137/***************************************

138Metodo do Gradiente Conjugado

139***************************************/

140void pGC(const gsl_matrix *a,

141 const gsl_vector *b,

142 gsl_vector *x)

143{

144 gsl_vector *r = gsl_vector_alloc(n);

145 gsl_vector *d = gsl_vector_alloc(n);

146 gsl_vector *ad = gsl_vector_alloc(n);

147

148 // x = 0

149 gsl_vector_set_zero(x);

150

151 // r = Ax - b

152 residuo(a, x, b, r);

153

154 // d = r

155 gsl_vector_memcpy(d, r);

156

157 for (int t=0; t<tmax; t++) {

158 // r.d, Ad, dAd

159 double rd = 0;

160 double dAd = 0;

161 #pragma omp parallel for reduction(+:rd,dAd)

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

163 rd += gsl_vector_get(r, i) *

164 gsl_vector_get(d, i);

165 double adi = 0;

166 for (int j=0; j<n; j++)

167 adi += gsl_matrix_get(a, i, j) *

168 gsl_vector_get(d, j);

169 gsl_vector_set(ad, i, adi);

170 dAd += gsl_vector_get(d, i) * adi;

171 }

172

173 // alpha

174 double alpha = rd/dAd;

175

176 // x = x - alpha*d

177 #pragma omp parallel for

178 for (int i=0; i<n; i++)

179 gsl_vector_set(x, i,

180 gsl_vector_get(x, i) -

181 alpha *

182 gsl_vector_get(d, i));

183

184 // residuo

185 residuo(a, x, b, r);

186

187 // rAd

188 double rAd = 0;

189 #pragma omp parallel for reduction(+:rAd)

190 for (int i=0; i<n; i++)

191 rAd += gsl_vector_get(r, i) *

192 gsl_vector_get(ad, i);

193

194 // beta

195 double beta = rAd/dAd;

196

197 // d

198 #pragma omp parallel for

199 for (int i=0; i<n; i++)

200 gsl_vector_set(d, i,

201 beta *

202 gsl_vector_get(d, i) -

203 gsl_vector_get(r, i));

204

205 // criterio de parada

206 // ||r||_2 < tol

207 double crt = 0;

208 #pragma omp parallel for reduction(+: crt)

209 for (int i=0; i<n; i++)

210 crt += gsl_vector_get(r, i) *

211 gsl_vector_get(r, i);

212 crt = sqrt(crt);

213 printf('Iter. %d: %1.1e\n', t, crt);

214 if (crt < tol)

215 break;

216 }

217

218 gsl_vector_free(r);

219 gsl_vector_free(d);

220 gsl_vector_free(ad);

221

222}

223/*************************************/

2.10.4 Exercícios Resolvidos

Em remoção

ER 2.10.1.

Faça uma implementação MP para computar a inversa de uma matriz $A$ usando o Método de Gauss-Seidel. Assuma que $A$ seja uma matriz estritamente diagonal dominante de dimensões $n\times n$ ( $n$ grande).

Solução 0.

A inversa da matriz $A$ é a matriz $B$ de dimensões $n\times n$ tal que

AB=I

(2.38)

Denotando por $b_{k}$ , $k=0,1,\dotsc,n$ , as colunas da matriz $B$ , temos que o problema de calcular $B$ é equivalente a resolver os seguintes $n$ sistemas lineares

Ab_{k}=i_{k},\quad,k=0,1,\dotsc,n,

(2.39)

onde $i_{k}$ é a j-ésima coluna da matriz identidade $I$ . Podemos usar o método de Gauss-Seidel para computar a solução de cada um destes sistemas lineares. Embora o método não seja paralelizável, os sistemas são independentes um dos outros e podem ser computados em paralelo. O pseudocódigo pode ser escrito como segue:

1.

Alocar a matriz A.
2.
(início da região paralela)
1. (a)
  Para $k=0,1,\dotsc,n$ (laço em paralelo):
  1. i.
    
    Alocar $i_{k}$ .
  2. ii.
    
    Inicializar $b_{k}$ .
  3. iii.
    
    Resolver pelo Método de Gauss-Seidel
    
    $Ab_{k}=i_{k}$ (2.40)

A implementação fica como Exercício E 2.10.2.

ER 2.10.2.

Faça uma implementação MP do método de sobre-relaxação de Jacobi (método JOR) para computar a solução de um sistema linear $Ax=b$ , com $A$ matriz estritamente diagonal dominante de dimensões $n\times n$ ( $n$ grande).

Solução 0.

O método JOR é uma variante do método de Jacobi. A iteração JOR é

	$\displaystyle x_{i}(0)=\text{aprox. inicial},$		(2.41)
	$\displaystyle x_{i}(t+1)=(1-\gamma)x_{i}(t)-\frac{\gamma}{a_{i,i}}\left[\sum_{% j\neq i}a_{i,j}x_{j}(t)-b_{i}\right],$		(2.42)

para cada $i=0,1,\dotsc,n-1$ e $t=0,1,2,\ldots$ , com $0<\gamma<1$ . Note que se $\gamma=1$ , então temos o Método de Jacobi.

A implementação MP do Método JOR pode ser feita de forma análoga a do Método de Jacobi (veja o código pJacobi.cc na Subseção 2.10.1). A implementação fica como exercício E 2.10.1.

2.10.5 Exercícios

Em remoção

E. 2.10.1.

Complete o ER 2.10.2.

E. 2.10.2.

Complete o ER 2.10.1.

E. 2.10.3.

O Método de Richardson para o cálculo da solução de um sistema linear $Ax=b$ de dimensões $n\times n$ tem a seguinte iteração

	$\displaystyle x(0)=\text{aprox. inicial},$		(2.43)
	$\displaystyle x(t+1)=x(t)-\gamma\left[Ax(t)-b\right],$		(2.44)

onde $\gamma$ é uma parâmetro escalar de relaxação e $t=0,1,2,\ldots$ . Faça uma implementação MP deste método.

E. 2.10.4.

O Método das Sucessivas Sobre-relaxações (SOR) é uma variante do Método de Gauss-Seidel. A iteração SOR é

	$\displaystyle x_{i}(0)=\text{aprox. inicial},$		(2.45)
	$\displaystyle x_{i}(t+1)=(1-\gamma)x_{i}(t)-\frac{\gamma}{a_{i,i}}\left[\sum_{% j<i}a_{i,j}x_{j}(t+1)+\sum_{j>i}a_{i,j}x_{j}(t)-b_{i}\right],$		(2.46)

onde $0<\gamma<1$ , $i=0,1,\dotsc,n-1$ e $t=0,1,2,\ldots$ .

Este método não é paralelizável, mas ele pode ser adaptado pela distribuição paralela do cálculo das incógnitas a cada iteração conforme o Método tipo Gauss-Seidel apresentado na Subseção 2.10.2. Faça a adaptação do Método SOR e implemente em MP.

E. 2.10.5.

Faça a implementação do método do Gradiente Conjugado para computar a inversa de uma matriz $A$ simétrica positiva definida de dimensões $n\times n$ ( $n$ grande).

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.9 Decomposição LU 2.11 Métodos iterativos para problemas não lineares

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

2.10 Métodos iterativos para Sistemas Lineares

Em remoção

Nesta seção, vamos discutir sobre a paralelização MP de alguns métodos iterativos para sistemas lineares

Ax=b

(2.31)

com $A=[a_{i,j}]_{i,j=0}^{n-1}$ , $x=(x_{i})_{i=0}^{n-1}$ e $b=(b_{i})_{i=0}^{n-1}$ .

2.10.1 Método de Jacobi

Em remoção

Nós podemos escrever a $i$ -ésima equação do sistema $Ax=b$ como

\sum_{j=1}^{n}a_{i,j}x_{j}=b_{i}.

(2.32)

Isolando $x_{i}$ , obtemos

x_{i}=-\frac{1}{a_{i,i}}\left[\sum_{j\neq i}a_{i,j}x_{j}-b_{i}\right].

(2.33)

	$\displaystyle x_{i}(0)=\text{aprox. inicial},$		(2.34)
	$\displaystyle x_{i}(t+1)=-\frac{1}{a_{i,i}}\left[\sum_{j\neq i}a_{i,j}x_{j}(t)% -b_{i}\right],$		(2.35)

para cada $i=0,1,\dotsc,n-1$ e $t=0,1,2,\ldots$ . O número máximo de iterações $t_{\text{max}}$ e o critério de parada podem ser escolhidos de forma adequada.

O pseudocódigo serial para o método de Jacobi pode ser escrito como segue

1.

Alocar a aproximação inicial $x^{0}$ .
2.
Para $t=0,1,2,\cdots,t_{\text{max}}$ :
1. (a)
  Para $i=0,1,2,\cdots,n$ :
  1. i.
    
    $\displaystyle x_{i}=-\frac{1}{a_{i,i}}\left[\sum_{j\neq i}a_{i,j}x^{0}_{j}-b_{% i}\right]$ .
2. (b)
  
  Verificar o critério de parada.
3. (c)
  
  $x0=x$ .

\|x-x^{0}\|_{2}<\text{tol},

(2.36)

onde tol é a tolerância.

Código: pJacobi.cc

⬇

1#include <omp.h>

2#include <stdio.h>

3#include <time.h>

5#include <gsl/gsl_matrix.h>

6#include <gsl/gsl_vector.h>

7#include <gsl/gsl_blas.h>

8#include <gsl/gsl_rng.h>

10// random +/- 1

11double randsig(gsl_rng *rng);

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

15 int n = 999;

16 int tmax = 50;

17 double tol = 1e-8;

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

20 gsl_vector *b = gsl_vector_alloc(n);

22 gsl_vector *x = gsl_vector_alloc(n);

23 gsl_vector *x0 = gsl_vector_alloc(n);

25 // gerador randomico

26 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

27 gsl_rng_set(rng, time(NULL));

29 // Inicializacao

30 // Matriz estritamente diagonal dominante

31 printf('Inicializacao ... \n');

32 double sig;

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

34 double s = 0;

35 for (int j=0; j<n; j++) {

36 double aux = gsl_rng_uniform(rng);

37 gsl_matrix_set(a, i, j,

38 randsig(rng)*aux);

39 s += aux;

40 }

41 gsl_matrix_set(a, i, i,

42 randsig(rng) * s);

43 gsl_vector_set(b, i,

44 randsig(rng) *

45 gsl_rng_uniform(rng));

46 gsl_vector_set(x0, i,

47 randsig(rng) *

48 gsl_rng_uniform(rng));

49 }

50 printf('feito.\n');

52 // Jacobi

53 for (int t=0; t<tmax; t++) {

54 #pragma omp parallel for

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

56 double s = 0;

57 for (int j=0; j<i; j++)

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

59 gsl_vector_get(x0, j);

60 for (int j=i+1; j<n; j++)

61 s += gsl_matrix_get(a, i, j) *

62 gsl_vector_get(x0, j);

63 gsl_vector_set(x, i,

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

65 gsl_matrix_get(a, i, i));

66 }

67 // criterio de parada

68 // ||x-x0||_2 < tol

69 gsl_blas_daxpy(-1.0, x, x0);

70 double e = gsl_blas_dnrm2(x0);

71 printf('Iter. %d: %1.0e\n', t, e);

72 if (e < tol)

73 break;

74 gsl_vector_memcpy(x0, x);

75 }

77 gsl_matrix_free(a);

78 gsl_vector_free(b);

79 gsl_vector_free(x);

80 gsl_vector_free(x0);

81 gsl_rng_free(rng);

83 return 0;

84}

86double randsig(gsl_rng *rng)

87{

88 double signal = 1.0;

89 if (gsl_rng_uniform(rng) >= 0.5)

90 signal = -1.0;

91 return signal;

92}

2.10.2 Método tipo Gauss-Seidel

Em remoção

1.

Alocar a aproximação inicial $x^{0}$ .
2.
Para $t=0,1,2,\cdots,t_{\text{max}}$ :
1. (a)
  Para $i=0,1,2,\cdots,n$ :
  1. i.
    
    $\displaystyle x_{i}=-\frac{1}{a_{i,i}}\left[\sum_{j<i}a_{i,j}x_{j}+\sum_{j>i}a% _{i,j}x^{0}_{j}-b_{i}\right]$ .
2. (b)
  
  Verificar o critério de parada.
3. (c)
  
  $x0=x$ .

x_{i}=-\frac{1}{a_{i,i}}\left[\sum_{\hat{j}\neq i}a_{i,\hat{j}}x_{\hat{j}}+% \sum_{j\neq i}a_{i,j}x^{0}_{j}-b_{i}\right],

(2.37)

1.

Alocar a aproximação inicial $x$ .
2.
Para $t=0,1,2,\cdots,t_{\text{max}}$ :
1. (a)
  
  $x^{0}=x$ .
2. (b)
  distribuição de laço em paralelo:
  1. i.
    
    Para $i=0,1,2,\cdots,n$ :
    
    A.
    
    $\displaystyle x_{i}=-\frac{1}{a_{i,i}}\left[\sum_{j\neq i}a_{i,j}x_{j}-b_{i}\right]$ .
3. (c)
  
  Verificar o critério de parada.

Este método tipo Gauss-Seidel converge mais rápido que o método de Jacobi em muitos casos. Veja [Bertsekas2015a, p. 151–153], para alguns resultados sobre convergência.

Código: pGSL.cc

⬇

1#include <omp.h>

2#include <stdio.h>

3#include <time.h>

5#include <gsl/gsl_matrix.h>

6#include <gsl/gsl_vector.h>

7#include <gsl/gsl_blas.h>

8#include <gsl/gsl_rng.h>

10// random +/- 1

11double randsig(gsl_rng *rng);

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

15 int n = 999;

16 int tmax = 50;

17 double tol = 1e-8;

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

20 gsl_vector *b = gsl_vector_alloc(n);

22 gsl_vector *x = gsl_vector_alloc(n);

23 gsl_vector *x0 = gsl_vector_alloc(n);

25 // gerador randomico

26 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

27 gsl_rng_set(rng, time(NULL));

29 // Inicializacao

30 // Matriz estritamente diagonal dominante

31 printf('Inicializacao ... \n');

32 double sig;

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

34 double s = 0;

35 for (int j=0; j<n; j++) {

36 double aux = gsl_rng_uniform(rng);

37 gsl_matrix_set(a, i, j,

38 randsig(rng)*aux);

39 s += aux;

40 }

41 gsl_matrix_set(a, i, i,

42 randsig(rng) * s);

43 gsl_vector_set(b, i,

44 randsig(rng) *

45 gsl_rng_uniform(rng));

46 gsl_vector_set(x, i,

47 randsig(rng) *

48 gsl_rng_uniform(rng));

49 }

50 printf('feito.\n');

52 // Random Gauss-Seidel

53 for (int t=0; t<tmax; t++) {

54 gsl_vector_memcpy(x0, x);

55 #pragma omp parallel for

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

57 double s = 0;

58 for (int j=0; j<i; j++)

59 s += gsl_matrix_get(a, i, j) *

60 gsl_vector_get(x, j);

61 for (int j=i+1; j<n; j++)

62 s += gsl_matrix_get(a, i, j) *

63 gsl_vector_get(x, j);

64 gsl_vector_set(x, i,

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

66 gsl_matrix_get(a, i, i));

67 }

68 // criterio de parada

69 // ||x-x0||_2 < tol

70 gsl_blas_daxpy(-1.0, x, x0);

71 double e = gsl_blas_dnrm2(x0);

72 printf('Iter. %d: %1.0e\n', t, e);

73 if (e < tol)

74 break;

75 }

77 gsl_matrix_free(a);

78 gsl_vector_free(b);

79 gsl_vector_free(x);

80 gsl_vector_free(x0);

81 gsl_rng_free(rng);

83 return 0;

84}

86double randsig(gsl_rng *rng)

87{

88 double signal = 1.0;

89 if (gsl_rng_uniform(rng) >= 0.5)

90 signal = -1.0;

91 return signal;

92}

2.10.3 Método do Gradiente Conjugado

Em remoção

O pseudocódigo deste método é apresentado como segue:

1.

Alocar a aproximação inicial $x$ .
2.

Calcular o resíduo $r=Ax-b$ .
3.

Alocar a direção $d=r$ .
4.
Para $t=0,1,\cdots,t_{\text{max}}$ :
1. (a)
  
  $\displaystyle\alpha=-\frac{r\cdot d}{d\cdot Ad}$ .
2. (b)
  
  $x=x+\alpha d$ .
3. (c)
  
  $r=Ax-b$ .
4. (d)
  
  $\displaystyle\beta=\frac{r\cdot Ad}{d\cdot Ad}$ .
5. (e)
  
  $d=-r+\beta d$

Código: pGC.cc

⬇

1#include <omp.h>

2#include <stdio.h>

3#include <time.h>

4#include <math.h>

6#include <gsl/gsl_matrix.h>

7#include <gsl/gsl_vector.h>

8#include <gsl/gsl_blas.h>

9#include <gsl/gsl_rng.h>

11int n = 999;

12int tmax = 50;

13double tol = 1e-8;

15// inicializacao

16void init(gsl_matrix *a,

17 gsl_vector *b);

19// random +/- 1

20double randsig(gsl_rng *rng);

22// residuo

23void residuo(const gsl_matrix *a,

24 const gsl_vector *x,

25 const gsl_vector *b,

26 gsl_vector *r);

28// Metodo do Gradiente Conjugado

29void pGC(const gsl_matrix *a,

30 const gsl_vector *b,

31 gsl_vector *x);

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

35 // sistema

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

37 gsl_vector *b = gsl_vector_alloc(n);

39 // incognita

40 gsl_vector *x = gsl_vector_alloc(n);

42 // inicializacao

43 init(a, b);

45 // Metodo do Gradiente Conjugado

46 pGC(a, b, x);

48 gsl_matrix_free(a);

49 gsl_vector_free(b);

50 gsl_vector_free(x);

52 return 0;

53}

55/******************************************

56Inicializacao

57******************************************/

58void init(gsl_matrix *a,

59 gsl_vector *b)

60{

61 printf('Inicializacao ... \n');

62 // gerador randomico

63 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

64 gsl_rng_set(rng, time(NULL));

66 // C - Matriz estritamente diagonal positiva

67 double sig;

68 gsl_matrix *c = gsl_matrix_alloc(n,n);

69 #pragma omp parallel for

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

71 double aux;

72 double s = 0;

73 for (int j=0; j<n; j++) {

74 aux = gsl_rng_uniform(rng);

75 gsl_matrix_set(c, i, j,

76 randsig(rng) * aux);

77 s += aux;

78 }

79 gsl_matrix_set(c, i, i,

80 randsig(rng) * s);

81 gsl_vector_set(b, i,

82 randsig(rng) *

83 gsl_rng_uniform(rng));

84 }

85 // A = C'C: Simetrica positiva definida

86 #pragma omp parallel for

87 for (int i=0; i<n; i++)

88 for (int j=0; j<n; j++) {

89 double s;

90 gsl_vector_const_view ci =

91 gsl_matrix_const_column(c, i);

92 gsl_vector_const_view cj =

93 gsl_matrix_const_column(c, j);

94 gsl_blas_ddot(&ci.vector, &cj.vector, &s);

95 gsl_matrix_set(a, i, j, s);

96 }

98 gsl_rng_free(rng);

99 gsl_matrix_free(c);

100

101 printf('feito.\n');

102}

103/*****************************************/

104

105/*****************************************

106Sinal randomico

107*****************************************/

108double randsig(gsl_rng *rng)

109{

110 double signal = 1.0;

111 if (gsl_rng_uniform(rng) >= 0.5)

112 signal = -1.0;

113 return signal;

114}

115/***************************************/

116

117/***************************************

118residuo

119***************************************/

120void residuo(const gsl_matrix *a,

121 const gsl_vector *x,

122 const gsl_vector *b,

123 gsl_vector *r)

124{

125 #pragma omp parallel for

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

127 double s = 0;

128 for (int j=0; j<n; j++)

129 s += gsl_matrix_get(a, i, j) *

130 gsl_vector_get(x, j);

131 gsl_vector_set(r, i,

132 s - gsl_vector_get(b, i));

133 }

134}

135/*************************************/

136

137/***************************************

138Metodo do Gradiente Conjugado

139***************************************/

140void pGC(const gsl_matrix *a,

141 const gsl_vector *b,

142 gsl_vector *x)

143{

144 gsl_vector *r = gsl_vector_alloc(n);

145 gsl_vector *d = gsl_vector_alloc(n);

146 gsl_vector *ad = gsl_vector_alloc(n);

147

148 // x = 0

149 gsl_vector_set_zero(x);

150

151 // r = Ax - b

152 residuo(a, x, b, r);

153

154 // d = r

155 gsl_vector_memcpy(d, r);

156

157 for (int t=0; t<tmax; t++) {

158 // r.d, Ad, dAd

159 double rd = 0;

160 double dAd = 0;

161 #pragma omp parallel for reduction(+:rd,dAd)

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

163 rd += gsl_vector_get(r, i) *

164 gsl_vector_get(d, i);

165 double adi = 0;

166 for (int j=0; j<n; j++)

167 adi += gsl_matrix_get(a, i, j) *

168 gsl_vector_get(d, j);

169 gsl_vector_set(ad, i, adi);

170 dAd += gsl_vector_get(d, i) * adi;

171 }

172

173 // alpha

174 double alpha = rd/dAd;

175

176 // x = x - alpha*d

177 #pragma omp parallel for

178 for (int i=0; i<n; i++)

179 gsl_vector_set(x, i,

180 gsl_vector_get(x, i) -

181 alpha *

182 gsl_vector_get(d, i));

183

184 // residuo

185 residuo(a, x, b, r);

186

187 // rAd

188 double rAd = 0;

189 #pragma omp parallel for reduction(+:rAd)

190 for (int i=0; i<n; i++)

191 rAd += gsl_vector_get(r, i) *

192 gsl_vector_get(ad, i);

193

194 // beta

195 double beta = rAd/dAd;

196

197 // d

198 #pragma omp parallel for

199 for (int i=0; i<n; i++)

200 gsl_vector_set(d, i,

201 beta *

202 gsl_vector_get(d, i) -

203 gsl_vector_get(r, i));

204

205 // criterio de parada

206 // ||r||_2 < tol

207 double crt = 0;

208 #pragma omp parallel for reduction(+: crt)

209 for (int i=0; i<n; i++)

210 crt += gsl_vector_get(r, i) *

211 gsl_vector_get(r, i);

212 crt = sqrt(crt);

213 printf('Iter. %d: %1.1e\n', t, crt);

214 if (crt < tol)

215 break;

216 }

217

218 gsl_vector_free(r);

219 gsl_vector_free(d);

220 gsl_vector_free(ad);

221

222}

223/*************************************/

2.10.4 Exercícios Resolvidos

Em remoção

ER 2.10.1.

Solução 0.

A inversa da matriz $A$ é a matriz $B$ de dimensões $n\times n$ tal que

AB=I

(2.38)

Denotando por $b_{k}$ , $k=0,1,\dotsc,n$ , as colunas da matriz $B$ , temos que o problema de calcular $B$ é equivalente a resolver os seguintes $n$ sistemas lineares

Ab_{k}=i_{k},\quad,k=0,1,\dotsc,n,

(2.39)

1.

Alocar a matriz A.
2.
(início da região paralela)
1. (a)
  Para $k=0,1,\dotsc,n$ (laço em paralelo):
  1. i.
    
    Alocar $i_{k}$ .
  2. ii.
    
    Inicializar $b_{k}$ .
  3. iii.
    
    Resolver pelo Método de Gauss-Seidel
    
    $Ab_{k}=i_{k}$ (2.40)

A implementação fica como Exercício E 2.10.2.

ER 2.10.2.

Solução 0.

O método JOR é uma variante do método de Jacobi. A iteração JOR é

	$\displaystyle x_{i}(0)=\text{aprox. inicial},$		(2.41)
	$\displaystyle x_{i}(t+1)=(1-\gamma)x_{i}(t)-\frac{\gamma}{a_{i,i}}\left[\sum_{% j\neq i}a_{i,j}x_{j}(t)-b_{i}\right],$		(2.42)

para cada $i=0,1,\dotsc,n-1$ e $t=0,1,2,\ldots$ , com $0<\gamma<1$ . Note que se $\gamma=1$ , então temos o Método de Jacobi.

A implementação MP do Método JOR pode ser feita de forma análoga a do Método de Jacobi (veja o código pJacobi.cc na Subseção 2.10.1). A implementação fica como exercício E 2.10.1.

2.10.5 Exercícios

Em remoção

E. 2.10.1.

Complete o ER 2.10.2.

E. 2.10.2.

Complete o ER 2.10.1.

E. 2.10.3.

O Método de Richardson para o cálculo da solução de um sistema linear $Ax=b$ de dimensões $n\times n$ tem a seguinte iteração

	$\displaystyle x(0)=\text{aprox. inicial},$		(2.43)
	$\displaystyle x(t+1)=x(t)-\gamma\left[Ax(t)-b\right],$		(2.44)

onde $\gamma$ é uma parâmetro escalar de relaxação e $t=0,1,2,\ldots$ . Faça uma implementação MP deste método.

E. 2.10.4.

O Método das Sucessivas Sobre-relaxações (SOR) é uma variante do Método de Gauss-Seidel. A iteração SOR é

	$\displaystyle x_{i}(0)=\text{aprox. inicial},$		(2.45)
	$\displaystyle x_{i}(t+1)=(1-\gamma)x_{i}(t)-\frac{\gamma}{a_{i,i}}\left[\sum_{% j<i}a_{i,j}x_{j}(t+1)+\sum_{j>i}a_{i,j}x_{j}(t)-b_{i}\right],$		(2.46)

onde $0<\gamma<1$ , $i=0,1,\dotsc,n-1$ e $t=0,1,2,\ldots$ .

E. 2.10.5.

Faça a implementação do método do Gradiente Conjugado para computar a inversa de uma matriz $A$ simétrica positiva definida de dimensões $n\times n$ ( $n$ grande).

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

| | | |