Решите матричную систему уравнений с полосами

Мне нужно решить двумерное уравнение Пуассона, то есть систему уравнений для AX=B, где A — матрица n на n, а B — вектор n на 1. Будучи матрицей дискретизации для двумерной задачи Пуассона, я знаю, что только 5 диагоналей не будут нулевыми. Lapack не предоставляет функций для решения этой конкретной задачи, но имеет функции для решения системы уравнений с ленточной матрицей, а именно DGBTRF (для факторизации LU) и DGBTRS. Теперь 5 диагоналей: главная диагональ, первые диагонали выше и ниже главной и две диагонали выше и ниже на m диагоналей относительно главной диагонали. Прочитав документацию lapack о хранении диапазонов, я узнал, что мне нужно создать матрицу (3 * m + 1) на n для хранения A в формате хранения диапазонов, назовем эту матрицу AB. Теперь вопросы:

1) в чем разница между dgbtrs и dgbtrs_? Intel MKL предоставляет оба, но я не могу понять, почему

2) dgbtrf требует, чтобы матрица хранения полос была массивом. Должен ли я линеаризовать AB по строкам или по столбцам?

3) это правильный способ вызова двух функций?

int n, m;
double *AB;
/*... fill n, m, AB, with appropriate numbers */

int *pivots;
int nrows = 3 * m + 1, info, rhs = 1;
dgbtrf_(&n, &n, &m, &m, AB, &nrows, pivots, &info);
char trans = 'N';
dgbtrs_(&trans, &n, &m, &m, &rhs, AB, &nrows, pivots, B, &n, &info);

person Patrik    schedule 06.06.2012    source источник


Ответы (2)


  1. Он также предоставляет DGBTRS и DGBTRS_. Это fortran administrativa, о которых вам не стоит беспокоиться. Просто вызовите dgbtrs (причина в том, что на некоторых архитектурах имена подпрограмм на Фортране имеют добавленное подчеркивание, а на других нет, и имена могут быть как в верхнем, так и в нижнем регистре - Intel MKL #define подходит для dgbtrs).

  2. Подпрограммы LAPACK ожидают, что основные матрицы столбцов (т. е. в стиле Fortran): сохраняйте столбцы один за другим. Групповое хранилище, которое вы должны использовать, не сложно: http://www.netlib.org/lapack/lug/node124.html.

  3. Мне это кажется хорошим, но, пожалуйста, попробуйте его заранее на небольших задачах (кстати, это всегда хорошая идея). Также убедитесь, что вы обрабатываете ненулевые info (так сообщается об ошибках).

Лучший стиль - использовать MKL_INT вместо простого int, это typedef для правильного типа (может отличаться на некоторых архитектурах).

Также не забудьте выделить память для pivots перед вызовом dgbtrf.

person Alexandre C.    schedule 06.06.2012
comment
о № 2, не испортит ли это факторизацию LU? Я имею в виду, что нет возможности сообщить dgbtrf, что AB является основным. Кроме того, как мне установить главные размеры в двух случаях? - person Patrik; 07.06.2012
comment
Да, вы будете решать A'x = B. Позвольте мне исправить мой ответ по этому поводу. Для ведущего измерения они обычно равны количеству строк, если только вы не делаете что-то необычное (например, передаете подматрицы в качестве аргументов). Опять же, лучший способ решить проблемы LAPACK — это попробовать решить небольшие проблемы, которые вы можете решить вручную. - person Alexandre C.; 07.06.2012

Это может быть не по теме. Но для уравнения Пуассона решение на основе БПФ намного быстрее. Просто выполните 2D FFT вашего потенциального поля, разделенного на -(k ^ 2 + lambda ^ 2), затем выполните IFFT. лямбда - это небольшое число, чтобы избежать расхождения при k = 0. Уравнение с 5 диагональю представляет собой ограниченную по полосе аппроксимацию уравнения Пуассона, которая аппроксимирует дифференциальный оператор конечной разностью.

http://en.wikipedia.org/wiki/Screened_Poisson_equation

person fchen    schedule 28.08.2013