OpenMP, Linux и немного фана

Попался вот такой кусок кода для вычисления тривиального числа Pi

n = 1000000000
62    h   = 1.0 / (double) n;
63    sum = 0.0;
64
65    for (i = 1; i <= n; i++)
66    {
67           x = h * ((double)i - 0.5);
68           sum +=(4.0 / (1.0 + x*x));
69    }
70
71    mypi = h * sum;

изначально все строилось для MPI , правда захотелось перенести на OpenMP

казалось бы все просто:

 53 #pragma omp parallel for reduction(+: sum)
 54 for (int i = 1; i <= n; i++)
 55    {
 56
 57       x = h * ((double)i - 0.5);
 58       sum +=(4.0 / (1.0 + x*x));
 59    }

(гусары - молчать , пока что ... ).

(компилируем как gcc ./pim.c -o ./pim -O3 -fopenmp -mfpmath=sse)

Практически все (окромя точности, но об этом позже) устраивает и почти красиво, пример Thread Profiler for Linux:

tprofile_cl ./pim
Building project
Instrumenting
 16% pim             ( API Imports ):..

Running:  /mpi/pi/pi/pim

n = 1000000000
Num threads before = 1
myid: 1 i: 0: 0 c:0
mypi: 3.1415926535898211

Application finished

Intel® Thread Profiler  3.1  Summary Report
application:          /mpi/pi/pi/pim
collection:           Wed Jan  7 17:52:15 2009
runtime:              3.43215s
# of processors:      4
# of threads:         4
# of waits:           3
wait frequency:       0.874088
average concurrency:  3.96682
Concurrency:
   0 [..........] 0%      0
   1 [..........] 0.913%  0.0301119
   2 [..........] 0.0051% 0.000168224
   3 [..........] 0.349%  0.0114959
   4 [##########] 98.7%   3.25637

Паралельность просто фантастическая!

однако сюрпризы только начинаються если присмотреться подольше и поближе :) ,

вот пример нескольких последовательных, и не только, запусков:

[maxd@fc9i pi]$ time ./pim
n = 1000000000
mypi:   3.1415926535898211

real    0m3.345s
user    0m13.125s
sys     0m0.032s
[maxd@fc9i pi]$ time ./pim
n = 1000000000
mypi:   3.1415926535898211
real 0m55.433s
user 3m38.935s
sys 0m0.033s

В итоге - иногда время работы ~3 сек и очень часто - ~55 сек - крута, 18х slowdown. Что то мы не подкрутили!

посмотрим поглубже и видим что то типа

  400710:       f2 0f 2a c0             cvtsi2sd %eax,%xmm0
  400714:       83 c0 01                add    $0x1,%eax
  400717:       66 0f 28 ea             movapd %xmm2,%xmm5
  40071b:       39 c2                   cmp    %eax,%edx
  40071d:       f2 0f 5c c4             subsd  %xmm4,%xmm0
  400721:       f2 41 0f 59 44 24 08    mulsd  0x8(%r12),%xmm0
  400728:       f2 41 0f 11 44 24 10    movsd  %xmm0,0x10(%r12)
  40072f:       f2 0f 59 c0             mulsd  %xmm0,%xmm0
  400733:       f2 0f 58 c3             addsd  %xmm3,%xmm0
  400737:       f2 0f 5e e8             divsd  %xmm0,%xmm5
  40073b:       f2 0f 58 cd             addsd  %xmm5,%xmm1
  40073f:       7d cf                   jge    400710 <main.omp_fn.0+0x70>

(btw AT&T синтаксис вовсе не дурен)

запустив VTune/PTU увидим впечатляющий CPI (циклов на инструкцию) ~ 45 :) при желаемом -как можно ближе к 0.25 но меньше 1 уже хорошо . Кстати - саму ошибку уже видно с разных сторон. Для разнообразия запустим Thread Checker for Linux и уже наглядно видим что то типа

 tcheck_cl ./pim
Intel® Thread Checker 3.1 command line instrumentation driver (26185)
Copyright (c) 2007 Intel Corporation. All rights reserved.
Building project
Instrumenting
 16% pim             ( All Functions ):..

Running:  /mpi/pi/pi/pim

n = 1000000000
Num threads before = 1
myid:   1       i:      0: 0    c:0
mypi:   3.1415926535898211

Application finished

_______________________________________________________________________________
|ID|Short De|Seve|C|Contex|Description                        |1st Acc|2nd Acc|
|  |scriptio|rity|o|t[Best|                                   |ess[Bes|ess[Bes|
|  |n       |Name|u|]     |                                   |t]     |t]     |
|  |        |    |n|      |                                   |       |       |
|  |        |    |t|      |                                   |       |       |
_______________________________________________________________________________
....


|3 |Write ->|Erro|2|[pim, |Memory write at [pim, 0x728]       |[pim,  |[pim,  |
|  |Write da|r   |1|0x6a0]|conflicts with a prior memory write|0x728] |0x728] |
|  |ta-race |    |4|      |at [pim, 0x728] (output dependence)|       |       |
|  |        |    |8|      |                                   |       |       |
|  |        |    |7|      |                                   |       |       |
|  |        |    |6|      |                                   |       |       |
|  |        |    |1|      |                                   |       |       |
|  |        |    |1|      |                                   |       |       |

 ....

 

(форматирование несколько "поехало") - кстати Thread Checker оказался самым быстрым способом борьбы.

Тут даже придумано красивое название - False Sharing: несколько потоков мешают друг другу работая с одним ресурсом. Как другой вариан объяснения - читайте blog недавно были интересные блоги насчет красивых и иногда жизненых примеров для общих проблем.

На четырех ядерной машинке 18х ухудшения - интересно.

Как вылечить ? Например так :

 55 #pragma omp parallel for private(x) reduction(+: sum)
 56 for (i = 1; i <= n; i++)
 57    {
 58       x = h * ((double)i - 0.5);
 59       sum +=(4.0 / (1.0 + x*x));
 60    }

добавленое - выделено.

Кстати - исходный код был С где declaration для переменной x было вне цикла, как можно было сделать для С++ и тогда проблемы не было бы вообще.

в итоге наш цикл превращаеться в более простую конструкцию типа:

  400710:       f2 0f 2a c0             cvtsi2sd %eax,%xmm0
  400714:       83 c0 01                add    $0x1,%eax
  400717:       66 0f 28 f3             movapd %xmm3,%xmm6
  40071b:       39 c2                   cmp    %eax,%edx
  40071d:       f2 0f 5c c5             subsd  %xmm5,%xmm0
  400721:       f2 0f 59 c2             mulsd  %xmm2,%xmm0
  400725:       f2 0f 59 c0             mulsd  %xmm0,%xmm0
  400729:       f2 0f 58 c4             addsd  %xmm4,%xmm0
  40072d:       f2 0f 5e f0             divsd  %xmm0,%xmm6
  400731:       f2 0f 58 ce             addsd  %xmm6,%xmm1
  400735:       7d d9                   jge    400710 <main.omp_fn.0+0x70>

которая выполняеться и быстрее и стабильнее.

Насчет точности:

Используемый алгоритм подвержен проблеме округления в частности при использовании OpenMP.

например для использования SSE можно попробывать использовать

gcc ./pim.c -o ./pim -O3 -fopenmp -mfpmath=sse

для FPU:

gcc ./pim.c -o ./pim -O3 -fopenmp -mfpmath=387

результаты будут разные:

sse: 3.1415926535898211
387: 3.1415926535897936
org: 3.14159265358979323846
где org - просто Pi взятое из какой либо википедии, кстати FPU более точнее но об этом я уже как то писал.

Время работы FPU кстати неплохое

sse: real 0m3.280s

fpu: real    0m3.894s

 

Присутствующая векторизация в gcc для этого цикла не сработала и пришлось бы менять код.

Intel Compiler сделал вещи шустрее (используя packed arguments instead of scalar ~ векторизация , loop unrolling etc )  и с уже приведенным C кодом :

icc ./pim.c -O3 -xSSSE3 -openmp -vec-report5 -o ./pim.icc

и результат неплох

real    0m1.857s
user    0m6.528s
sys     0m0.033s

За сим прощаюсь с успехом поздравив Вас в удачной оптимизации :)

Nähere Informationen zur Compiler-Optimierung finden Sie in unserem Optimierungshinweis.