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

Обратим внимание на то, что, при разбиении локальной порции на мелкие блоки, каждый из которых умещается в сопроцессоре, блоки эти придется нарезать с перекрытиями (теневыми гранями), подобно тому, как нарезается на локальные порции распределенный массив при распараллеливании приложения. Очевидно, что размер перекрытия при нарезке мелких блоков такой же, как и при нарезке локальных порций в параллельной реализации (в данном случае - одна строка). А что будет, если нарезать мелкие блоки с перекрытием в две строки? Каждый такой мелкий блок можно будет, единожды скопировав в сопроцессор, обработать два раза. Правда, при обратном копировании обработанных данных в память, придется "выбросить" области перекрытия удвоенного размера, но, по сравнению с двойной экономией на копировании данных, это не страшно. Легко показать, что при увеличении размера перекрытия в N раз по сравнению с требованиями исходного численного метода, единожды скопированный в сопроцессор мелкий блок может быть обработан внутри сопроцессора N раз. Значит, следующий шаг преобразования нашей программы - это включение в нее нарезки локальной порции на мелкие блоки с произвольной шириной перекрытия? Нет, к этому преобразованию переходить еще рано. Для того, чтобы многократная обработка мелких блоков данных была корректной, необходимо, чтобы была корректной многократная (без обменов краями с соседями) обработка всей локальной порции. Для этого, в свою очередь, потребуется нарезать уже сами локальные порции с расширенными в N раз перекрытиями, а обмены такими расширенными краями сделать в N раз более редкими. С этого преобразования программы и начнем. Текст преобразованной таким образом программы приводится ниже.

Текст функции main():

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "vsda.h"
#include "dimension.h"
// Change here the number of steps, the cell geometry, etc
#define STEPITER 1000
#define delx 0.5
#define dely 0.25
#define NITER 10002
#define STW 3
// end change here.
/***/
       extern void itstep( int mx, int my, void *f, void *newf, void *r, 
                                           double rdx2, double rdy2, double beta, int niters );
       int main( int argc, char **argv )
{
       int i, j, k, n, mx1, mx2, my1, my_number, n_of_nodes, 
                                                   totalmx, partmx, leftmx, mx, my, niter, stw;
       FILE *fp;
       double howlong;
       double rdx2, rdy2, beta;
       VSDA vsda;
       long *distr_out;
       MPI_Request rq[5];
/***/
       MPI_Init( &argc, &argv );
       MPI_Comm_rank( MPI_COMM_WORLD, &my_number );
       MPI_Comm_size( MPI_COMM_WORLD, &n_of_nodes );
       if ( argc != 3 )
        {
           if ( !my_number )
            {
              fprintf( stderr, "Usage: %s <nrows> <ncolumns>\n", argv[0] );
              fflush( stderr );
            }
           return (-1);
        }
       totalmx = mx = (int)atol( argv[1] );
       my = (int)atol( argv[2] );
       if ( my < 1 )
        {
          if ( !my_number )
           {
             fprintf( stderr, "Number of columns (%d) should be positive\n", my );
             fflush( stderr );
           }
          return (-1); 
        }
       distr_out = (long*)malloc( n_of_nodes*sizeof( *distr_out ) );
       if ( !distr_out )
        {
         fprintf( stderr, "No memory\n" );
         fflush( stderr );
         return( -1 );
        }
       totalmx += 2;    
       my += 2;
       stw = STW;
       niter = NITER;
/* Make STW odd, NITER divisible by STW: */
       if ( !(stw%2) ) stw--;
       niter -= (niter%stw);
/* Here we know the array width, so make the arrays themselves: */
       {
        DIM2( double, f, my );
        DIM2( double, newf, my );
        DIM2( double, r, my );
        void *tmpf;
/* Make the distributed array, discover the local portion length: */
        VSDA_Create( MPI_COMM_WORLD, VSDA_BLOCKED, 0, 
                        (long)totalmx, (long)sizeof(*f), (long)stw, NULL, distr_out, &vsda );
        mx = distr_out[my_number];
        if ( n_of_nodes > 1 )
         {
          if ( my_number > 0 ) mx += stw;
          if ( my_number < (n_of_nodes-1) ) mx += stw;
         }
        f = (typeof(f))malloc( mx*sizeof(*f) );
        newf = (typeof(newf))malloc( mx*sizeof(*newf) );
        r = (typeof(r))malloc( mx*sizeof(*r) );
        if ( (!f) || (!newf) || (!r) )
         {
          fprintf( stderr, "Cannot allocate, exiting\n" );
          exit( -1 );
         }
        rdx2 = 1./delx/delx;
        rdy2 = 1./dely/dely;
        beta = 1.0/(2.0*(rdx2+rdy2)); 
        if ( !my_number )
         {
           printf( "Solving heat conduction task on %d by %d grid by %d processors\n",
                  totalmx-2, my-2, n_of_nodes );     
           fflush( stdout );
         }
        for ( i = 0; i < mx; i++ )
         {
          for ( j = 0; j < my; j++ )
           {
            if ( ((i==0)&&(my_number==0)) || (j==0) || 
                  ((i==(mx-1))&&(my_number==(n_of_nodes-1))) || (j==(my-1)) )
             newf[i][j] = f[i][j] = 1.0;
            else
             newf[i][j] = f[i][j] = 0.0;
            r[i][j] = 0.0;
           }
         }
/* Iteration loop: */
        howlong = MPI_Wtime();
        for ( n = 0; n < niter/stw; n++ )
         {
          if ( !my_number )
           {
             if ( !(n%STEPITER) )
             printf( "Iteration %d\n", n );
           }
/* Step of calculation starts here: */
          itstep( mx, my, f, newf, r, rdx2, rdy2, beta, stw );
/* Do all the transfers: */
          VSDA_Update_begin( &vsda, rq, newf );
          VSDA_Update_end( rq );
/* swap the halves: */      
          tmpf = newf;
          newf = f;
          f = (typeof(f))tmpf;
         }
        if ( !my_number )
         {
          printf( "Elapsed time: %f sec\n", (float)(MPI_Wtime()-howlong) );
          fp = fopen( "output_par_2.dat", "w" );
          fclose( fp );
         }
        for ( j = 0; j < n_of_nodes; j++ )
         {
          MPI_Barrier( MPI_COMM_WORLD );
          if ( j == my_number )
           {
            int istart, isize;
            istart = 1;
            isize = mx-1; 
            if ( n_of_nodes > 1 )
             {
              if ( my_number > 0 ) istart = stw;
              if ( my_number < (n_of_nodes-1) ) isize = mx-stw;
             }
            fp = fopen( "output_par_2.dat", "a" );
            for ( i = istart; i < isize; i++ )
             {
              fwrite( &(f[i][1]), my-2, sizeof(f[0][0]), fp );
             } 
            fclose( fp );
           }
         }
       }
       MPI_Finalize();
       return 0;
}

Обратим внимание на то, что не только усложнять, но и вообще как-либо изменять ту часть логики программы, которая организует обмены данными между MPI - процессами, нам не пришлось, что свидетельствует о действительной пригодности библиотеки VSDA для упрощения разработки задач этого класса.

Текст функции itstep():

#include "dimension.h"
      void itstep( int mx, int my, void *pf, void *pnewf,
                                 void *pr, double rdx2, double rdy2, double beta, int niter )
{
      DIM2( double, f, my ) = (typeof(f))pf;
      DIM2( double, newf, my ) = (typeof(newf))pnewf;
      DIM2( double, r, my ) = (typeof(r))pr;
      int i, j, k, mx1, my1;
      void *tmpf;
/***/
      mx1 = mx - 1;
      my1 = my - 1;
      for ( k = 0; k < niter; k++ )
       {
        for ( i = 1; i < mx1; i++ )
         {
          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;
           }
         }
        tmpf = newf;
        newf = f;
        f = (typeof(f))tmpf;
       }
}

Приведем необходимые пояснения. Коэффициент расширения теневых граней задан константой STW (STencil Width, ширина шаблона) в тексте функции main(). Чтобы избежать еще большего усложнения программы, нам пришлось наложить довольно противоестественные ограничения на исходные данные. Мы требуем, чтобы число итераций было кратно STW, и чтобы само STW было нечетным. Последнее требование связано с тем, что и в функции main(), и в вызываемой из нее функции itstep(), независимо друг от друга меняются местами указатели f и newf. Для корректной работы программы необходимо, чтобы они менялись местами согласованно, то есть чтобы после вызова itstep() результаты счета были в том массиве, на которую в вызывающей программе указывает newf. А для этого необходимо и достаточно, чтобы число замен местами f и newf, которое равно STW, было нечетным.

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