Ошибка производительности в методе конечных разностей

Question

Я написал фрагмент кода C, который использует метод конечных разностей для оценки значений. Это метод усреднения. Я профилировал код и обнаружил, что одна функция iterate() является самой медленной.

void iterate(double data[][ARRAY_SIZE], int nx, int ny, int dx, int dy) 
{ 
    for (int i = 0; i < nx; ++i) 
    { 
     for (int j = 0; j < ny; ++j) 
     { 
      if (i % (dx + 1) == 0 && j % (dy + 1) == 0) 
       continue; 
      else if (i == 0 && 0 < j && j < ny) 
       data[i][j] = (data[i][j - 1] + data[i][j + 1] + data[i + 1][j])/3; 
      else if (j == 0 && 0 < i && i < nx) 
       data[i][j] = (data[i - 1][j] + data[i + 1][j] + data[i][j + 1])/3; 
      else if (i == nx - 1 && 0 < j && j < ny) 
       data[i][j] = (data[i][j - 1] + data[i][j + 1] + data[i - 1][j])/3; 
      else if (j == ny - 1 && 0 < i && i < nx) 
       data[i][j] = (data[i - 1][j] + data[i + 1][j] + data[i][j - 1])/3; 
      else 
       data[i][j] = (data[i - 1][j] + data[i + 1][j] + data[i][j - 1] + data[i][j + 1])/4; 
     } 
    } 
} 

Этот цикл работает медленно, и я не уверен, что здесь отсутствует, что замедляет работу. Есть ли лучший способ сделать то же самое?

2000 итераций с массивом 400x400double занимает

real 0m1.950s 
user 0m1.940s 
sys 0m0.004s 

Answer 1

Вот некоторые идеи:

1. Оно кажется, что ny должна равняться ARRAY_SIZE. Вы можете также опустить его как параметр и просто использовать константу времени компиляции.
2. Все предложения if/else, кроме последнего, применимы только к определенной строке или столбцу. Так вытащите их. Например, вы можете обрабатывать первую строку и столбец как 1D-циклы, прежде чем делать всю матрицу за пределами краев, а затем обработать крайний правый столбец и нижнюю строку.

В конце концов, ваш основной цикл должен быть больше, как это:

for (int i = 1; i < nx - 1; ++i) 
{ 
    for (int j = 1; j < ARRAY_SIZE - 1; ++j) 
    { 
     data[i][j] = (data[i - 1][j] + data[i + 1][j] + data[i][j - 1] + data[i][j + 1])/4; 
    } 
} 

Answer 2

Рассмотрим реализацию:

void iterate(double data[][ARRAY_SIZE], int nx, int ny, int dx, int dy) 
{ 
    // because nx - 1 and ny - 1 are used 
    nx--; 
    ny--; 
    // because dx + 1 and dy + 1 are used 
    dx++; 
    dy++; 

    int i = 0; 
    int j = 0; 

    // case i == 0 && 0 < j && j < ny 
    for (j = 1; j < ny; ++j) 
    { 
     if (j % dy) 
      data[0][j] = (data[i][j - 1] + data[i][j + 1] + data[i + 1][j])/3.0; 
    } 

    j = 0; 

    // case j == 0 && 0 < i && i < nx 
    for (i = 1; i < nx; ++i) 
    { 
     if (i % dx) 
      data[i][0] = (data[i - 1][j] + data[i + 1][j] + data[i][j + 1])/3.0; 
    } 

    // default case 
    for (i = 1; i < nx; ++i) 
    { 
     for (j = 1; j < ny; ++j) 
     { 
      if (i % dx || j % dy) 
       data[i][j] = (data[i - 1][j] + data[i + 1][j] + data[i][j - 1] + data[i][j + 1]) * 0.25; 
     } 
    } 

    // case i == nx && 0 < j && j < ny 
    for (j = 1; j < ny; ++j) 
    { 
     if (nx % dx || j % dy) 
      data[nx][j] = (data[i][j - 1] + data[i][j + 1] + data[i - 1][j])/3.0; 
    } 

    // case j == ny && 0 < i && i < nx 
    for (i = 1; i < nx; ++i) 
    { 
     if (ny % dy || i % dx) 
      data[i][ny] = (data[i - 1][j] + data[i + 1][j] + data[i][j - 1])/3.0; 
    } 
} 

Основные три точки являются:

1. уменьшить количество операций во внутреннем контуре для двойного цикла for
2. уменьшить количество тривиальных операций, делая их только один раз
3. не следует смешивать типы данных и силы принуждения (использование / 3.0 и * 0.25)

Единственное в своем коде не объяснил, что i % dx || j % dy равно до !(i % dx == 0 && j % dy == 0).