Быстрый алгоритм для краевой задачи

Question

Я ищу самый быстрый способ решения следующей задачи:

Учитывая некоторый объем точек решетки в трехмерной сетке, некоторые точки b_i (граница) удовлетворяют f(b_i)=0, тогда как другая точкаa_0 удовлетворяет f(a_0)= 1.

Все остальные точки (не граничные) представляют собой некоторую линейную комбинацию окружающих 26 точек.Например, я мог бы хотеть

f(i, j, k)= .05*f(i+1, j, k)+.1*f(i, j+1, k)+.07*f(i, j, k+1)+...

Сумма коэффициентов .05+.1+.07+... составит в сумме 1.Моя цель - найти для f(x_i) для всех x_i в объеме.

В настоящее время я использую метод последовательной избыточной релаксации (SOR), который в основном инициализирует границу объема, присваивает каждой точке средневзвешенное значение 26 окружающих точек и повторяет.Метод SOR принимает комбинацию f(x_i) после самой последней итерации и f(x_i) после предыдущей итерации.

Мне было интересно, если кто-нибудь знает какие-либо более быстрые способы решить эту проблему для 3D-сетки размером примерно 102x102x48.SOR в настоящее время требуется около 500-1000 итераций, чтобы достичь желаемого уровня (варьируется в зависимости от используемых коэффициентов).Я больше всего хочу использовать Matlab, IDL и C ++.Кто-нибудь имеет представление о том, как быстро SOR сравнивается с преобразованием задачи в линейную систему и использованием матричных методов (например, BCGSTAB)?Кроме того, какой метод будет наиболее эффективно (и легко) распараллеливать?У меня есть доступ к 250-ядерному кластеру, и я пытался сделать метод SOR параллельным, используя mpi и c ++, но не наблюдал такого увеличения скорости, как хотелось бы (идеал был бы порядка 100 раз).Буду очень признателен за любые идеи по ускорению решения этой проблемы.Благодарю.

Answer 1

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

Методы сопряженного градиента обычно сходятся быстрее, чем SOR (если я правильно помню,примерно на порядок).Я никогда не использовал BCGSTAB, но я помню, как GMRES хорошо работал над несимметричными задачами, и они, вероятно, могут извлечь выгоду из предварительной обработки.

Что касается возможностей для распараллеливания, большинству методов типа CG требуется только вычислениематрично-векторное произведение A*x, поэтому вам никогда не нужно формировать полную матрицу.Вероятно, это самая большая стоимость каждой итерации, и именно здесь вы должны взглянуть на многопоточность.

Две хорошие ссылки на это: Голуб и Ван Лоан и Трефетен и Бау .Последнее гораздо читабельнее ИМХО.

Надеюсь, это поможет ...

Answer 2

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

Допустим, у вас есть N баллов в вашем объеме. Давайте назовем p1...pN эти точки. То, что вы ищете, это f(p1)...f(pN)

Давайте возьмем точку pi. У него 26 соседей. Давайте назовем их pi-1...pi-26. Я полагаю, у вас есть доступ для каждой точки pi к коэффициентам линейной комбинации. Назовем их ci-1...ci-j...ci-26 (для "коэффициента линейной комбинации для точки pi ее j-го соседа)

Давайте сделаем еще лучше, вы можете считать, что pi - это линейная комбинация всех точек в вашем пространстве, за исключением того, что большинство из них (кроме 26) равны 0. У вас есть коэффициенты ci-1...ci-N.

Теперь вы можете построить большую матрицу N*N из этих ci-j коэффициентов:

+---------------------    -------+   
|  0  | c1-2 | c1-3 | ... | c1-N |   |f(p1)|    |f(p1)|
+---------------------    -------+   |     |    |     |
| c2_1|   0  | c2-3 | ... | c1-N |   |f(p2)|    |f(p2)|
+---------------------    -------+ * .     . =  .     .
|                                |   .     .    .     .
.                                .   .     .    .     .
+---------------------    -------+   .     .    .     .
| cN-1| cN-2 | cN-3 | ... |   0  |   |f(pN)|    |f(pN)|
+---------------------    -------+

Удивительно! Решение, которое вы ищете, это один из собственных векторов, соответствующих собственному значению 1!

Используйте оптимизированную матричную библиотеку, которая эффективно вычисляет собственные векторы (для разреженной матрицы) и надеюсь, что она уже распараллелена!

Редактировать: Забавно, я просто перечитал ваш пост, кажется, я только что дал вам ваше решение BCGSTAB ... Извините! ^^

Повторное редактирование: На самом деле, я не уверен, вы говорите о "Biconjugate Gradient Stabilized Method"? Потому что я не вижу линейного метода, о котором вы говорите, если вы делаете градиентный спуск ...

Повторное редактирование: я люблю математику =) Я даже могу доказать, что при условии sum of the ci = 1, что 1 на самом деле является собственным значением. Я не знаю размерности соответствующего пространства, но оно как минимум 1!