Что касается картинок. Я написал программу и запускал ее в учебном классе на 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-1802. Tomas Petricek. Accelerator and F# (II.): The Game of Life on GPU
3. Последняя версия исходного кода к статье на bitbucket
