Parallelization of the ADI method exploring vector computing in GPUs

Abstract

The 2D convection-diffusion is a well-known problem in scientific simulation that often uses a direct method to solve a system of N linear equations, which requires N3 operations. This problem can be solved using a more efficient computational method, known as the alternating direction implicit (ADI). It solves a system of N linear equations in 2N times with N operations each, implemented in two steps, one to solve row by row, the other column by column. Each N operation is fully independent in each step, which opens an opportunity to an embarrassingly parallel solution. This method also explores the way matrices are stored in computer memory, either in row-major or column-major, by splitting each iteration in two. The major bottleneck of this method is solving the system of linear equations. These systems of linear equations can be described as tridiagonal matrices since the elements are always stored on the three main diagonals of the matrices. Algorithms tailored for tridiagonal matrices, can significantly improve the performance. These can be sequential (i.e. the Thomas algorithm) or parallel (i.e. the cyclic reduction CR, and the parallel cyclic reduction PCR). Current vector extensions in conventional scalar processing units, such as x86-64 and ARM devices, require the vector elements to be in contiguous memory locations to avoid performance penalties. To overcome these limitations in dot products several approaches are proposed and evaluated in this work, both in general-purpose processing units and in specific accelerators, namely NVidia GPUs. Profiling the code execution on a server based on x86-64 devices showed that the ADI method needs a combination of CPU computation power and memory transfer speed. This is best showed on a server based on the Intel manycore device, KNL, where the algorithm scales until the memory bandwidth is no longer enough to feed all 64 computing cores. A dual-socket server based on 16-core Xeon Skylakes, with AVX-512 vector support, proved to be a better choice: the algorithm executes in less time and scales better. The introduction of GPU computing to further improve the execution performance (and also using other optimisation techniques, namely a different thread scheme and shared memory to speed up the process) showed better results for larger grid sizes (above 32Ki x 32Ki). The CUDA development environment also showed a better performance than using OpenCL, in most cases. The largest difference was using a hybrid CR-PCR, where the OpenCL code displayed a major performance improvement when compared to CUDA. But even with this speedup, the better average time for the ADI method on all tested configurations on a NVidia GPU was using CUDA on an available updated GPU (with a Pascal architecture) and the CR as the auxiliary method.

O problema da convecção-difusão é utilizado em simulaçãos cientificas que regularmente utilizam métodos diretos para solucionar um sistema de N equações lineares e necessitam de N3 operações. O problema pode ser resolvido utilizando um método computacionalmente mais eficiente para resolver um sistema de N equações lineares com N operações cada, implementado em dois passos, um solucionando linha a linha e outro solucionando coluna a coluna. Cada par de N operações são independentes em cada passo, havendo assim uma oportunidade de utilizar uma solução em baraçosamente paralela. Este método também explora o modo de guardar as matrizes na memória do computados, sendo esta por linhas ou em colunas, dividindo cada iteração em duas, este método é conhecido como o método de direção alternada. O maior bottleneck deste problema é a resolução dos sistemas de equações lineares criados pelo ADI. Estes sistemas podem ser descritos como matrizes tridiagonais, visto todos os seus elementos se encontrarem nas 3 diagonais interiores e a utilização de métodos estudados para este caso é necessário para conseguir atingir a melhor performance possível. Esses métodos podem ser sequenciais (como o algoritmo de Thomas) ou paralelos (como o CR e o PCR) As extensões vectoriais utilizadas nas atuais unidades de processamento, como dispositivos x86-64 e ARM, necessitam que os elementos do vetor estejam em blocos de memória contíguos para não sofrer penalizações. Algumas abordagens foram estudadas neste trabalho para as ultrapassar, tanto em processadores convencionais como em aceleradores de computação. Os registos do tempo em servidores baseado em dispositivos x86-64 mostram que o ADI necessitam de uma combinação de poder de processamento assim como velocidade de transferência de dados. Isto é demonstrado especialmente no servidor baseado no dispositivo KNL da Intel, no qual o algoritmo escala até que a largura de banda deixe de ser suficiente para o problema. Um servidor com dois sockets em que cada é composto por um dispositivo com 16 cores baseado na arquitetura Xeon Skylake, com acesso ao AVX-512, mostrou ser a melhor escolha: o algoritmo faz as mesmas operações em menos tempo e escala melhor. Com a introdução de computação com GPUs para melhorar a performance do programa mostrou melhores resultados para problemas de maiores dimensões (tamanho acima de 32Ki x 32Ki celulas). O desenvolvimento em CUDA também mostrou melhores resultados que em OpenCL na maioria dos casos. A maior divergência foi observada ao utilizar o método CR-PCR, onde o OpenCL mostrou melhor performance que em CUDA. Mas mesmo com este método sendo mais eficaz que o mesmo em CUDA, o melhor performance com o método ADI foi observado utilizando CUDA no GPU mais recente estudado com o método CR.