FFTW

      Commentaires fermés sur FFTW

Description

FFTW est une bibliothèque C open source dédiée au calcul rapide de la transformée de Fourier discrète (DFT) en une ou plusieurs dimensions.

Elle fournit des fonctions permettant de calculer des transformées de Fourier complexes ou réelles, ainsi que leurs inverses, en optimisant automatiquement les algorithmes utilisés en fonction de la taille des données et de l’architecture matérielle.

FFTW s’appuie sur un mécanisme de planification (« planning ») qui évalue différentes stratégies de calcul afin de sélectionner la plus performante pour un problème donné, offrant ainsi d’excellentes performances sur un large éventail de plateformes.

L’objectif est d’offrir une bibliothèque portable, robuste et hautement optimisée, adaptée aux applications de traitement du signal, de traitement d’images, de simulation scientifique, d’analyse spectrale et, plus généralement, à tous les domaines nécessitant des calculs efficaces de transformées de Fourier.

Mise en place de l’environnement

ml fftw
  • Version(s) disponible(s) :
    • 3.3.10
    • 3.3.11 (Construite avec GCC 15.2.0) – par défaut
    • 3.3.11-intelmpi-2021.16 (Construite avec GCC 15.2.0 et Intel MPI 2021.16)

Tutoriel

Cas d’usage : FFT 1D simple

  • Créer le fichier fftw_simple.c
#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>

int main()
{
    int N = 8;

    // Tableau d'entrée réel
    double *in;

    // Tableau de sortie complexe
    fftw_complex *out;

    // Allocation mémoire alignée FFTW
    in  = fftw_malloc(sizeof(double) * N);
    out = fftw_malloc(sizeof(fftw_complex) * N);

    // Exemple de signal
    for (int i = 0; i < N; i++)
        in[i] = i;

    // Création du plan FFT
    fftw_plan plan = fftw_plan_dft_r2c_1d(
        N,
        in,
        out,
        FFTW_ESTIMATE
    );

    // Exécution
    fftw_execute(plan);

    // Affichage
    for (int i = 0; i < N/2 + 1; i++)
    {
        printf("%d : %f + %fi\n",
               i,
               out[i][0],
               out[i][1]);
    }

    // Libération
    fftw_destroy_plan(plan);
    fftw_free(in);
    fftw_free(out);

    return 0;
}
  • Compiler
ml gcc fftw
gcc fftw_simple.c -o fftw_simple -lfftw3 -lm -o fftw_simple
  • Exécuter dans un job
srun ./fftw_simple

Cas d’usage : Calcul parallèle d’une FFT 1D avec FFTW et MPI

Créer le fichier fftw_mpi.c

#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
#include <fftw3-mpi.h>


int main(int argc, char **argv)
{
    MPI_Init(&argc, &argv);

    fftw_mpi_init();


    int rank, size;

    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);


    ptrdiff_t N = 1024;


    ptrdiff_t local_n0;
    ptrdiff_t local_0_start;
    ptrdiff_t local_n1;
    ptrdiff_t local_1_start;


    ptrdiff_t alloc_local =
        fftw_mpi_local_size_1d(
            N,
            MPI_COMM_WORLD,
            FFTW_FORWARD,
            FFTW_ESTIMATE,
            &local_n0,
            &local_0_start,
            &local_n1,
            &local_1_start
        );


    fftw_complex *data;

    data = fftw_malloc(sizeof(fftw_complex) * alloc_local);


    if (!data)
    {
        fprintf(stderr,
                "Erreur allocation mémoire rang %d\n",
                rank);

        MPI_Finalize();
        return 1;
    }


    /*
       Initialisation locale
    */
    for (ptrdiff_t i = 0; i < local_n0; i++)
    {
        ptrdiff_t global_i = local_0_start + i;

        data[i][0] = (double)global_i;
        data[i][1] = 0.0;
    }


    /*
       Affichage distribution MPI
    */
    MPI_Barrier(MPI_COMM_WORLD);

    for (int r = 0; r < size; r++)
    {
        if (rank == r)
        {
            printf("\nProcessus %d/%d\n", rank, size);

            printf("  Indices globaux : %td -> %td\n",
                   local_0_start,
                   local_0_start + local_n0 - 1);

            printf("  Nombre de points locaux : %td\n",
                   local_n0);

            printf("  Exemple avant FFT :\n");

            for (int i = 0; i < 5 && i < local_n0; i++)
            {
                printf("    data[%d] = %f + %fi\n",
                       i,
                       data[i][0],
                       data[i][1]);
            }

            fflush(stdout);
        }

        MPI_Barrier(MPI_COMM_WORLD);
    }



    /*
       Création plan FFT
    */
    fftw_plan plan =
        fftw_mpi_plan_dft_1d(
            N,
            data,
            data,
            MPI_COMM_WORLD,
            FFTW_FORWARD,
            FFTW_ESTIMATE
        );


    if (!plan)
    {
        printf("Erreur plan FFT rang %d\n", rank);

        fftw_free(data);
        MPI_Finalize();
        return 1;
    }



    double t1 = MPI_Wtime();


    fftw_execute(plan);


    double t2 = MPI_Wtime();



    /*
       Affichage résultat FFT
    */
    MPI_Barrier(MPI_COMM_WORLD);


    for (int r = 0; r < size; r++)
    {
        if (rank == r)
        {
            printf("\nRésultats FFT rang %d\n", rank);

            for (int i = 0; i < 5 && i < local_n0; i++)
            {
                printf("    FFT[%td] = %f + %fi\n",
                       local_0_start+i,
                       data[i][0],
                       data[i][1]);
            }

            printf("Temps calcul local : %f secondes\n",
                   t2-t1);

            fflush(stdout);
        }

        MPI_Barrier(MPI_COMM_WORLD);
    }



    if(rank == 0)
    {
        printf("\nFFT globale terminée avec %d processus MPI\n",
               size);
    }



    fftw_destroy_plan(plan);

    fftw_free(data);


    MPI_Finalize();

    return 0;
}
  • Compiler
ml gcc mpi fftw/3.3.11-intelmpi-2021.16 
mpicc fftw_mpi.c -o fftw_mpi -lfftw3_mpi -lfftw3 -lm
  • Exécuter dans un job
srun --ntasks-per-node=4 mpirun -n 4 ./fftw_mpi

Documentation