1. No category

Transformadas de Fourier
Clase 11, 13/05/2015
http://fisica.cab.cnea.gov.ar/gpgpu/index.php/en/icnpg/clases
Pero antes... streams with Thrust
Streams con Thrust
Ejemplo de clase de streams reciclado ...
● cp -r /share/apps/codigos/alumnos_icnpg2015/streams_thrust .
Incluimos thrust 1.8 headers (no necesario con cuda 7 toolkit)
● nvcc async_thrust.cu -I /share/apps/codigos/thrust-master
●
●
Ejecutamos nvprof -o a.prof ./a.out
qsub submit_gpu.sh
Verificar concurrencia de streams (overlaps copy-kernel, kernelkernel...)
Nvpp → import a.prof
Transformadas de Fourier
Clase 11, 13/05/2015
http://fisica.cab.cnea.gov.ar/gpgpu/index.php/en/icnpg/clases
Librerias hasta ahora
C++ template library en
GPU/CPU
Random Numbers en GPU/CPU
Random123:
a Library of Counter-Based
Random Number Generators
CUFFT Transformada
Rápida de Fourier en GPU
Algebra lineal densa
en GPU/CPU
¿ Que es la CUFFT ?
The FFT is a divide-and-conquer algorithm for efficiently computing discrete Fourier
transforms of complex or real-valued data sets. It is one of the most important and
widely used numerical algorithms in computational physics and general signal
processing. The CUFFT library provides a simple interface for computing FFTs on an
NVIDIA GPU, which allows users to quickly leverage the floating-point power and
parallelism of the GPU in a highly optimized and tested FFT library.
¿ Que se puede hacer con CUFFT ?
Repaso de Transformada de Fourier
“Descomposición en senos y cosenos”
h(t)
t
Simetrías
Ver Por ej:
Press etal, Numerical Recipes
Repaso de Transformada de Fourier
h(t)
t
Repaso de Transformada de Fourier
“Descomposición en senos y cosenos”
h(t)
t
Porque es tan importante la transformada de Fourier en física computacional ?
Convolution theorem
Operador...
Interacción...
Correlation theorem
Patrones magnéticos (en 2D)
Lineal pero No-local
Discretización + Transformada de Fourier
Local y lineal
Requiere transformar/antitransformar matrices en cada paso de tiempo....
Naive:
O(N*N)
PseudoEspectral:
O(N log(N))
Repaso de Transformada de Fourier
“Descomposición en senos y cosenos”
h(t)
t
Porque es tan importante la transformada de Fourier en física computacional ?
Convolution theorem
Operador...
Interacción...
Correlation theorem
Transformada de Fourier discreta
h(t)
Sampling theorem
Transformada
discreta
Transformada de Fourier discreta
h(t)
Transformada
discreta
AntiTransformada
Transformada de Fourier discreta
Transformada discreta 1D
Transformada discreta 2D
Transformada discreta 3D, etc
Transformada de Fourier discreta
h(t)
Transformada
discreta
AntiTransformada
Transformada de Fourier discreta
h(t)
Transformada
discreta
AntiTransformada
Como calcularlas eficientemente?
Algoritmos
Secuenciales
Naïve approach [hasta ~1960]: O(N*N)
Algoritmos de Fast Fourier transform (FFT)
[Cooley–Tukey o Gauss?] : O(N*log2 N) → recursivo
Fast Fourier
Transforms
●
La librería pública de FFT mas popular:
“La mas rápida del oeste”
Our benchmarks, performed on on a variety of platforms, show that FFTW's performance is typically superior to that of
other publicly available FFT software, and is even competitive with vendor-tuned codes. In contrast to vendor-tuned
codes, however, FFTW's performance is portable: the same program will perform well on most architectures without
modification. Hence the name, "FFTW," which stands for the somewhat whimsical title of "Fastest Fourier Transform in
the West."
●
Vendor-tuned codes:
MKL: Math Kernel Library
Intel's Math Kernel Library (MKL) is a library of optimized, math routines for science, engineering, and financial applications.
Core math functions include BLAS, LAPACK, ScaLAPACK, Sparse Solvers, Fast Fourier Transforms, and Vector Math.
The library supports Intel and compatible processors and is available for Windows, Linux and Mac OS X operating systems.
CUFFT
CUFFT
CUFFT Samples
Ejemplo de CUFFT
3D Complex-to-Complex Transforms
#define NX 64
#define NY 64
#define NZ 128
cufftHandle plan;
cufftComplex *data1, *data2;
cudaMalloc((void**)&data1, sizeof(cufftComplex)*NX*NY*NZ);
cudaMalloc((void**)&data2, sizeof(cufftComplex)*NX*NY*NZ);
/* Create a 3D FFT plan. */
cufftPlan3d(&plan, NX, NY, NZ, CUFFT_C2C);
/* Transform the first signal in place. */
cufftExecC2C(plan, data1, data1, CUFFT_FORWARD);
/* Transform the second signal using the same plan. */
cufftExecC2C(plan, data2, data2, CUFFT_FORWARD);
/* Destroy the cuFFT plan. */
cufftDestroy(plan);
cudaFree(data1); cudaFree(data2);
#include <cufft.h>
typedef cufftDoubleReal REAL;
typedef cufftDoubleComplex COMPLEX;
// includes de thrust
CUFFT y
Thrust
int main(){
// declaro/aloco arrays input, y su transformada, output
thrust::device_vector<REAL> input(N);
thrust::device_vector<COMPLEX> output(N/2+1);
// llenar input[] con alguna señal ...
// cast a punteros normales a memoria del device
REAL * d_input_raw = thrust::raw_pointer_cast(&input[0]);
COMPLEX * d_output_raw = thrust::raw_pointer_cast(&output[0]);
cufftHandle planfft; // declara un plan ...
cufftPlan1d(&planfft,N,CUFFT_R2C,1); // lo inicializa ...
cufftExecR2C(planfft, d_input_raw, d_output_raw); // y lo ejecuta ...
// copio la transformada del device al host, y la imprimo
thrust::host_vector<REAL> output_host = output;
for(i=0; i < N/2+1; i++) cout << output_host[i].x << “ ” << output_host[i].y << endl;
}
Ecuación de difusión o “del calor” en
dos dimensiones con fuentes arbitrarias
Fuentes dependientes
del espacio y del tiempo
Campo escalar bidimensional
Fluídos [Navier-Stokes]
Magnetismo
[Ginzburg-Landau]
Física Estadística
[Fokker-Planck]
Mecánica Cuántica
[Schrodinger]
Etc...
Ecuación de difusión o “del calor” en
dos dimensiones con fuentes arbitrarias
en el espacio de Fourier
Transformando Fourier
Modos desacoplados
Solución analítica!
Ecuación de difusión con fuentes
arbitrarias
Diferencias finitas espaciales
Esquema explícito
Esquema implícito [Crank-Nicholson]
Ecuación de difusión con fuentes
arbitrarias
Diferencias finitas
Stencil h
Convolucion en 2D
Se puede resolver en el espacio real,
Pero lo haremos en espacio de Fourier
0
C'
0
C'
-4C'+1
C'
0
C'
0
Ecuación de difusión con fuentes
arbitrarias
Stencil h
0
C'
0
C'
-4C'+1
C'
0
C'
0
Todo desacoplado :-)
Aplicar el stencil o hacer la convolucion 2D en el espacio real
se reduce a multiplicar por un campo escalar constante la transformada
Ecuación de difusión con fuentes
arbitrarias
Esto se calcula una sola vez al principio
Loop temporal
for(t=0;t<TRUN;t++)
{
}
Euler Step
Trivialmente
Paralelizable
Ecuación de difusión con fuentes
arbitrarias
Euler Step
Trivialmente
Paralelizable
// Paso de Euler en el espacio de Fourier
__global__ void euler_step
(cufftComplex *a, const cufftComplex *b,
const float *qqmodule, int nx, int ny)
{
// idx and idy, the location of the element in the original NxN array
int idx = blockIdx.x*blockDim.x+threadIdx.x;
Conviene trabajar
int idy = blockIdx.y*blockDim.y+threadIdx.y;
Con grillas 2D!
if ( idx < nx && idy < ny){
int index = idx*ny + idy;
const float fac=-(qqmodule[index]);
a[index].x = (1.f+C_BETA*fac*DT)*a[index].x + DT*b[index].x;
a[index].y = (1.f+C_BETA*fac*DT)*a[index].y + DT*b[index].y;
}
}
Ecuación de difusión con fuentes
arbitrarias
for(t=0;t<TRUN;t++)
{
Para dibujar (cada tanto) el campo en
espacio real → CUFFT → memcpy
…
Las fuentes son dadas instantaneamente
en el espacio real, pero las necesitamos
en Fourier → CUFFT
…
}
Kernel Paso de Euler
Hands-on!
Ecuación de difusión con fuentes
arbitrarias
●
●
●
●
cp -r /share/apps/codigos/alumnos_icnpg2015/heat .
make
qsub script.sh
Check .ppm images!
Ecuación de difusión con fuentes
arbitrarias
for(t=0;t<TRUN;t++)
{
Para dibujar (cada tanto) el campo en
espacio real → CUFFT → memcpy
…
Las fuentes son dadas instantaneamente
en el espacio real, pero las necesitamos
en Fourier → CUFFT
…
}
Kernel Paso de Euler
Hands-on:
Ecuación de difusión con fuentes
arbitrarias
Loop temporal principal
main.cu
En este problema:
Solo para hacer imágenes
Controla fuentes
en espacio real...
Un paso de tiempo en el
Espacio de Fourier
Hands-on:
Ecuación de difusión con fuentes
arbitrarias
Files: main.cu, heat.cu, ppmimages.cuh, Makefile, script.sh
Loop temporal
(alto nivel)
Kernels, CUFFT
Imágenes ppm
a partir del campo
en espacio real [host]
10 FIXED FIXME:
● Declaración de planes CUFFT 2D [forward y backward, R2C, C2R]
● Alocación de memoria para arrays en espacio real y complejo
● Ejecución de planes de transformada/antitransformada
● Destrucción de planes, liberación de memoria
● Varios Kernels con grillas 2D
Make & Run
● Producir imágenes / hacer película →
● Jugar con parámetros y fuentes...
● Donde y como se puede optimizar?
● Como agregar términos no lineales?
Hands-on:
Ecuación de difusión con fuentes
arbitrarias
Una configuración de fuentes divertida
“Placa
Metálica”
[Condiciones
de Contorno
Periódicas]
Temp(x,y,t)=???
“Fuente
fría”
“Fuente
caliente”
Hands-on:
Ecuación de difusión con fuentes
arbitrarias
Una configuración de fuentes divertida
“Placa
Metálica”
[Condiciones
de Contorno
Periódicas]
Referencias
●
CUFFT Library User Guide NVIDIA.
●
CUDA Libraries SDK Code Samples (simpleCUFFT, PDEs, etc).
FFTW home page: http://www.fftw.org/ [interface parecida
para CPU + openmp, mpi].
●
●
Numerical Recipes in C [algoritmo FFT explicado].
●
Libros de Numerical PDEs.... [como discretizar las PDEs?]