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

Модельное приложение выберем то же самое, что в «Десяти простых шагах». Приведем здесь несколько модифицированную выдержку из Шага 3:

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dimension.h"
// Change here the number of steps, the cell geometry, etc
#define NITER 5000
#define STEPITER 1000
#define delx 0.5
#define dely 0.25
// end change here.
       int main( int argc, char **argv )
{
       int i, j, n, mx1, my1, mx, my;
       FILE *fp;
       double rdx2, rdy2, beta;
/***/
       if ( argc != 3 )
        {
         fprintf( stderr, "Usage: %s <nrows> <ncolumns>\n", argv[0] );
         return (-1);
        }
       mx = (int)atol( argv[1] );
       my = (int)atol( argv[2] );
       mx += 2;
       my += 2;
/* Here we know the array sizes, so make the arrays themselves: */
       {
        DIM2( double, f, my );
        DIM2( double, newf, my );
        DIM2( double, r, my );
        void *tmpf;
/***/
        f = (typeof(f))malloc( mx*sizeof(*f) );
          newf = (typeof(newf))malloc( mx*sizeof(*f) );
        r = (typeof(r))malloc( mx*sizeof(*r) );
        if ( (!f) || (!newf) || (!r) )
         {
          fprintf( stderr, "Cannot allocate, exiting\n" );
          return( -1 );
         }
        rdx2 = 1./delx/delx;
        rdy2 = 1./dely/dely;
        beta = 1.0/(2.0*(rdx2+rdy2)); 
        printf( "Solving task on %d by %d grid\n", mx-2, my-2 );     
        fflush( stdout );
        for ( i = 0; i < mx; i++ )
         {
          for ( j = 0; j < my; j++ )
           {
            if ( (i==0) || (j==0) || (i==(mx-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;        
           }
         }    
        mx1 = mx - 1;
        my1 = my - 1; 
/* Iteration loop: */
        for ( n = 0; n < NITER; n++ )
         {
           if ( !(n%STEPITER) ) printf( "Iteration %d\n", n );
/* Step of calculation starts here: */
           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;
              }
            }
/* swap the halves: */
	      tmpf = newf;      
          newf = f;
	     f = (typeof(f))tmpf;
         }    
        fp = fopen("output_seq.dat", "w" );
        for ( i = 1; i < (mx-1); i++ )
         {
          fwrite( &(f[i][1]), my-2, sizeof(f[0][0]), fp );
         } 
       fclose( fp );    
       }     
       return 0;
}

Несколько кратких пояснений.

Температура нормализована к диапазону от 0.0 до 1.0. Уравнение решается с правой частью, представленной массивом r, который в данном расчете заполняется нулями, но в пересчете значений температуры, тем не менее, участвует.

Имеются два указателя, f и newf, которые указывают на два разных массива значений температуры, меняясь ролями на каждой итерации (если на текущей итерации f указывает на первый массив значений температуры, а newf – на второй, то на следующей итерации будет наоборот). При этом тот массив значений температуры, на которую на данной итерации указывает f, играет роль «старых» значений, а тот, на которую указывает newf, роль «новых».

В остальном программа вполне тривиальна. Отныне она является для нас руководством к действию. Естественно, мы будем стремиться к параллельной, а затем — и гибридно-параллельной, реализации в точности тех вычислений, которые представлены приведенной выше последовательной программой.

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