пакетное решение CUDA с разреженной лентой Ax = b для различных b's

Question

У меня есть разреженная матрица A, и я бы хотел (направлять) решение Ax = b. У меня около 500 векторов b, поэтому я бы хотел решить для соответствующих 500 x. Я новичок в CUDA, поэтому я немного смущен относительно того, какие варианты у меня есть.

cuSOLVER имеет пакетный прямой решатель cuSolverSP для разреженного A_i x_i = b_i, используя QR here. (С тобой все будет хорошо, так как А прилично обучен.) Однако, насколько я могу судить, я не могу использовать тот факт, что все мои A_i одинаковы.

Может ли альтернативный вариант сначала определить разреженную факторизацию LU (QR) на процессоре или графическом процессоре, а затем параллельно выполнять обратную подстановку (соответственно, backsub и matrix mult) на GPU? Если cusolverSp< t >csrlsvlu() для одного b_i, есть ли стандартный способ выполнить эту операцию для нескольких b_i?

Наконец, поскольку у меня нет интуиции для этого, следует ли ожидать ускорения на GPU для любого из этих вариантов, учитывая необходимые накладные расходы? x имеет длину ~ 10000-100000. Благодарю.

Answer 1

В настоящее время я работаю над чем-то подобным. Я решил в основном обернуть конъюгированный градиент и 0-й уровень неполных холеских предустановленных сопряженных градиентов-решателей, которые были собраны вместе с CUDA SDK в небольшой класс.

Вы можете найти их в каталоге CUDA_HOME под пути: samples/7_CUDALibraries/conjugateGradient и /Developer/NVIDIA/CUDA-samples/7_CUDALibraries/conjugateGradientPrecond

В принципе, вы бы загрузить матрицу в память устройства один раз (и ICCG, вычислить соответствующий кондиционер/матричный анализ), то вызовите ядро ​​с различными векторами b.

Я не знаю, что вы ожидаете от своей структуры матричной полосы, но если она симметрична и диагонально доминирует (по диагональным полосам вдоль каждой строки и столбца противоположный знак диагонали, а их сумма меньше, чем диагональная запись) или положительно определенный (без собственных векторов с собственным значением 0), тогда CG и ICCG должны быть полезными. В качестве альтернативы, различные многосеточные алгоритмы - это еще один вариант, если вы хотите пройти их кодирование.

Если ваша матрица является только положительной полуопределенной (например, имеет по крайней мере один собственный вектор с собственным значением нуля), вы все равно можете уйти с использованием CG или ICCG, если вы гарантируете, что: 1) Правая рука side (b векторы) ортогональны нулевому пространству (пустое пространство, означающее собственные векторы с собственным значением нуля). 2) Решение, которое вы получаете, ортогонально нулевому пространству.

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

И, наконец, если ваша матрица имеет Circulant Band structure, вы можете рассмотреть возможность использования решения на основе быстрого преобразования Фурье (FFT). Численные решатели на основе FFT часто могут обеспечить превосходную производительность в тех случаях, когда они применимы.

Answer 2

Если вы не возражаете идти с библиотекой с открытым исходным кодом, вы можете также проверить бугор: CUSP Quick Start Page

Он имеет довольно приличный набор решателей, в том числе нескольких способов: предварительно кондиционированных CUSP Preconditioner Examples

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