Основы MPI

Что такое MPI? Смотрим википедию. Здесь же мы придержимся нашего жанра "шпоргалки".
Шесть "главных" MPI команд:

• MPI_Init(int *argc, char ***argv)
• MPI_Finalize()
• MPI_Comm_size(MPI_Comm comm, int *size)
• MPI_Comm_rank(MPI_Comm comm, int *rank)
• MPI_Get_processor_name(char *name, int *resultlen)
• MPI_Send(void* buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm)
• MPI Recv(void* buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Status *status)

Основные понятия MPI:

• типы MPI:
  MPI_CHAR, MPI_SHORT, MPI_INT, MPI_LONG, MPI_UNSIGNED_CHAR, MPI_UNSIGNED_SHORT, MPI_UNSIGNED, MPI_UNSIGNED_LONG, MPI_FLOAT, MPI_DOUBLE, MPI_LONG_DOUBLE
• MPI коммуникатор
• MPI tag
• MPI Processes

Проиллюстрируем это обилие информации классическим приветствием миру. Создаем файл helloworld.c следующего содержания:

#include
#include
int main(int argc, char ** argv) {
    int size,rank;
    int length;
    char name[80];

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD,&rank);
    MPI_Comm_size(MPI_COMM_WORLD,&size);
    MPI_Get_processor_name(name,&length);
    printf("Hello MPI! Process %d of %d on %s\n",size,rank,name);
    MPI_Finalize();
    return 0;
}

Компилируем и запускаем на выполнение:
#mpicc helloworld.c -o helloworld
#mpirun -np 2 ./helloworld
Запустить на выполнение можно так:
#mpirun -np 2 -machinefile hostfile ./helloworld
где hostfile - файл с таким, например, содержанием:
aa slots=4
bb slots=4
cc slots=4
(взято с мануала к команде mpirun) aa, bb, cc - имена узлов. Поварируйте число после ключа -np, узнаете много интересного.
Рассмотрим передачу сообщений между процессами. Отредактируем файл helloworld.c следующим образом:

#include
#include
int main(int argc, char ** argv) {
    int size,rank;
    int length;
    char name[80];
    int dest = 0;
    int tag = 999;

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD,&rank);
    MPI_Comm_size(MPI_COMM_WORLD,&size);
    MPI_Get_processor_name(name,&length);
    printf("Hello MPI! Process %d of %d on %s\n",size,rank,name);
    if (rank != 0 ) { /* I’m a client */
        MPI_Send(name,80,MPI_CHAR,dest,tag,MPI_COMM_WORLD);
    }else{ /* I’m the server (rank == 0) */
        MPI_Status status;
        int source;
        for(source = 1; source < size; source++) {
 MPI_Recv(name,80,MPI_CHAR,source,tag,MPI_COMM_WORLD,&status);
            printf(" mesg from %d of %d on %s\n",source,size,name);
        }
    }   
    MPI_Finalize();
    return 0;
}

Компилируем и запускаем на выполнение теми же командами, что и первый пример.
Оба примера компилировались и тестировались на компьютере с Debian Wheezy:
#uname -a
Linux myhost 3.2.0-4-686-pae #1 SMP Debian 3.2.39-2 i686 GNU/Linux

Дискретное преобразование Фурье в fftw3 и в CUDA

Если вам вдруг понадобилось вычислить дискретное преобразование Фурье, то для этой цели прекрасно подойдут два инструмента: http://www.fftw.org/ и CUDA. К достоинствам первого можно отнести его открытость, специализация именно на преобразовании Фурье, предельная простота установки. 

Рассмотрим для начала пример с использованием fftw.

Создадим файл fftw3.c следующего содержания:

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

int
main ( void )
{
    int i;
    fftw_complex *in;
    fftw_complex *in2;
    int n = 10;
    fftw_complex *out;
    fftw_plan plan_backward;
    fftw_plan plan_forward;

    //Создаем входные данные
    in = fftw_malloc ( sizeof ( fftw_complex ) * n );
    printf ( "Input Data:\n" );
    for ( i = 0; i < n; i++ ){
        in[i][0] = i;in[i][1] = 0;
        printf ( "  %3d  %12f  %12f\n", i, in[i][0], in[i][1] );
    }

    //Выполняем прямое преобразование Фурье
    out = fftw_malloc ( sizeof ( fftw_complex ) * n );
    plan_forward = fftw_plan_dft_1d ( n, in, out, FFTW_FORWARD, FFTW_ESTIMATE );
    fftw_execute ( plan_forward );

    printf ( "Output FFT Coefficients:\n" );
    for ( i = 0; i < n; i++ )
        printf ( "  %3d  %12f  %12f\n", i, out[i][0], out[i][1] );

    //Выполняем обратное преобразование Фурье
    in2 = fftw_malloc ( sizeof ( fftw_complex ) * n );
    plan_backward = fftw_plan_dft_1d ( n, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE );
    fftw_execute ( plan_backward );

    printf ( "Recovered input data:\n" );
    for ( i = 0; i < n; i++ )
        printf ( "  %3d  %12f  %12f\n", i, in2[i][0]/(double)(n), in2[i][1]/(double)(n) );

    //Прибираем за собой
    fftw_destroy_plan ( plan_forward );
    fftw_destroy_plan ( plan_backward );
    fftw_free ( in );
    fftw_free ( in2 );
    fftw_free ( out );

  return 0;
}

Компилируем его с помощью следующей команды:

#gcc fftw3.c -o fftw3 -lfftw3 -lm

 Как видите, если убрать код, выводящий данные на стандартный вывод, то остальной код, выполняющий главную функциональность, будет смотреться предельно просто и ясно. Всего несколько строчек. Для завершения картины запускаем программу:

#./fftw3

и изучаем вывод. Должны увидеть, что обратное преобразование восстанавливает входные данные для прямого преобразования.

Далее рассмотрим решение этой же задачи средствами CUDA. К недостаткам этого метода можно отнести не всегда простую установку CUDA  и обязательное наличие видеокарты, поддерживающую эту технологию. К достоинствам отнесем всю мощь данной технологии. Код взят из примеров, поставляемых вместе с тулкитом с сайта NVIDIA. Называется он (код, с помощью которого была написана программа, исходники которой далее приводятся) иронично simpleCUFFT. Хотя простоты там, на мой взгляд нет совсем, потому что пример, показывающий как пользоваться преобразованием Фурье, завален кодом, делающим конволюцию сигнала с фильтром (зачем это людям, изучающим CUDA?), и выковырять пару строчек, которые и делают прямую задачу примера, было для меня не просто.  За деревьями леса не видно. Поэтому, думаю, этот пост будет полезен.

Итак, создаем файл truelySimpleFFT.cu следующего содержания:

#include<stdlib.h>
#include<stdio.h>
#include<string.h>
#include<math.h>

#include <cuda_runtime.h>
#include <cufft.h>

typedef float2 Complex;
#define N_ELEM 10
int
main(int argc, char **argv)
{
    Complex *h_signal;
    Complex *d_signal;
    int mem_size;
    cufftHandle plan;
   
    printf("[truely simpleCUFFT] is starting...\n");
    //Выделяем память для входных данных на хосте
    h_signal = (Complex *)malloc(sizeof(Complex) * N_ELEM);

    //Инициализируем входные данные
    for (unsigned int i = 0; i < N_ELEM; ++i){
        h_signal[i].x = (float)i;
        h_signal[i].y = 0;
        printf("%f\t%f\n", h_signal[i].x, h_signal[i].y);
    }  
   
    //Выделяем память для входных данных на видеокарте и копируем их туда
    mem_size = sizeof(Complex) * N_ELEM;     
    cudaMalloc((void **)&d_signal, mem_size);                                  
    cudaMemcpy(d_signal, h_signal, mem_size, cudaMemcpyHostToDevice);                 
     
    //Собственно выполняем преобразование Фурье
    printf("Transforming signal cufftExecC2C\n");
    cufftPlan1d(&plan, N_ELEM, CUFFT_C2C, 1);                        
    cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_FORWARD);
   
    //Копируем результат вычисления с видеокарты на хост
    cudaMemcpy(h_signal, d_signal, mem_size, cudaMemcpyDeviceToHost);
    for (unsigned int i = 0; i < N_ELEM; ++i)
        printf("%f\t%f\n", h_signal[i].x, h_signal[i].y);
   
    //Прибираем за собой
    cufftDestroy(plan);
    free(h_signal);
    cudaFree(d_signal);
    cudaDeviceReset();
   
    exit(0);
}

 

Компилируем его с помощью следующих двух команд:

#nvcc  -o simpleCUFFT.o -c truelySimpleFFT.cu

#g++ -o simpleCUFFT simpleCUFFT.o -lcudart -lcudart -lcufft
Необходимо отметить следующее. Если вы устанавливали SDK с сайта NVIDIA "руками", то вам надо вставить в выше приведенные команды ключи -L и -I с соответствующими путями к заголовочным файлам и библиотекам, которые вы указывали при установке (а также nvcc с полным путем). Я устанавливал с помощью репозитария (пользуюсь Дебиан) и в моем случае надо вводить команды как написано выше.

Также обратите внимание на подобие кода, выполняющего преобразование Фурье в обоих примерах.

Запускаем на выполнение командой:

 #./simpleCUFFT

Обе программы компиллировались и ваполнялись на Debian GNU/Linux 7.0 (wheezy).

Простейшая генерация случайных чисел в языке Си

Приведем простейшую функцию генерации случайных чисел в языке Си в диапазоне от нуля до единицы:

double
frand ( void ){
  double value;
  value = ( ( double ) rand ( ) / ( RAND_MAX ) );
  return value;
}

У Кернигана и Ричи используется конструкция ( ( double ) rand ( ) / ( RAND_MAX + 1 ) ), но она может привести к отрицательным величинам у случайных чисел.

Простейшая работа со временем в языке Си

Допустим нам необходимо получить строку выводящую время на "человеческом" языке. Как исходный шаблон можем использовать следующую функцию:

void
timestamp ( void )
{
    static char time_buffer[100];
    const struct tm *tm;
    size_t len;
    time_t now;

    now = time ( NULL );
    tm = localtime ( &now );
    len = strftime ( time_buffer, 100, "%d %B %Y %I:%M:%S %p", tm );
    printf ( "%s\n", time_buffer );

    return;
}

Пусть теперь нам надо посмотреть, как много времени забирает выполнение того или иного участка кода. Для решения этой задачи будет удобным следующий пример:

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
int
main(int argc, char* argv[], char* envv[]){
        clock_t startClock, finishClock;
        double diffClock;
        startClock = clock();
        ...//Здесь находится кусок кода, время выполнения которого мы хотим узнать
        finishClock = clock();
        diffClock = ((double)(finishClock - startClock))/CLOCKS_PER_SEC;
        printf("Passed time is: %lfsec\n", diffClock);
        exit(0);
}

Расширения языка Си в CUDA

Имеется четыре типа расширений языка Си в CUDA. Это:

1. Квалификаторы типа функции, определяющие, где эта функция будет выполняться - в центральном процессоре или на видеокарте, и возможно ли вызывать эту функцию из центрального процессора или видеокарты.
2. Квалификаторы типа переменных, определяющие распределение памяти для этих переменных на видеокарте.
3. Директива, определяющая как ядро выполняется на видеокарте.
4. Четыре встроенных переменных, определяющих размерность грида, размерность блока, индекс блока и индекс нити.

Рассмотрим эти расширения подробней.

I.Квалификаторы типа функций

1. __device__ квалификатор говорит, что функция выполняется на видеокарте и может быть вызвана только из видеокарты.
2. __global__ квалификатор объявляет функцию выполняемой в  центральном процессоре и вызываемой только процессором. Квалифицирование функции только __host__ аналогично объявлению функции без какого-то ни было квалификатора. Однако функцию можно квалифицировать комбинацией __host__ __device__. Тогда функция будет откомпилирована и для центрального процессора и для видеокарты.  

 Ограничения:

•  функции квалифицированные __device__ и __global__ не поддерживают рекурсии.
•  функции квалифицированные __device__ и __global__ не могут декларировать внутри себя статические переменные.
• функции квалифицированные __device__ и __global__ не могут иметь переменное количество аргументов.
• функции квалифицированные __device__ не могут вызываться по их адресу.
• функции не могут быть одновременно квалифицированы как __global__ и __host__.
• функции квалифицированные как __global__ обязаны возвращать пустой тип.
• функции квалифицированные как __global__ должны специфицировать  свою конфигурацию выполнения.
• функции квалифицированные как __global__ выполняются синхронно.
• параметры функций квалифицированных как __global__ передаются через общую память в видеокарту и ограничены 256-ю байтами.

II. Квалификаторы переменных.

1. __device__ квалификатор говорит, что переменная размещается на видеокарте. Может использоваться совместно с другими квалификаторами типа переменных. Если таковых нет, то это означает, что переменная размещена в глобальной памяти, имеет время жизни равное времени жизни приложения и доступна из нитей в гриде и с центрального процессора через библиотеку времени исполнения.
2. __constant__ квалификатор опционально используемый предыдущим квалификатором объявляет, что переменная размещается в области постоянной памяти, имеет время жизни равное времени жизни приложения и доступна из нитей в гриде и с центрального процессора через библиотеку времени исполнения.
3. __shared__ квалификатор опционально используемый с  __device__ квалификатором объявляет, что перменная размещается в общей памяти блока нитей, имеет время жизни равное времени жизни блока и доступна только из нитей данного блока.

Ограничения:

• Эти квалификаторы не допустимы для членов struct и union, в формальных параметрах функций и локальныхперменных внутри функций, что выполняются на центральном процессоре.
• __shared__ и __shared__ не могут использоваться вместе.
• __shared__ и __shared__ переменные подразумеваются статическими.
• __constant__ переменные не могут использоваться внутри видеокарты, только на центральном процессоре. Поэтому они имеют область видимости файла.
• __shared__ переменные не могут иметь инициализацию как часть их декларации.

III. Конфигурация выполнения.

Любой вызов функции __global__ должен специфицировать конфигурацию выполнения этой функции. Конфигурация выполнения функции определяет размерность грида и блоков, которые будут использоваться при вызове функции на видеокарте. Эта конфигурация задается аргументами между тройными угловыми скобками. <<>>. Эта конструкция помещается между именем функци и списком передаваемым в функцию списком "обычных" аргументов.

Описание аргументов в угловых скобках следующее:

1. Dg - имеет тип dim3 и определяет размерность и длину грида так что выражение Dg.x*Dg.y равно числу запущенных блоков.
2. Db - имеет тип dim3 и определяет размерность и длину каждого блока, так что выражение Db.x*Db.y*Db.z равно числу нитей в каждом блоке.
3. Ns - имеет тип size_t и определяет количество байтов в общей памяти, которая выделяется динамически при данном вызове функции дополнительно к ее статической памяти. Ns - необязательный аргумент со значением по умолчанию раным нулю.

Аргументы конфигурации выполнения вычисляются раньше  "обычных" аргументов функции.

IV. Встроенные переменные

1. gridDim - переменная типа dim3. Содержит размерность грида.
2. blockIdx - переменная типа uint3. Содержит индекс блока в гриде.
3. blockDim - переменная типа dim3. Содержит размерность блока.
4. threadIdx - переменная типа uint3. Содержит индекс нити в блоке.

Ограничения:

• не допустима операция получения адреса от любой встроенной переменной.
• не допустима операция присваивания для любой встроенной переменной.