Como discutido previamente, a iteração de Gauss-Seidel é estritamente sequencial, mas, como uma ordem de atualização diferente, podemos atualizar variáveis simultaneamente e atribuir grupos de variáveis diferentes a threads, e posteriormente, processos, diferentes.
O primeiro para isto é a implementação da atualização, sequencial ainda, das variáveis na ordem decorrente da pintura do grafo com duas cores.
Para facilitar a implementação, vamos criar uma nova classe que se comporta como um range, só que retorna primeiros os números que correspondem à cor vermelha, depois à preta, ignorando as extremidades do intervalo, onde são aplicadas as condições de contorno, lembrando que só aplicamos condições de Dirichlet.
O arquivo de apoio foi modificado para o que é mostrado abaixo.
""" gssup03.py This just organizes a few support functions for the GS iteration. This version introduces the Red-Black order range. """ import logging import numpy as np from math import pi from itertools import cycle, chain __format = "%(asctime)s: %(message)s" logging.basicConfig(format=__format, level=logging.INFO, datefmt="%H:%M:%S") #%% def zero_rhs(x, L): """Sets up a null vector for the RHS of the equation""" return np.zeros_like(x) def sin_rhs(x, L, A): """Sets up a rhs that produces a single period of a senoidal function with amplitude A as the solution.""" f = 2*pi/L return -f*f*A*np.sin(f*x) def cub_rhs(x, L, A): """Sets up rhs to have a cubic solution with zeros on 0, L/2, L/3 and value A on L. Don't forget to set up the correct boundary conditions.""" return -9*(A*L - 3*A*x)/(L*L*L) #%% class Problem: def __init__(self, L, nnode, U0=None, UL=None, calc_rhs=zero_rhs, **kwargs): """ L: float -- Length of the domain nnode: int -- Number of nodes U0, UL: float -- Dirichlet boundary conditions calc_rhs function(x, L, *args) -- computes the rhs of the equation """ self.x = np.linspace(0, L, nnode) self.dx = self.x[1] - self.x[0] self.dx2 = self.dx*self.dx self.u = np.zeros_like(self.x) self.L = L self.nn = nnode if U0 is not None: self.u[0] = U0 if UL is not None: self.u[-1] = UL self.f = calc_rhs(self.x, self.L, **kwargs) self.hist = None #%% class SymmRange: """ A range that goes forward and backward. This really is a bit of overengineering, and it's taylored to the symmetric GS iteration. Also, never tested with step different from 1, so beware. """ def __init__(self, start=0, stop=0, step=1): self.ranges = cycle([range(start, stop, step), range(stop-1, start-1, -step)]) self.current_range = iter(next(self.ranges)) def __iter__(self): return self def __next__(self): try: return next(self.current_range) except StopIteration: self.current_range = iter(next(self.ranges)) raise StopIteration #%% class RBRange: """ Returns a range **with step 1** in the red-black ordering, skipping over 0 and n-1, pretty much just useful for GS iterations. """ def __init__(self, size): self._s = size self.reset() def __iter__(self): return self def reset(self): r = range(self._s) self.c = chain(r[1:-1:2], r[2:-1:2]) def __next__(self): try: return next(self.c) except StopIteration: self.reset() raise StopIteration def info(*args, **kwargs): """Directly logs info calls, 99% of what we're doing.""" logging.info(*args, **kwargs)
Abaixo está a implementação da atualização Red-Black, comparada com a implementação simétrica. Não generalizem conclusões baseados neste único exemplo!
Este programa é essencialmente igual ao anterior, a modificação na ordem de atualização é feita simplesmente trocando o objeto que fornece a ordem de atualização.
""" Comparison of convergence rates between forward, backward and symmetric GS iterations """ import numpy as np import gssup03 as gs import matplotlib.pyplot as plt gs.info("Setting up problems.") A = 4.0 L = 3.0 nn = 21 maxiter = 1000 tol = 5e-4 probs_sym = { 'linear': gs.Problem(L, nn, UL=A, calc_rhs=gs.zero_rhs), 'sinusoidal': gs.Problem(L, nn, calc_rhs=gs.sin_rhs, A=A), 'cubic': gs.Problem(L, nn, UL=A, calc_rhs=gs.cub_rhs, A=A)} probs_rb = { 'linear': gs.Problem(L, nn, UL=A, calc_rhs=gs.zero_rhs), 'sinusoidal': gs.Problem(L, nn, calc_rhs=gs.sin_rhs, A=A), 'cubic': gs.Problem(L, nn, UL=A, calc_rhs=gs.cub_rhs, A=A)} ## Run GS iterarions symrange = gs.SymmRange(1, nn-1) rbrange = gs.RBRange(nn) orders = {'symmetric': (probs_sym, symrange), 'rb': (probs_rb, rbrange)} for order, (prob, prob_range) in orders.items(): gs.info("Running {} ordering.".format(order)) for name, p in prob.items(): gs.info("Solving {} problem.".format(name)) err = np.zeros(maxiter, dtype=np.float64) converged = False for it in range(maxiter): uold = p.u.copy() for i in prob_range: p.u[i] = 0.5*(p.u[i-1] + p.u[i+1] - p.dx2*p.f[i]) err[it] = np.linalg.norm(uold - p.u)/np.linalg.norm(p.u) if err[it] <= tol: converged = True gs.info('Problem {} converged after {} iterations'.\ format(name, it)) break if not converged: gs.info('Problem {} failed to converge after {} iterations'.\ format(name, it)) err = np.resize(err, it+1) p.hist = err figs = { 'linear': (0,3), 'sinusoidal': (1,4), 'cubic': (2,5)} styles = {'symmetric': 'r+', 'rb': 'b'} # Format graphs for name in probs_sym.keys(): plt.figure(figs[name][0]) plt.title("Solution for problem {}".format(name)) plt.figure(figs[name][1]) plt.yscale('log') plt.title("{} problem convergence history".format(name)) # Plot results for order, (prob, _) in orders.items(): for name, p in prob.items(): plt.figure(figs[name][0]) plt.plot(p.x, p.u, styles[order]) plt.figure(figs[name][1]) plt.plot(p.hist, styles[order])
O principal problema com esta implementação, comum à todas as anteriores, é que o acesso aos elementos individuais
de um vetor, por exemplo u[i], é extremamente lento! Vamos fazer um pequeno teste de desempenho para ver
qual o custo de usar o Python como uma linguagem compilada convencional.
Em primeiro lugar, vamos reimplementar a ordenação Red-Black usando operações vetoriais do numpy.