Разбор примера по использованию подпрограммы для решения системы уравнений с плотной матрицей

Рассмотрим подробнее пример по использованию подпрограммы PDGESV, которая предназначена для решения линейной системы алгебраических уравнений с плотной матрицей общего вида  A. Данный пример приводится в описании этой подпрограммы в разделе "Пример использования" ( PDGESV.htm ).

Пример демонстрирует решение системы уравнений  A * X = B с матрицей  A размера 9 * 9. Ниже показана матричная форма записи этой системы.

  19   3    1   12    1   16   1   3  11 x1
0
-19   3    1    12    1   16   1   3  11 x2
0
-19  -3    1    12    1   16   1   3  11 x3
1
-19  -3  -1    12    1   16   1   3  11 x4
0
-19  -3  -1  -12    1   16   1   3  11 x5
= 0
-19  -3  -1  -12  -1   16   1    3  11 x6
0
-19  -3  -1  -12  -1  -16    1   3  11 x7
0
-19   3  -1  -12  -1  -16  -1   3  11 x8
0
-19  -3  -1  -12  -1  -16  -1  -3  11 x9 0

Для простоты решается система с одной правой частью ( NRHS = 1), т.е. правая часть  B системы - вектор (матрица из одного столбца).

Сначала мы должны задать решетку процессов и распределить по ней исходные глобальные матрицы  A и  B.

Предположим, что матрица  A разделена на блоки размером  MB * NB, где  MB = NB = 2. Выбрана решетка процессов размером  2 * 3, т.е.  NPROW = 2 и  NPCOL = 3. Таким образом, имеет место случай распределения матрицы  A, приведенный в разделе документации "Пример блочно - циклического распределения плотной матрицы по решетке процессов". После подстановки конкретных значений элементов матрицы  A получаем следующую картинку распределения ее элементов по процессам. ( Вверху указаны номера столбцов решетки процессов, слева - номера строк решетки процессов.)

0 1 2
0  19   3
-19   3
 1   3
 1   3
 1   12
 1   12
11
11
 1   16
 1   16
-19  -3
-19  -3
 1   3
 1   3
-1  -12
-1  -12
11
11
 1   16
-1   16
-19  -3 -1  -3 -1  -12 11 -1  -16
1 -19  -3
-19  -3
 1   3
 1   3
 1   12
-1   12
11
11
 1   16
 1   16
-19  -3
-19   3
 1   3
-1   3
-1  -12
-1  -12
11
11
-1  -16
-1  -16

Из последнего рисунка видно, что процесс с координатами (0, 0) содержит локальный массив (часть матрицы  A) размера (5, 4).

На следующем рисунке показано разбиение на блоки и распределение по той же решетке процессов исходного вектора правых частей  B.

0 1 2
0 b1
b2
b5
b6
b9
1 b3
b4
b7
b8

Как видим, блоки состоят из двух соседних элементов вектора  B и распределяются только в столбце решетки процессов с координатой 0. Процессы в других столбцах решетки не содержат никаких частей исходного глобального вектора  B.

После вызова подпрограммы  PDGESV и выполнения всех необходимых вычислений процесс (0, 0) будет содержать на месте локальной части вектора  B локальную часть глобального вектора решений  X. Ниже изображено, на месте каких элементов локального вектора  B (обозначенных как ~bi ) с локальными индексами, какие элементы глобального вектора  X  с глобальными индексами оказываются расположенными по окончании вычислений. Это расположение соответствует расположению элементов исходного глобального вектора правых частей  B, показанному на предыдущем рисунке. Ниже показаны также вычисленные значения элементов вектора  X.

x1 -> ~b1 0
x2 -> ~b2 -1/6
x5 -> ~b3 = 0
x6 -> ~b4 0
x9 -> ~b5 0

Аналогичное соответствие и полученные результаты изображены далее для процесса с координатами (1, 0).

x3 -> ~b1 1/2
x4 -> ~b2 0
x7 -> ~b3 = -1/2
x8 -> ~b4 1/6

Таким образом глобальный выходной вектор решений  X имеет следующий вид:

x1 0
x2 -1/6
x3 1/2
x4 0
x5 = 0
x6 0
x7 -1/2
x8 1/6
x9 0

Кроме того, выполняется проверка точности полученных результатов вычислением нормализованной невязки по формуле:

                            || A * x - b ||
                     ______________________   .
                     ( || x || * || A || * eps * N )

Здесь  eps  означает машинное эпсилон, которое вычисляется с помощью вспомогательной подпрограммы  PDLAMCH.