Исследования › RefNUMA – библиотека для организации виртуальной общей памяти в программах, использующих MPI.

RefNUMA – библиотека для организации виртуальной общей памяти в программах, использующих MPI.

Пример 8. То же модельное приложение для гибридного суперкомпьютера.

В состав каждого вычислительного узла гибридного суперкомпьютера, помимо процессоров, входит, как минимум, один сопроцессор-ускоритель. Некоторые функции в составе параллельной программы переносятся на сопроцессор с целью ускорения счета. Далее будем говорить, для краткости, о процессорных и сопроцессорных функциях в составе программы. Также введем понятие процессорной и сопроцессорной памяти и, соответственно, процессорных и сопроцессорных адресов, переменных и массивов. Дело в том, что сопроцессор оснащен собственной памятью, в которую должны время от времени копироваться исходные данные, и из которой, соответственно, должны копироваться результаты.

Модель взаимодействия процессорной и сопроцессорной частей программы следующая.

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

Таким образом, простейший цикл использования сопроцессора включает следующие действия:

— процессорная часть программы запрашивает необходимые массивы в памяти сопроцессора, и запоминает предоставленные ей в ответ сопроцессорные адреса этих массивов,

— процессорная часть программы копирует исходные данные во входные сопроцессорные массивы,

— процессорная часть программы вызывает сопроцессорную функцию, передавая ей в качестве аргументов, в частности, указатели на массивы, ранее выделенные в памяти сопроцессора,

— процессорная часть программы копирует результаты счета из выходных сопроцессорных массивов.

Гибридные расширения RefNUMA не ставят задачи реализовать полную программную поддержку какого-либо конкретного сопроцессора-ускорителя. Автоматизируются только некоторые рутинные действия по копированию данных между процессорной и сопроцессорной памятью, которые являются общими для всех сопроцессоров-ускорителей. В механизм собственно вызова сопроцессорной функции RefNUMA никак не вмешивается. Считается, что вызов сопроцессорной функции с передачей аргументов может быть выполнен. Как конкретно – знает программист, и использует для этого соответствующую базовую низкоуровневую технологию, специфичную для конкретного сопроцессора: например, OpenCL для GPGPU.

Автоматизация копирования данных, в свою очередь, затрагивает два аспекта: собственно копирование данных из блочно распределенного разделяемого массива (процессорного) в ускорительный одномерный массив и обратно, и автоматическое выполнение пересылок данных между ускорителями, управляемыми из разных процессов параллельной программы. Остановимся на втором аспекте немного подробнее.

Как мы отмечали выше, обработка данных в сопроцессоре включает в себя копирование входных данных в память сопроцессора, вызов сопроцессорной функции и обратное копирование результатов. Очень часто – например, в нашем модельном приложении, да и вообще в итерационных численных методах, в промежутке между начальной загрузкой исходных данных и конечной выгрузкой результатов выполняется не один, а много шагов обработки. При этом, после каждого шага требуется обмен данными между сопроцессорами разных процессов параллельной программы. Традиционно такой обмен данными записывается довольно громоздко, поскольку посредниками при таком обмене вынуждены выступать процессорные части соответствующих программ. Ведь данные, подлежащие пересылке из сопроцессора в сопроцессор, надо прочитать из сопроцессора, послать процессу – хозяину другого сопроцессора, а там – принять и записать по месту назначения.

Такие «транзитные» обмены данными между сопроцессорами разных процессов, оказывается, при использовании RefNUMA легко выполнить автоматически, конечно, не в самом общем случае, а при некоторых условиях, накладываемых на прикладную программу. Забегая вперед, отметим, что условия эти соблюдаются для нашего модельного приложения, и соответствующая автоматизация в примере из настоящего раздела используется в полной мере.

Рассмотрение текста примера начнем с текста главной программы. Для большей наглядности, все изменения (преимущественно, вставки) текста по сравнению с предыдущим примером, посвященные обслуживанию сопроцессорной части программы, оформлены как условно компилируемые по ключу FORACCEL. Если этот ключ не определен, то компилируется чисто процессорный вариант, совпадающий с предыдущим примером.

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <coarray.h>
#include <shared.h>
#ifdef FORACCEL
#include <accelmem.h>
#include <accelshared.h>
#endif
    COARRAY_MEM_ALLOC( 2000000000l );
    static void usage( void )
{
    printf( "Usage: progrev_coarr <nrows> <ncols> <niter>\n" );
    exit( -1 );
}
    void itstep( int mx, int my, void *pf, void *pnewf, void *pr,
                                                  double rdx2, double rdy2, double beta );
    int main( int argc, char *argv[] )
{
    int mx, my, mx1, niter, n_nodes, my_node;
    void *ra;
    long im, j, k;
    register int i;
    FILE *fp;
    double rdx2, rdy2, beta, t;
#ifdef FORACCEL
    long accelstart, accelfinish, accelsize;
#endif
/***/
    COARRAY_Init( &argc, &argv );
    if ( argc != 4 ) usage();
    rdx2 = 1.0;
    rdy2 = 1.0;
    beta = 0.25;
    my_node = coarray_my_node();
    n_nodes = coarray_n_nodes();
    mx = (int)atol( argv[1] ) + 2;
    my = (int)atol( argv[2] ) + 2;
    niter = (int)atol( argv[3] );
     {
      COARRAY_Create_vectorized( double, f, mx, my, 1.0, fo );
      coarray_set_name( (void**)fo, "Array_f" );
      COARRAY_Create_vectorized( double, newf, mx, my, 1.0, newfo );
      coarray_set_name( (void**)newfo, "Array_newf" );
      COARRAY_Create_vectorized( double, r, mx, my, 1.0, ro );
      coarray_set_name( (void**)ro, "Array_r" );
#ifdef FORACCEL
      void *af, *anewf, *ar;
#endif
      if ( (!f) || (!newf) || (!r) )
       {
        printf( "No memory\n" );
        exit( -1 );
       }
#ifdef FORACCEL
      accelstart = SHARED_START( 1, mx );
      accelfinish = SHARED_FINISH( mx-1, mx );
      accelsize = accelfinish-accelstart;
      if ( accelsize <= 0 )
       {
        accelsize = 0;
       }
      else
       {
        accelstart--;
        accelfinish++;
        accelsize += 2;
       }
//      printf( "Process %d: accelstart=%ld accelfinish=%ld accelsize=%ld\n", my_node,
                     accelstart, accelfinish, accelsize );
      if ( accelsize )
       {
        af = calloc_accel( accelsize, my*sizeof(double) );
        anewf = calloc_accel( accelsize, my*sizeof(double) );
        ar = calloc_accel( accelsize, my*sizeof(double) );
        if ( (!af) || (!anewf) || (!ar) )
         {
          printf( "No device memory\n" );
          exit( -1 );
         }
       }
#endif
/***/
      PARALLEL_DO( i, 0, mx, 1, mx )
       {
        for ( j = 0; j < my; j++ )
         {
          f[i][j] = newf[i][j] = 0.0;
          r[i][j] = 0.0;
          if ( (i == 0) || (i == (mx-1)) || (j == 0) || (j == (my-1)) )
                                                               f[i][j] = newf[i][j] = 1.0;
         }
       }
#ifdef FORACCEL
      vectorized_to_accel( af, (void **)f, mx, my, accelstart, accelsize,
                                                             sizeof(double), REALLYCOPY );
      vectorized_to_accel( anewf, (void **)newf, mx, my, accelstart, accelsize,
                                                             sizeof(double), REALLYCOPY );
      vectorized_to_accel( ar, (void **)r, mx, my, accelstart, accelsize, sizeof(double),
                                                             sizeof(double), REALLYCOPY );
#endif
      coarray_barrier();
      t = MPI_Wtime();
      for ( k = 0; k < niter; k++ )
       {
#ifdef FORACCEL
	if ( accelsize )
	 {
	  vectorized_to_accel( af, (void **)f, mx, my, accelstart, accelsize,
                                                 sizeof(double), JUSTUPDATE );
	  itstep( accelsize, my, af, anewf, ar, rdx2, rdy2, beta );
	  vectorized_from_accel( (void **)newf, mx, my, accelstart, anewf, accelsize,
                                                       sizeof(double), JUSTUPDATE );
         }
#else
	itstep( mx, my, (void*)f, (void*)newf, (void*)r, rdx2, rdy2, beta );
#endif
        coarray_barrier();
        ra = (void*)f;
        f = newf;
        newf = (typeof(newf))ra;
#ifdef FORACCEL
        ra = (void*)af;
        af = anewf;
        anewf = (typeof(anewf))ra;
#endif
       }
      t = MPI_Wtime() - t;
      printf( "Time: %f speed %f flops\n", t,
                  (((double)(mx-2))*((double)(my-2))*7.0*niter)/t );
#ifdef FORACCEL
      vectorized_from_accel( (void **)f, mx, my, accelstart+1, ((double*)af)+my,
                                              accelsize-2, sizeof(double), REALLYCOPY );
#endif
      if ( !my_node )
       {
        fp = fopen( "progrev_coarr.dat", "w" );
        fclose( fp );
       }
      NODE_BY_NODE_BEGIN( k, 0, n_nodes )
      fp = fopen( "progrev_coarr.dat", "a" );
      PARALLEL_DO( i, 1, mx-1, 1, mx )
       {
//        fwrite( f[i]+1, sizeof(double), my-2, fp );
	for ( j = 1; j < my-1; j++ ) fprintf( fp, "%f  ", f[i][j] ); fprintf( fp, "\n" );
       }
      fclose( fp );
      coarray_report( (void**)fo, &im, &j );
      printf( "In node %d for %s: requested %ld used %ld\n", my_node,
                   coarray_get_name((void**)fo), im, j );
      coarray_report( (void**)newfo, &im, &j );
      printf( "In node %d for %s: requested %ld used %ld\n", my_node,
                  coarray_get_name((void**)newfo), im, j );
      coarray_report( (void**)ro, &im, &j );
      printf( "In node %d for %s: requested %ld used %ld\n", my_node,
                  coarray_get_name((void**)ro), im, j );
      fflush( stdout );
      NODE_BY_NODE_END
     }
    COARRAY_Finalize();
    return 0;
}
#ifndef FORACCEL
    void itstep( int mx, int my, void *pf, void *pnewf, void *pr, double rdx2, double rdy2,
                        double beta )
{
    int mx1, my1;
    register int i, j;
/***/
    double **f = (typeof(f))pf;
    double **newf = (typeof(newf))pnewf;
    double **r = (typeof(r))pr;
    mx1 = mx - 1;
    my1 = my - 1;
    PARALLEL_DO( i, 1, mx1, 1, mx )
     {
      for ( j = 1; j < my1; j++ )
       {
        newf[i][j] = ((f[i-1][j]+f[i+1][j])*rdx2+(f[i][j-1]+f[i][j+1])*rdy2-r[i][j])*beta;
       }
     }
}
#endif

Поясним кратко смысл фрагментов приведенного текста, начинающихся со строки

#ifdef FORACCEL

Файлы заголовков accelmem.h и accelshared.h обязательны для включения в программы, использующие гибридные расширения RefNUMA.

Помимо хорошо знакомых нам массивов общей памяти f, newf и r, вводятся указатели af, anewf и ar. Это указатели на массивы в памяти сопроцессора. В них хранятся копии тех частей массивов f, newf и r, к которым осуществляет доступ данный процесс.

Фрагмент, начинающийся с обращения к SHARED_START(), выясняет, сколько модифицируемых (таких, что процесс что-то им присваивает) строк распределенных массивов f и newf досталось данному процессу. Каждый процесс в нашей программе модифицирует свои строки, за исключением нулевой и последней строк массива. Из-за того, что число строк в массиве может не делиться нацело на число процессов, и размер блока строк, принадлежащего процессу, может поэтому быть взят с запасом, последнему (или последним) по порядку номеров процессу (процессам) может вообще не достаться обрабатываемых строк. Если в чисто процессорном, не гибридном варианте программы этот случай отрабатывался автоматически, то в нашем, гибридном, варианте такую ситуацию в некоторых местах программы приходится учитывать явно. В дальнейшем значение переменной accelsize равно либо нулю, если это последний процесс, которому строк не досталось, либо числу строк распределенного массива, подлежащих копированию в сопроцессор. Чтобы выяснить это число, надо к числу модифицируемых (то есть своих для данного процесса, не нулевых и не последних) строк прибавить 2, поскольку в сопроцессор копируются модифицируемые строки, одна строка перед ними и одна строка после них.

Зная вычисленное таким образом значение accelsize, можно запросить память в сопроцессоре для массивов af, anewf и ar. Память запрашивается обращением к системной функции calloc_accel(), первый аргумент которой – размер требуемой памяти в строках распределенного массива, второй – размер строки в байтах.

После этого можно, при помощи системной функции vectorized_to_accel(), скопировать в сопроцессор начальное состояние необходимой этому сопроцессору части распределенного массива. Длина этой части массива – accelsize строк. При копировании из процессорного распределенного массива в сопроцессорный одномерный массив, копируемая часть распределенного массива укладывается в сопроцессорный одномерный массив плотно, строка за строкой.

Вызов функции:

vectorized_to_accel( af, (void **)f, mx, my, accelstart, accelsize,
                                                             sizeof(double), REALLYCOPY );

В итерационном цикле, как и в процессорном (не гибридном) варианте примера, вызывается функция itstep(). Текст этой функции в гибридном варианте находится в отдельном файле, который мы рассмотрим немного позже. Вызов функции получает в качестве аргументов, в том числе, и сопроцессорные массивы, из чего можно сделать вывод (правильный), что функция эта внутри себя вызывает некоторую сопроцессорную функцию, которая эти сопроцессорные массивы обрабатывает.

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

vectorized_from_accel( (void **)newf, mx, my, accelstart, anewf, accelsize,
                                                             sizeof(double), REALLYCOPY );

Однако, в обоих вызовах функций копирования, которые мы сейчас обсуждаем, последний аргумент задан не константой REALLYCOPY, как при первоначальном копировании исходных данных, а константой JUSTUPDATE. Кроме того, несколько настораживает тот факт, что при каждой итерации выполняется копирование входного массива в сопроцессор (перед итерацией), и копирование выходного массива из сопроцессора (после итерации). Казалось бы, по смыслу алгоритма достаточно копировать только края обрабатываемой каждым сопроцессором порции данных, точнее, читать из сопроцессора первую и предпоследнюю строки, а записывать нулевую и последнюю. Связаны ли между собой эти две «странности» в тексте программы? Безусловно.

Смысл копирования данных между процессором и сопроцессором перед каждой итерацией и после нее состоит в организации передачи крайних строк обрабатываемых данных между сопроцессорами разных процессов. Передачу процессорных данных между разными процессами RefNUMA выполняет автоматически – для этого она и предназначена. Однако, на сопроцессорную память ее (RefNUMы) способности не распространяются. Казалось бы, тогда следует потребовать от программиста, чтобы он вычислил в явном виде, какие именно строки распределенного массива, изменившиеся на данной итерации, понадобятся другим процессам, и скопировал их из сопроцессора в распределенный массив. Затем следует выполнить барьер, после чего скопировать обратно в сопроцессор именно те чужие строки процессорной версии массива, которые изменились в результате барьера. Легко понять, что возложение на программиста всех этих нетривиальных и кропотливых манипуляций губит идею виртуальной общей памяти на корню – ведь мы так старались уйти от понятия явной передачи данных из процесса в процесс, а теперь вынуждены вновь к нему вернуться. Нельзя ли «попросить» RefNUMу самостоятельно сообразить, какие именно части массивов f и newf куда и как копировать до и после очередной итерации? Можно и нужно. Именно это и делается при задании последнего параметра функций копирования равным предопределенной системной константе JUSTUPDATE. При таком значении последнего параметра происходит копирование не всего указанного объема данных, а только тех его частей, которые у разных сопроцессоров перекрываются.

Насколько произвольной по организации коммуникационной логики может быть программа, для которой RefNUMA в состоянии будет правильно вычислить подлежащие копированию области перекрытия? Требований к программе, для которой описанная только что автоматизация будет работать, всего два, и в нашей модельной программе они оба выполняются.

1). В сопроцессоре поддерживается копия строго определенной для данного процесса части каждого из используемых в нем массивов общей памяти. При каждом вызове любой сопроцессорной функции каждая такая копия части массива общей памяти, переданная в качестве аргумента, должна быть, относительно данного вызова, строго входным или строго выходным параметром. Параметры типа «входной, модифицируемый» не допускаются,

2). Каждый процесс параллельной программы (в совокупности, включая как процессорную, так и сопроцессорную часть) должен записывать только свои данные блочно распределенных разделяемых массивов, а чужие – только читать.

В заключение рассмотрим часть главной программы, в которой происходит запись результатов в файл. Вплоть до этого места программы, мы всякий раз, когда нам требовалось скопировать какие-то данные между процессорной и сопроцессорной памятью, указывали один и тот же диапазон копируемых строк, раз и навсегда посчитанный в начале программы. Однако, при выполнении итогового копирования результатов счета из сопроцессора, мы указали отдельно вычисленный, более узкий диапазон копируемых строк. Почему? Это – довольно тонкое место, которое лучше всего понять, но, в крайнем случае, можно просто запомнить.

Диапазоны копируемых строк, которыми мы пользовались до сих пор, в разных процессах перекрываются, поскольку каждый процесс, в дополнение к своим строкам, записывает в сопроцессор еще одну строку перед своими, и одну – после. Используя для итоговой выгрузки из сопроцессора результатов счета именно такой диапазон строк, мы будем загружать перекрывающиеся строки в распределенный массив дважды, в разных процессах. Можно показать, что в тех процессах, где загружаемые из сопроцессора строки являются чужими, будут загружены значения с прошлой итерации. В итоге распределенный массив окажется местами испорченным. Конечно, выстроив правильную комбинацию из нескольких вызовов vectorized_from_accel(), некоторые из которых будут в режиме JUSTUPDATE, значения загруженных перекрывающихся строк можно синхронизировать. Однако, объяснение и обоснование такой программы настолько многословно и вычурно, что здесь не приводится. Много проще написать программу так, чтобы каждая строка считывалась из сопроцессора один раз, из гарантированно правильного места, что мы и сделали в нашем примере.

Перейдем к рассмотрению второго файла с исходным текстом данной программы, в котором содержится гибридный вариант функции itstep(). Первое, что нам надлежит сделать при таком рассмотрении – это определиться, какой именно конкретный сопроцессор мы используем. Решение вопроса с конкретным сопроцессором состоит из двух частей:

1). На каком языке пишутся, как транслируются и вызываются сопроцессорные функции?

2). Как рассмотренные нами выше функции работы с памятью сопроцессора настраиваются на специфику конкретного сопроцессора?

На первый вопрос ответить проще всего. RefNUMA, даже в гибридном варианте, этим вопросом не занимается вообще. Программа с использованием гибридных расширений RefNUMA обязательно использует некоторую базовую технологию доступа к конкретному сопроцессору для запуска сопроцессорных программ. В нашем примере, в качестве сопроцессора используется устройство с архитектурой CUDA, и соответствующие инструментальные средства фирмы nVidia. Если бы использовался другой вид GPGPU, или вообще другой вид сопроцессора, пришлось бы применять другую базовую технологию. В любом случае, конкретный ответ на этот вопрос ортогонален к возможностям RefNUMA. Интересным частным случаем, который в нашем примере не рассматривается, является случай, когда сопроцессор вообще не существует физически, а использование гибридных расширений RefNUMA является просто приемом программной инженерии. Например, при переносе гибридной программы на суперкомпьютер из обычных многоядерных процессов общего назначения может оказаться удобным организовать «программный сопроцессор», написанный на OpenMP, который сможет эффективно использовать все ядра вычислительного узла для выполнения вычислительно емких фрагментов программы. Даже одноядерный «программный сопроцессор» может иметь некоторый смысл за счет того, что внутри сопроцессора обрабатываемые данные не имеют формы распределенного массива, а лежат в программно выделенной такому «сопроцессору» памяти плотно, и могут быть обработаны оптимизированными библиотечными функциями.

Ответ на второй вопрос несколько сложнее. В нашем примере в качестве сопроцессора используется устройство с архитектурой CUDA, и соответствующие базовые инструментальные средства. Значит, функция vectorized_to_accel() должна, так или иначе, вызвать функцию cudaMemcpy(). Если так, то что делать в случае использования вместо CUDA, например, OpenCL?

Вопрос решается путем введения базовых функций доступа к сопроцессору. Это функции, которые вызываются из функций и макросов гибридных расширений RefNUMA, имеющие совершенно определенный внешний интерфейс. Программист должен написать эти функции сам, опираясь на используемые им базовые инструментальные средства, и включить в состав своей программы. Для другого типа ускорителя их придется переписать. Программа пользователя не обязана (хотя может) вызывать эти функции напрямую. В нашем примере она этого не делает.

Ниже приводится исходный текст базовых функций доступа к сопроцессору, написанных с использованием базовых инструментальных средств технологии CUDA.

	void to_accel( void *acceladdr, void *cpuaddr, long bleng )
{
	cudaMemcpy( acceladdr, cpuaddr, bleng, cudaMemcpyHostToDevice );
}
	void from_accel( void *cpuaddr, void *acceladdr, long bleng )
{
	cudaMemcpy( cpuaddr, acceladdr, bleng, cudaMemcpyDeviceToHost );
}
	void *malloc_accel( long bleng )
{
	void *ret;
	cudaMalloc( &ret, bleng );
	return ret;
}

Функция

to_accel( void *acceladdr, void *cpuaddr, long bleng ) 

Функция

from_accel( void *cpuaddr, void *acceladdr, long bleng )

Функция

void *malloc_accel( long bleng )

Все упомянутые массивы являются сплошными в соответствующей памяти.

Теперь осталось рассмотреть текст функции itstep() в ее конкретном варианте для CUDA. Объяснение технологии CUDA, даже в ее простейшем варианте, не входит в задачу настоящего документа. Полный исходный текст файла приводится ниже.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <accelmem.h>
#define GRIDSIZE 32
#define BLOCKSIZE 1024
    extern __shared__ double pmiddle[];
    __global__ void iteration( double *f, double *fn, double *rhs, int nrows, int ncols,
			       double rdx2, double rdy2, double beta2 )
{
    int i, n;
    double *pf, *pfn, *prhs;
    double *middle = pmiddle;
    double upper, lower;
/***/
/* calculate this grid start & size: */
    n = nrows/gridDim.x;
    pf = f + n*ncols*blockIdx.x;
    pfn = fn + n*ncols*blockIdx.x;
    prhs = rhs + n*ncols*blockIdx.x;
    if ( blockIdx.x == (gridDim.x-1) ) n = nrows - n*(gridDim.x-1);
/* make the shadow edges: */
    if ( blockIdx.x == 0 )
     {
      n++;
     }
    else
     {
      pf -= ncols;
      pfn -= ncols;
      prhs -= ncols;
      if ( blockIdx.x == (gridDim.x-1) )
       n++;
      else
       n += 2;
     }
/* pipeline the grid portion through the shared memory "vector registers": */
    upper = pf[threadIdx.x];
    pf += ncols;
    pfn += ncols;
    middle[threadIdx.x] = pf[threadIdx.x];
    prhs += ncols;
    for ( i = 2; i < n; i++ )
     {
      lower = pf[threadIdx.x+ncols];
      __syncthreads();
      if ( (threadIdx.x > 0) && (threadIdx.x < (blockDim.x-1)) )
       {
        pfn[threadIdx.x] = beta2*(
    	       rdx2*(upper+lower)
	     + rdy2*(middle[threadIdx.x-1]+middle[threadIdx.x+1])
	     - prhs[threadIdx.x]);
       }
      prhs += ncols;	
      pf += ncols;
      pfn += ncols;
      __syncthreads();
      upper = middle[threadIdx.x];
      middle[threadIdx.x] = lower;
     }
}
       	void itstep( int nrows, int ncols, void *f, void *fn, void *r,
                              double rdx, double rdy, double beta )
{
	int n, nl;
/***/
	n = nrows/GRIDSIZE;
	nl = nrows - n*(GRIDSIZE-1);
	if ( (n < 1) || (nl < 1) )
	 {
	  fprintf( stderr,
   "Too many grid blocks (%d) for the number of rows per device (%d), exiting\n",
              GRIDSIZE, nrows );
	  exit( -1 );
	 }	
	n = ncols;
	while ( n > 2 )
	 {
	  nl = n;
	  if ( nl > BLOCKSIZE ) nl = BLOCKSIZE;
	  iteration<<<GRIDSIZE,nl,nl*sizeof(double)>>>(
                      (double*)f, (double*)fn, (double*)r, nrows, ncols, rdx, rdy, beta
                                                                                        );
	  f = ((double*)f)+(nl-2);
	  fn = ((double*)fn)+(nl-2);
	  r = ((double*)r)+(nl-2);
	  n -= (nl-2);
	 }
}
	void to_accel( void *acceladdr, void *cpuaddr, long bleng )
{
	cudaMemcpy( acceladdr, cpuaddr, bleng, cudaMemcpyHostToDevice );
}
	void from_accel( void *cpuaddr, void *acceladdr, long bleng )
{
	cudaMemcpy( cpuaddr, acceladdr, bleng, cudaMemcpyDeviceToHost );
}
	void *malloc_accel( long bleng )
{
	void *ret;
	cudaMalloc( &ret, bleng );
	return ret;
}

Подробное объяснение текста функций itstep() и iteration() можно прочитать в документе «Десять простых шагов к стройной прикладной программе для гибридного суперкомпьютера», в котором рассматривается пример гибридного модельного приложения для CUDA, очень похожий в своей сопроцессорной части на данный пример.

В заключение добавим, что, с целью упрощения процедуры трансляции и сборки гибридных суперкомпьютерных приложений, настоятельно рекомендуется разносить использование технологий межузловых коммуникаций (RefNUMA, MPI, shmem и т. п.), с одной стороны, и базовых технологий доступа к сопроцессору (CUDA, OpenCL), с другой стороны, по разным файлам с исходными текстами. В приведенном только что примере это сделано именно так.