Я являюсь автором проекта по математическому моделированию прикладной механики и в работе моей программы до 90% вычислительного времени уходит на решение системы линейных уравнений. Цель этой статьи сугубо практическая - найти оптимальный метод решения системы линейных уравнений с точки зрения производительность/трудозатрат для небольшого проекта и рассказать о результате.

В прошлом я уже несколько раз обращал внимание на вычисления на GPU, но всегда что-то останавливало. И вот у меня накопился достаточный практический опыт программирования на C/C++ и наконец дошли руки, чтобы протестировать OpenCL.

В результате изучения в поисковых системах и в различных списках, включая страницу Википедии, были найдены следующие проекты с линейной алгеброй:

• clMagma - его я открыл и закрыл. Насколько я понял, это проект переноса на OpenCL алгоритмов с языка C, которые в свою очередь переписаны с Fortran. Я уважаю Fortran, но посчитал нецелесообразным осваивать столь специфичный код для данного проекта;

• ArrayFire - многообещающий проект, который к сожалению не собрался на моей машине без проблем и я решил что не судьба;

• ViennaCL 1.7.1 - интересный проект с реализацией LU-алгоритма, и немалым количеством итеративных методов решения и написанный на относительно современном C++. Однако у него есть 2 минуса: во-первых реализован алгоритм LU, а не LUP и, во-вторых, нет поддержки комплексных чисел. К сожалению, проект ViennaCL остановился в развитии (текущий релиз вышел 8 лет назад, а последний коммит был сделан 5 лет назад). По реализации LUP даже создан Issue в 2013 году. С другой стороны, у них есть инфраструктура для запуска произвольных kernel's, возможно она будет удобна.

Изначально планировалось просто применить существующую библиотеку, но обзор имеющихся вариантов оставил негативное впечатление, а благодаря книге "OpenCL in Action" и серии статей на Habr'е, я понял что зверь не так страшен и можно попробовать переписать используемый алгоритм на OpenCL.

Реализация алгоритма

В качестве основы для собственной реализации алгоритма была использована обёртка OpenCL-Wrapper. Кое-что, на мой взгляд, там реализовано не совсем удачно, но для стартового проекта OpenCL отличное решение.

При реализации алгоритма возникали сложности. Так совсем не очевидным оказалось модель синхронизации потоков. Можно устанавливать барьер/семафор для потоков в определенной части программы, однако такой подход не поможет, если количество work item'ов больше чем количество вычислительных ядер. Причетельно что не возникает deadlock, а алгоритм просто тихо продолжается без полной синхронизации. Так на моей интегрированной видеокарте с 192 ядрами расхождение с эталонным решением начиналось на 193 шаге.

Вышесказанное заставило прибегнуть к другой модели синхронизации. Исходный алгоритм был разбит на отдельные kernels, которые запускались по очереди с различными параметрами. Алгоритм превратился в эдакий конвеер.

Алгоритм LUP-разложения без дополнительной памяти (in-place) распадается на 3 kernels. Во-первых ищем максимальный элемент в столбце:

__kernel void find_max_in_column_kernel(const uint k, const uint N, __global double* A, __global int* ip, __global int* ier) { 

    int m = k;
	for (int i = k + 1; i < N; ++i) {
		if (fabs(A[i * N + k]) > fabs(A[m * N + k])) m = i;
	}

	ip[k] = m ; // save row number of max element

	// swap rows if needed
	if (m != k) {
		ip[N-1] = -ip[N-1];

		// swap diagonal element A
		double t = A[m * N + k];
		A[m * N + k] = A[k * N + k];
		A[k * N + k] = t;

		if (t == 0.0) {
			*ier = k;
			ip[N-1] = 0;
			return;
		}
	}
}

Эта задача редукции и её эффективное распараллеливание это отдельный интересный вопрос. Но пока нас устроит если ею займётся одно вычислительное ядро GPU. Хоть CPU и быстрее GPU в однопоточном режиме, но данные уже загружены в видеопамять и переносить их туда-сюда для одного шага выглядит нецелесообразным. Затем делим все элементы в столбце ниже диагонального на величину этого диагонального элемента:

__kernel void identify_column_kernel(const uint k, const uint N, __global double* A, __global int* ip, __global int* ier) {
	
	const uint i = get_global_id(0);

    // equivalent to "for (uint i = k+1; i < N; i++)", but executed in parallel
	if (i >= k + 1 && i < N)
		A[i * N + k] /= - A[k * N + k];
}

И наконец вычитаем строку на диагонали из всех остальных строк

__kernel void columns_process_kernel(const uint k, const uint N, __global double* A, __global int* ip, __global int* ier) {
	
	const uint j = get_global_id(0);

    // equivalent to "for (uint i = k+1; i < N; i++)", but executed in parallel
	if (j >= k + 1 && j < N) {
		// swap row element
		double tmp = A[ip[k] * N + j];
		A[ip[k] * N + j] = A[k * N + j];
		A[k * N + j] = tmp;

		if (tmp != 0.0)
			for (int i = k + 1; i < N; ++i)
				A[i * N + j] += A[i * N + k] * tmp;
	}
}

Во всех kernels есть параметр k, обозначающий шаг алгоритма, все три kernels запускаются в цикле по k ∊ [0;N-1].

После разложения матрицы, её компоненты L, U и P используются для решения системы уравнений вида с плотной матрицей заполненной случайными значениями. Было выполнено несколько бенчмарков в которых сравнивалась производительность следущих методов:

• метод LUP-разложения на CPU, в данный момент используемый в проекте;

• метод LU-разложения на GPU в реализации ViennaCL;

• метод LUP-разложения на GPU собственной реализации.

Результаты

Дисклеймер: в процессе тестирования у меня возникли сложности с оборудованием, так что результаты нужно воспримать как ориентир. Так проект ViennaCL тестировался уже в самом конце и по нему можно сделать вывод только что их реализация имеет тот же порядок производительности что и моя реализация.

Результаты бенчмарков можно увидеть на рисунках ниже. Каждая точка запускалась по несколько раз, чтобы получить доверительный интервал. Как и ожидалось, решение системы уравнений на GPU требует меньше времени почти во всех случаях, за исключением маленьких матриц (меньше 1000х1000).

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

Исходный последовательный алгоритм имеет сложность O(N³). В параллельной версии, при N меньшем количества вычислительных ядер, один цикл исчезает (выполняется только одна итерации на каждом ядре) и сложность становиться O(N²). На интегрированной видеокарте с небольшим количеством ядер, когда каждое ядро выполняет несколько шагов цикла, была получена аппроксимация сложности ~О(N^2,5).

Для матриц большой размерности имеет смысл попробовать другой класс алгоритмов решения - итеративные решатели. В планах также было протестировать метод BiCGStab, однако по неизвестной пока причине он переставал сходиться уже на матрице 40х40 заполненной случайными значениями и в реалзиации ViennaCL и в реализации IML++ (Iterative Methods Library) v. 1.2a. Возможно требуется применение предобуславливателя (предварительные шаги, повышающие сходимость решения) или более продвинутых методов. А также я не обнаружил способа задания начального приближения в реализации ViennaCL, что снижает ценность этой реализации для моего проекта.

Заключение

Результаты эксперимента показали отличные возможности применения OpenCL. Дальнешие планы по развитию проекта предполагают рост размерности системы уравнений и без применения GPU или смены численного метода (что не факт что возможно) мне не обойтись.

Все найденные библиотеки (насколько я изучил) и полученные результаты касались систем уравнений с действительными числами. Потребуется также разработать реализацию алгоритма LUP-разложения матрицы с комплексными значениями.

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