Вычисления на GPU в .NET приложениях

В 1991 году мне в руки попался журнал Scientific American [1] с красивыми фрактальными картинками. В 90-х фракталы были в моде, даже растровый редактор был назван Fractal Design Painter, а в журналах обсуждались фрактальные алгоритмы сжатия изображений с умопомрачительной эффективностью (не менее 1000:1). В итоге оказалось что алгоритмы показывают обещанные результаты только для изображений типа листа папоротника.
Что касается картинок. Я написал программу и запускал ее в учебном классе на PC 8088, а преподаватель дал погонять на 80486 (с встроенным сопроцессором) в лаборатории в институте. Сейчас трудно вспомнить результаты, приблизительно 320x200 за 5 минут, это на 8088 но с сопроцессором.
Пару лет назад, с появлением DirectX 10, GPU, CUDA, захотелось посмотреть что же дал прогресс за последние 20 лет, а тут еще попалась пара интересных статей про MS Accelerator. Так как во фракталах Ляпунова значение каждой точки вычисляется независимо, то есть вычисления идеально распараллеливаются - задача как раз для GPU.

Первые цифры

Метрику для обозначения производительности возьмем такую - "миллион точка-итераций в секунду" (mega iterations per second, mis). Это отражает особенность алгоритма - в большинстве случаев сходимость достигается после 500 итераций, на значительная часть областей требует 2-4 тыс. итераций.
Исходный код выглядит так:

protected static double CalculateExponent(double[] pattern, double initial, int warmup, int iterations)
{
 var x = initial;
 var patternSize = pattern.Length;

 for (var i = 0; i < warmup; i++)
 {
  var r = pattern[i % patternSize];
  x *= r * (1 - x);
 }

 double total = 0;
 for (var i = warmup; i < iterations; i++)
 {
  var r = pattern[i % patternSize];
  var d = Math.Log(Math.Abs(r - 2 * r * x));

  total += d;
  if (Double.IsNaN(d) || Double.IsNegativeInfinity(d) || Double.IsPositiveInfinity(d))
  {
   return d;
  }

  x *= r * (1 - x);
 }

 return total / Math.Log(2) / (iterations - warmup);
}

Многопоточная реализация для CPU без использования SSE показала 42 mis на двухядерном Core2 6400 @3.0GHz. Немного, с учетом примерно 12 операций выборки/вычисления на цикл, или ~500 млн на два ядра. Около 10 тактов на операцию, что слишком плохо. Вероятно Log() и арифметика в таком режиме плохо конвейеризируются. Это отдельная задача для исследований.

MS Accelerator

Следом была написана реализация под Accelerator v.2, очень простая библиотека в идеологии SIMD, где формула выражение строится над большой матрицей. Библиотека содержит аналоги практически всех функций System.Math и перегружает операторы. Адресация соседних ячеек возможна с помощью функций сдвига матрицы. Библиотека действительно очень проста, изучается за один вечер и позволяет писать наглядный код типа:

var dimensions = new[] { 1000, 1000 };
var x = new FPA((float)settings.InitialValue, dimensions);

// creating A,B arrays on the fly a few times faster!
var fr = new Func<int, FPA>(i => Math.Replicate(new FPA(
 settings.Pattern[i % settings.Pattern.Length] == 'a' ? aArray : bArray), w, h));

// warmup cycle, no limit calculation
for (var i = 0; i < settings.Warmup; i++)
{
 var r = fr(i);
 x *= r * (1.0f - x);
}

var total = new FPA(0, dimensions);
for (var i = settings.Warmup; i < settings.Iterations; i++)
{
 var r = fr(i);

 total += Math.Log2(Math.Abs(r - 2 * r * x));
 x *= r - r * x;
}
total *= 1f / (settings.Iterations - settings.Warmup);

return total;

Результат в среднем 250 mis, но на некоторых разрешениях получается до 400 mis, то есть быстрее почти в 10 раз.
На этом можно было бы остановиться, но от видеокарты со 192 процессорами, работающими на частоте 1 ГГЦ ожидаешь хотя бы 100-кратного преимущества. Также были замечены разные странности, частично отраженные в вышеприведенном коде, и в частности, чрезмерное использование памяти видеокарты.

В результате, проштудировав соответствующую статью на Хабре, я решил попробовать библиотеку Brahma (она же за деньги под названием Tidepowered).

Brahma

Собственно код я написал и он доступен в исходниках. Но мне не хотелось переезжать на тормозную VS 2010, а на компиляторе C# 3.5 из-под VS2008 выражения (expression) получались, как выяснилось, неразборчивыми для Brahma. После пары часов упражнений я решил что мне эта прослойка не нужна и я быстрее напишу нативный код для OpenCL.
Смысл использования Brahma/Tidepowerd мне непонятен - документации мало, денег стоит много, синтаксис странный, кроссплатформенность неважная (например, вышеупомянутые проблемы на .NET 3.5), диагностика в случае ошибок никакая. Хотя код, в общем-то, нагляден:

(ds, a, b, m, t) =>
 from pt in ds
 let i = pt.GlobalID0
 let j = pt.GlobalID1
 let x = initialX
 let r = default(float32)
 let bv = b[i]
 let av = a[j]
 let warmup = engine.Loop(0, warmupCount, idxs =>
  idxs.Select(idx =>
   new Set[]
    {
     r <= (m[idx%maskLen] == 0 ? av : bv),
     x <= r*x - r*x*x
    })
  )
 let total = default(float32)
 let calculate = engine.Loop(warmupCount, iterationsCount, idxs =>
  idxs.Select(idx =>
   new Set[]
    {
     r <= (m[idx%maskLen] == 0 ? av : bv),
     total <= total + Math.Log(Math.Fabs(r - 2*r*x)),
     x <= r*x - r*x*x
    })
  )
 select new Set[]
  {
   t[j*columns + i] <= total*divider
  }

OpenCL/Cloo

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

Ссылки

1. "Leaping into Lyapunov Space", by A. K. Dewdney, Scientific American, Sept. 1991, pp. 178-180
2. Tomas Petricek. Accelerator and F# (II.): The Game of Life on GPU
3. Последняя версия исходного кода к статье на bitbucket