Главная
Вычислительные ресурсы
C чего начать
Вопросы - ответы
Документация
Исследования
Контакты

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

Более обстоятельные рассуждения показывают, что так поступать нельзя. Блоки обрабатываются последовательно, и содержат перекрытия. Обработав блок номер i, мы разместили в его краях новые (вычисленные на этой итерации) значения, в то время как в область перекрытия блока номер i+1 следовало бы, согласно итерационному процессу Якоби, поместить старые значения. Но старые значения краев блока номер i к моменту обработки блока номер i+1 уже утрачены, взять их неоткуда.

Здесь нам впору было бы пожалеть о том, что при прошлом преобразовании программы мы избавились от двух массивов значений температур, старого и нового, заменив их одним. Эта проблема при наличии двух массивов была бы решена автоматически. При еще более внимательном рассмотрении, впрочем, можно обнаружить, что жалеть нам не о чем, поскольку перед нами стоит еще одна похожая, но более сложная, проблема.

Вспомним, что нарезка обрабатываемых данных на мелкие блоки нам нужна не сама по себе, а в связи с возможностью обрабатывать каждый блок многократно (N > 1 раз). Это, в свою очередь, подразумевает расширение областей перекрытия в N раз. Однако, при N-кратной обработке в расширенных областях перекрытия обрабатываемого блока образуется "мусор", который после каждой N-кратной обработки блока должен быть отброшен. Если мы позволим этому мусору образовываться поверх краев соседних блоков, прямо в массиве новых значений, данные будут на первом же мелком блоке необратимо испорчены.

Значит, преобразование программы, которое нам сейчас предстоит, сложнее, чем кажется. Нам придется не просто реализовать нарезку блоков с перекрытиями, расширенными в N раз, но и изготовить программную модель внутренней памяти сопроцессора, в которую (и из которой) блоки обрабатываемых данных будут копироваться, с учетом того, что в области перекрытия перед началом обработки блока должны быть старые значения, а после конца обработки получившийся на их месте "мусор" не следует копировать обратно, в основной массив программы.

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

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

#include <stdio.h>
#include "dimension.h"
#define PORTION 15
    extern void itstep_accel( int mx, int my, double rdx2, double rdy2, double beta, 
                                                        int niter );
    extern void f_to_accel( int internaloffset, void *array, int length );
    extern void f_from_accel( int internaloffset, void *array, int length );
    extern void r_to_accel( int internaloffset, void *array, int length );
    extern int itstep_accel_test( int my_number, int n_of_nodes, int mx, int my, 
                                                             int stw );
    
    int itstep_test( int my_number, int n_of_nodes, int mx, int my, int stw )
{
    if ( stw > (PORTION/3) ) // See comments below
     {
      fprintf( stderr, 
"stencil width(%d) greater than 1/3 of the portion size (%d), this is not allowed\n", 
                              stw, PORTION );
      return 1;
     }
    else
     return itstep_accel_test( my_number, n_of_nodes, mx, my, stw );
}
    void itstep( int mx, int my, void *pf, void *pr, double rdx2, double rdy2, 
                                 double beta, int niter )
{
    DIM2( double, f, my ) = (typeof(f))pf;
    DIM2( double, r, my ) = (typeof(r))pr;
    int k, total, is, notlast, left;
/***/
    total = mx;
// "is" is the number of rows at the beginning not to save back (i. e. from accel to CPU):
    is = 0;
    while ( total > 0 )
     {
// "k" is the number of rows to save back:
      k = total;
// "notlast" is the number of rows at the end not to save back:
      notlast = 0;
      if ( k > PORTION )
       {
        k = PORTION;
        left = total - k;
// In the line below, we silently assume that PORTION is greater than niter:
        if ( left < niter ) k -= (niter-left);
        notlast = niter;
       }
// Copy all the portion of "r" to accel:       
      r_to_accel( 0, r-is, (k+is+notlast)*my );
// Copy portion of "f" without the starting "is" rows to accel:
      f_to_accel( is*my, f, (k+notlast)*my );  
// For the very first portion, "is" is surely zero, so processing 
// immediately after the copy is OK:
      itstep_accel( k+is+notlast, my, rdx2, rdy2, beta, niter );
// In case the portion is not last, we need to copy "is" starting rows 
// of the next portion to accel, to let them be "old", not "new":
      if ( notlast )
       {
// But, in case the portion is the very first, we cannot copy "is" starting rows 
// of the next portion, because the place for it 
// in accel is occupied by useful data not saved back to CPU yet, so 
// handle the case of "is == 0" separately:
        if ( is )
         {
// Do it the usual way:
          f_to_accel( 0, f+(k-is), is*my );
          f_from_accel( is*my, f, k*my ); 
         }
        else
         {
// Perform a trick:
          f_from_accel( 0, f, niter*my );
          f_to_accel( 0, f+(k-niter), niter*my );
          f_from_accel( niter*my, f+niter, (k-niter)*my );
// The trick above is OK iff niter is NOT GREATER THAN a half of the PORTION!!!
         } 
       }
      else 
       {
// In case the portion is last, just save back:
        f_from_accel( is*my, f, k*my ); 
       }
      is = niter;
      f += k;
      r += k;
      total -= k;
     }
}

Программа совсем невелика по размеру, но простой и очевидной ее назвать трудно. Обсуждать написанное будем очень подробно.

Константа PORTION имеет отношение к объему внутренней памяти "сопроцессора" (далее будем говорить о сопроцессоре без кавычек, хотя пока это лишь его программная модель). Ее значение - это максимальное число строк обрабатываемых массивов в мелком блоке, не считая области перекрытия (теневых граней).

Далее приводятся прототипы пяти функций, которые эта функция вызывает:

    extern void itstep_accel( int mx, int my, double rdx2, double rdy2, double beta, 
                                                        int niter );
    extern void f_to_accel( int internaloffset, void *array, int length );
    extern void f_from_accel( int internaloffset, void *array, int length );
    extern void r_to_accel( int internaloffset, void *array, int length );
    extern int itstep_accel_test( int my_number, int n_of_nodes, int mx, int my, 
                                                             int stw );

Функция itstep_accel() выполняет ту же роль, что в предыдущих программах выполняла функция itstep(), то есть обрабатывает мелкий блок данных. У нее нет аргументов - обрабатываемых массивов, поскольку обрабатываемые массивы предварительно копируются в сопроцессор, запуск которого как раз и осуществляет обращение к itstep_accel(). Функции f_to_accel(), f_from_accel(), r_to_accel() выполняют копирование данных в массив f внутренней памяти сопроцессора, из этого массива, и в массив r, соответственно. При копировании данных при помощи этих функций, на стороне сопроцессора указывается смещение от начала массива, начиная с которого копировать. Смещение задается в вещественных числах типа double. Наконец, itstep_accel_test() выполняет в начале работы однократные проверки корректности обращения к сопроцессору, подобно тому, как мы это видели выше, в Разделе 5. Мнемоника "accel" в именах функций означает "accelerator" (ускоритель).

Функция itstep_accel_test(), которую мы подробно разбирали выше, в Разделе 5, выполняет новую для нас проверку, не превосходит ли коэффициент искусственного расширения области перекрытия треть максмального числа строк в мелком блоке без учета перекрытия. О том, почему при нарушении этого условия программа работать не сможет, мы узнаем, разобрав ее до конца. Если это условие выполнено, выполняются проверки следующего уровня - itstep_accel_test().

Рассмотрим работу функции itstep().

Общая логика нарезки локальной порции на мелкие блоки такова. Нам надо нарезать всю локальную порцию "встык", без перекрытий, на блоки, не превосходящие PORTION строк каждый. Каждый такой блок надо обрамить расширенными областями перекрытия, по niter строк каждая, за исключением случаев, когда край блока прилегает к краю локальной порции (у нулевого по счету блока нет области перекрытия в начале, у последнего по счету - в конце, следовательно, у единственного блока областей перекрытия нет вообще).

В программе, которую мы рассматриваем, переменная k хранит размер текущего блока без учета перекрытий (то есть блоки длиной k, вычисляемые на разных оборотах цикла нарезки, идут друг за другом "встык"), переменная is хранит размер области перекрытия перед началом блока, переменная notlast - размер области перекрытия после конца блока. В процессе нарезки есть не вполне очевидная тонкость. Пытаясь обрамить предпоследний по счету блок с конца областью перекрытия размером niter строк, мы можем обнаружить, что столько строк в конце локальной порции уже просто не осталось. Иными словами, что последний по счету блок так мал, что область перекрытия предпоследнего блока в нем не умещается, "перехлестывает" через него за границу локальной порции. Это ломает всю логику, поскольку, для корректности многократной обработки блока, необходимо выбросить в конце обработанного блока область "мусора" размером именно niter строк. Чтобы избежать этой неприятности, нам приходится при нарезке блоков заглядывать вперед, и, если надо, заранее уменьшать размер предпоследнего блока. Отрезанные при этом строки предпоследнего блока переходят в последний блок, он становится больше, и область перекрытия предпоследнего блока начинает в нем умещаться. Это делается в следующих строках программы:

     while ( total > 0 )
     {
// "k" is the number of rows to save back:
      k = total;
// "notlast" is the number of rows at the end not to save back:
      notlast = 0;
      if ( k > PORTION )
       {
        k = PORTION;
        left = total - k;
// In the line below, we silently assume that PORTION is greater than niter:
        if ( left < niter ) k -= (niter-left);
        notlast = niter;
       }

Смысл приведенного фрагмента таков:

В цикле до исчерпания локальной порции, total строк которой еще не обработаны, нарезаем "встык" блоки размером k, не превосходящим PORTION. Если на следующий раз осталось меньше niter строк, уменьшаем k на столько, чтобы строк на следующий раз осталось niter. Значение k при этом равно PORTION, значение left может быть равно минимум единице. Вычитая из PORTION значение niter-1, мы должны получить положительное число (число строк в предпоследнем блоке без учета перекрытия). Значит, значение PORTION должно быть строго больше, чем niter-1. Запомним это для написания проверок, и пойдем дальше.

Осталось самое сложное - осмыслить логику копирования данных между программой и сопроцессором. Проблема касается только копирования массива f, поскольку массив r не изменяется по ходу счета, что делает логику его копирования из программы в сопроцессор тривиальной. С массивом f все сложнее - нам надо в полной мере учесть и разрешить все трудности, с разбора которых мы начали этот раздел.

Рассмотрим момент, когда блок номер i (не последний) скопирован в сопроцессор и уже обработан в нем, то есть момент сразу после i-го по счету вызова функции itstep_accel(). В этот момент в массиве f в программе еще находятся старые значения из рассматриваемого блока, а в памяти сопроцессора - уже новые.

Точнее, в памяти сопроцессора для массива f находятся, именно в таком порядке, как указано ниже:

is         строк области перекрытия с предыдущим, i-1-м, блоком;
k          строк, составляющих тело i-го блока (новые значения);
notlast    строк области перекрытия со следующим, i+1-м, блоком.

Если наш блок является первым по счету, то is равно нулю, в противном случае is равно niter. Значение notlast точно равно не нулю, а niter, поскольку мы договорились, что блок не последний.

В первых is и последних notlast строках за счет многократной обработки блока образовался "мусор", который не подлежит обратному копированию в память программы. Средние k строк, напротив, следует скопировать из сопроцессора в память программы, и нам как раз сейчас предстоит это сделать. Однако, последние niter из этих средних k строк потребуются следующему, i+1-му, блоку в качестве области перекрытия с нашим блоком, который тогда уже станет предыдущим. При этом потребуются ему, согласно итерационному процессу Якоби, старые значения этих строк. Они пока еще есть в массиве f в программе, но, как только мы скопируем наши k строк из сопроцессора в программу, мы их "затрем". Значит, последние niter из средних k строк текущего блока массива f в программе надо куда-то деть, прежде чем копировать из сопроцессора в программу новые значения средних k строк на их законное место.

Проще всего скопировать их "впрок" в сопроцессор. На том месте, куда эти niter строк должны попасть для обработки следующего блока, в сопроцессоре сейчас все равно лежит is строк "мусора", поверх которого грех не скопировать данные для следующего цикла обработки. Однако, это можно делать только в случае, когда is не равно нулю, то есть когда необходимое нам место действительно занято "мусором". Напомним, что is равно нулю только в случае обработки самого первого блока. Отложив пока в сторону случай, когда is равно нулю, рассмотрим часть программы, в которой реализуется описанная выше логика.

// In case the portion is not last, we need to copy "is" starting rows 
// of the next portion to accel, to let them be "old", not "new":
      if ( notlast )
       {
// But, in case the portion is the very first, we cannot copy "is" starting rows 
// of the next portion, because the place for it 
// in accel is occupied by useful data not saved back to CPU yet, so 
// handle the case of "is == 0" separately:
        if ( is )
         {
// Do it the usual way:
          f_to_accel( 0, f+(k-is), is*my );
          f_from_accel( is*my, f, k*my ); 
         }
        else
         {
// Perform a trick:
          f_from_accel( 0, f, niter*my );
          f_to_accel( 0, f+(k-niter), niter*my );
          f_from_accel( niter*my, f+niter, (k-niter)*my );
// The trick above is OK iff niter is NOT GREATER THAN a half of the PORTION!!!
         } 
       }
      else 
       {
// In case the portion is last, just save back:
        f_from_accel( is*my, f, k*my ); 
       }

Внешний "if" приведенного фрагмента программы выясняет, не является ли обрабатываемый блок последним по счету, и, если является, просто копирует из сопроцессора k строк в память программы на их законное место (ветка "else" в конце приведенного фрагмента текста программы).

Далее проверяется, не равно ли нулю is. Если не равно, niter последних строк тела текущего блока копируются из памяти программы в сопроцессор, и только затем все тело (k строк) текущего блока (новые значения) копируется из сопроцессора в память программы.

Наконец, если is равно нулю (блок первый по счету, но не последний), приходится действовать специальным образом. В приведенном фрагменте текста программы это место отмечено комментарием "Perform a trick" - "Делаем фокус".

Сложность этой ситуации в том, что в первых niter строках блока в сопроцессоре, вместо "мусора", поверх которого не жалко скопировать область перекрытия для следующего блока, находятся полезные данные. Значит, прежде чем действовать обычным образом, как для случая is, не равного нулю, эти полезные данные из начала блока в сопроцессоре ("голову" тела блока) надо куда-то деть. Проще всего скопировать их (и только их!) из сопроцессора в массив f в памяти программы, на свое законное место, а затем действовать обычным образом: niter последних строк тела блока скопировать из памяти программы в сопроцессор, после чего докопировать из сопроцессора в память программы оставшийся "хвост" тела блока. Игра здесь идет на том, что старые значения, необходимые для области перекрытия в начале следующего блока, расположены в текущем блоке именно в "хвосте", а не в "голове".

Чтобы описанный здесь "фокус" удался, между размерами "хвоста" и "головы" должны соблюдаться некоторые пропорции. Попросту говоря, они не должны пересекаться. Размер "головы", как мы знаем, равен niter строкам, как и число строк в конце "хвоста", которые мы копируем на место предварительно "спасенной" "головы". Размер тела не последнего блока равен PORTION, "голова" и "хвост" вместе должны уложиться в этот размер. Значит, условием корректной работы данной программы является требование, чтобы значение niter не превосходило половины значения PORTION. Ранее мы формулировали более слабое требование на допустимое значение niter, которое поглощается данным.

В приведенном только что рассуждении нельзя не заметить небольшой неточности. Как мы видели выше, обсуждая "нарезку" локальной порции на блоки, размер тела не последнего (точнее, предпоследнего) блока может быть и меньше PORTION. Значит, сформулированное нами требование к области допустимых значений niter все еще слишком слабое. Уточним его. Прежде всего, может ли блок быть одновременно первым и предпоследним? Безусловно, может, если локальная порция состоит всего из двух блоков. Вспомним, на сколько строк максимум может быть "урезан" предпоследний блок? Очевидно, не более чем на niter строк. Значит, правильное ограничение на niter - не более трети значения PORTION.

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

#include <stdio.h>
#include <string.h>
#include "dimension.h"
#define LOCALMEM 1000000
#define MAXMY 10000
// The arrays below mimics the accelerator's local memory.
    static double localf[LOCALMEM];
    static double localr[LOCALMEM];
    void f_to_accel( int internaloffset, void *array, int length )
{
    memcpy( localf+internaloffset, array, length*sizeof(localf[0]) );
}
    void f_from_accel( int internaloffset, void *array, int length )
{
    memcpy( array, localf+internaloffset, length*sizeof(localf[0]) );
}
    void r_to_accel( int internaloffset, void *array, int length )
{
    memcpy( localr+internaloffset, array, length*sizeof(localr[0]) );
}
    int itstep_accel_test( int my_number, int n_of_nodes, int mx, int my, int stw )
{
    if ( my > MAXMY )
     {
      fprintf( stderr, " my (%d) greater than maxmy (%d)\n", my, MAXMY );
      return 1; 
     }
    else return 0;
}
    void itstep_accel( int mx, int my, double rdx2, double rdy2, double beta, int niter )
{
    DIM2( double, f, my ) = (typeof(f))localf;
    DIM2( double, r, my ) = (typeof(r))localr;
    static double  lower[MAXMY];
    static double middle[MAXMY];
    int i, j, k, mx1, my1;
/***/
    mx1 = mx - 1;
    my1 = my - 1;
    for ( k = 0; k < niter; k++ )
     {
      for ( i = 0; i < mx1; i++ )
       {
        for ( j = 0; j < my; j++ )
         {
          lower[j] = middle[j];
          middle[j] = f[i][j];
         } 
        if ( i )
         {
          for ( j = 1; j < my1; j++ )
           {
             f[i][j] = ((lower[j]+f[i+1][j])*rdx2+(middle[j-1]+middle[j+1])*rdy2-r[i][j])
                                        *beta;
           }
         }
       }
     }
}

Тексты эти, как мы видим, достаточно просты, и в специальных пояснениях не нуждаются. Обратим внимание на то, что проверка корректности использования этих функций в itstep_accel_test() не полна. Отсутствует проверка того, что в локальной памяти сопроцессора способны уместиться блоки массивов f и r размером PORTION+2*stw строк. И это совсем не удивительно: выбирая набор аргументов для itstep_accel_test(), мы просто забыли в качестве одного из аргументов передать значение PORTION. В следующих версиях программы мы эту ошибку исправим.

 
 
 
 
 
 
 
  Тел. +7(499)220-79-72; E-mail: inform@kiam.ru