Решение системы A * X = B или AT * X = B с матрицей общего вида методом Гаусса с выбором ведущего элемента по столбцу, с итерационным уточнением решения и оценкой границ ошибок
Подпрограмма вычисляет решение вещественной системы линейных уравнений
(1) sub(A) * X = sub(B) или
(2) sub(A)T * X = sub(B) ,
где sub(A) = A ( IA : IA+N-1, JA : JA+N-1) - квадратная распределенная матрица
порядка N,
X - распределенная матрица размеров N на NRHS,
sub(B) = B ( IB : IB+N-1, JB : JB+NRHS-1) - распределенная матрица
размеров N на NRHS.
Используется LU - разложение методом Гаусса с выбором ведущего элемента по столбцу для факторизации матрицы sub(A) в виде
sub(A) = P * L * U ,
где P - матрица перестановок,
L - нижняя треугональная матрица с единичными диагональными элементами,
U - верхняя треугональная матрица.
Матрицы L и U помещаются в памяти на место матрицы sub (AF). Полученное LU - разложение используется для решения системы sub (A) * X = sub (B) или sub (A) T * X = sub (B).
Итерационное уточнение решения и невязка вычисляются с той же точностью, которую имеют входные данные.
Матрицы A, AF, B и X должны быть заданы в виде массивов двойной точности. Все вычисления проводятся в режиме DOUBLE PRECISION.
Литература
http://www.netlib.org/scalapack/slug/index.html
CALL PDGESV3 ( TRANS, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF,
DESCAF, IPIV, B, IB, JB, DESCB, X, IX, JX, DESCX, FERR,
BERR, WORK, LWORK, IWORK, LIWORK, INFO)
Параметры
| TRANS - |
если TRANS = ' N ', решается система (1), если TRANS = ' T ', решается система (2), (глобальный входной параметр, тип символьный); |
| N - | порядок распределенной подматрицы sub (A), N ≥ 0 (глобальный входной параметр, тип целый); |
| NRHS - | число правых частей, т.е. число столбцов распределенной подматрицы sub (B), NRHS ≥ 0 (глобальный входной параметр, тип целый); |
| A - |
указатель на локальную память, занимаемую массивом двойной
точности размерности
( LLD_A, LOCc ( JA + N - 1)); это - локальные части факторизуемой распределенной матрицы sub(A) порядка N (локальный входной параметр); |
| IA - | номер строки в глобальном массиве A, указывающий на первую строку подматрицы sub (A) (глобальный входной параметр, тип целый); |
| JA - | номер столбца в глобальном массиве A, указывающий на первый столбец подматрицы sub (A) (глобальный входной параметр, тип целый); |
| DESCA - | дескриптор распределенной матрицы A (одномерный массив длины DLEN_); подробное описание см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый); |
| AF - | указатель на локальную память, занимаемую массивом двойной точности размерности ( LLD_AF, LOCc ( JA + N - 1)); на выходе - массив AF содержит локальные части распределенных множителей L и U, полученных при разложении (факторизации) подматрицы в виде sub (A) = P * L * U; единичные диагональные элементы матрицы L в памяти не хранятся (локальный выходной параметр); |
| IAF - | номер строки в глобальном массиве AF, указывающий на первую строку подматрицы sub (AF) (глобальный входной параметр, тип целый); |
| JAF - | номер столбца в глобальном массиве AF, указывающий на первый столбец подматрицы sub (AF) (глобальный входной параметр, тип целый); |
| DESCAF - | дескриптор распределенной матрицы AF (одномерный массив длины DLEN_); подробное описание см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый); |
| IPIV - |
одномерный массив длины ( LOCr (N_AF) + NB_AF),
содержащий информацию о выборе ведущего элемента на каждом шаге
исключения Гаусса; элемент массива IPIV ( i ) содержит глобальный номер строки, которая переставляется со строкой, имеющей локальный номер i; этот массив содержит информацию о перестановках строк распределенной матрицы A (локальный выходной параметр, тип целый); |
| B - | указатель на локальную память, занимаемую массивом двойной точности размерности ( LLD_B, LOCc ( JB + NRHS - 1)), в котором задаются локальные части распределенной матрицы sub (B) правой части (локальный входной параметр); |
| IB - | номер строки в глобальном массиве B, указывающий на первую строку подматрицы sub (B) (глобальный входной параметр, тип целый); |
| JB - | номер столбца в глобальном массиве B, указывающий на первый столбец подматрицы sub (B) (глобальный входной параметр, тип целый); |
| DESCB - | дескриптор распределенной матрицы B (одномерный массив длины DLEN_); подробное описание см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый); |
| X - | указатель на локальную память, занимаемую массивом двойной точности размерности ( LLD_X, LOCc ( JX + NRHS - 1)); на выходе (если INFO = 0) - это распределенная матрица X, столбцы которой образуют NRHS векторов решений исходной системы, полученных в результате итерационного уточнения (локальный выходной параметр); |
| IX - | номер строки в глобальном массиве X, указывающий на первую строку подматрицы sub ( X ) (глобальный входной параметр, тип целый); |
| JX - | номер столбца в глобальном массиве X, указывающий на первый столбец подматрицы sub ( X ) (глобальный входной параметр, тип целый); |
| DESCX - | дескриптор распределенной матрицы X (одномерный массив длины DLEN_); подробное описание см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый); |
| FERR - | одномерный массив длины LOCc ( N_B), каждый элемент FERR ( j ) которого содержит оценку границы относительной ошибки вычисленного вектора решений X ( j ) (т.е. j - ого столбца матрицы X ( IX : IX + N - 1, JX : JX + NRHS - 1) ), полученную в соответствии с описанием в http://www.netlib.org/scalapack/slug/, раздел 4.6 (см. также [12, 25, 27, 28]) (локальный выходной параметр, тип DOUBLE PRECISION ); |
| BERR - | одномерный массив длины LOCc ( N_B); содержащий покомпонентные оценки границ при обратном анализе относительных ошибок для каждого вектора решений X ( j ), т.е. наименьшее относительное возмущение в любом элементе матрицы A ( IA : IA + N - 1, JA : JA + N - 1) или B ( IB : IB + N - 1, JB : JB + NRHS - 1), которое делает X ( j ) точным решением (см. описание в http://www.netlib.org/scalapack/slug/, раздел 4.6 а также [12, 25, 27, 28]) (локальный выходной параметр, тип DOUBLE PRECISION ); |
| WORK - | одномерный рабочий массив двойной точности длины LWORK; на выходе - WORK (1) содержит минимальную требуемую величину LWORK (локальное рабочее пространство, локальный выходной параметр); |
| LWORK - |
размер рабочего пространства (массива) WORK; LWORK - локальный
входной параметр и должен быть не меньше, чем LWORK = PDGERFS(LWORK) + LOCr(N_A), где PDGERFS(LWORK) означает величину рабочего массива WORK, которая требуется при работе вызываемой подпрограммы PDGERFS (см. ниже); смысл обозначений LOCr (...) и LOCc (...) можно найти в разделах документации "Схема размещения в локальной памяти и блочно - циклическое отображение плотных матриц" и "Пример блочно - циклического распределения плотной матрицы по решетке процессов"; минимальная требуемая величина LWORK может быть вычислена самой подпрограммой, если к ней обратиться со значением LWORK = - 1; при этом подпрограмма не производит никаких других вычислений, а требуемое значение LWORK возвращается в элементе массива WORK (1); при этом, сообщений об ошибках не выдается (локальный или глобальный входной параметр); |
| IWORK - | одномерный рабочий массив целого типа длины LIWORK; на выходе IWORK (1) содержит минимальную требуемую величину LIWORK (локальное рабочее пространство, локальный выходной параметр); |
| LIWORK - |
размер рабочего пространства (массива) IWORK;
LIWORK - локальный входной параметр и должен быть не меньше,
чем LIWORK = LOCr ( N_A); минимальная требуемая величина LIWORK может быть вычислена самой подпрограммой, если к ней обратиться со значением LIWORK = - 1; при этом подпрограмма не производит никаких других вычислений, а требуемое значение LIWORK возвращается в элементе массива IWORK (1); при этом, сообщений об ошибках не выдается (локальный или глобальный входной параметр, тип INTEGER); |
| INFO - | целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр) |
| = 0 - | успешное завершение работы; |
| < 0 - |
если i - ый фактический параметр является массивом
и его j - ый элемент имеет недопустимое значение, то INFO = - ( i * 100 + j ), если i - ый фактический параметр является скаляром и имеет недопустимое значение, то INFO = - i; |
| > 0 - | если INFO = K, то элемент U( IA+K-1, JA+K-1) является точным нулем; разложение завершено, но матрица U вырождена и поэтому решение системы и оценки границ ошибок не вычисляются. |
Версии
| PSGESV3 - PCGESV3 PZGESV3 | решение системы с матрицей общего вида методом Гаусса с выбором ведущего элемента по столбцу для вещественных данных одинарной точности, комплексных данных одинарной точности и комплексных данных двойной точности соответственно |
Вызываемые подпрограммы
Здесь указаны только базовые подпрограммы (2 - ого уровня), которые вызываются из целевой подпрограммы (1 - ого уровня).
| PDGETRF - PZGETRF | LU - разложение матрицы общего вида методом Гаусса с выбором ведущего элемента по столбцу |
| PDGETRS - PZGETRS | Решение системы A X = B, AT X = B или AH X = B на основе LU - разложения, полученного подпрограммой PDGETRF |
| PDGERFS - PZGERFS | Выполнение итерационного уточнения решения системы A X = B, AT X = B или AH X = B, полученного подпрограммой PDGETRS (PZGETRS), и оценка границ ошибок полученного решения |
Замечания по использованию
| 1. | В подпрограммах PSGESV3, PCGESV3, PZGESV3 параметры A, B, AF, X, FERR и BERR имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно | |
| 2. | В версии подпрограммы PZGESV3 вместо параметров IWORK и LIWORK целого типа используются параметры с именами RWORK и LRWORK. При этом параметр RWORK имеет тип DOUBLE PRECISION, а параметр LRWORK ≥ 2 * LOCc (N_A). | |
| 3. | Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS), LSAME ( из пакета BLAS), PXERBA ( из пакета PBLAS), ICEIL,CHK1MAT, PCHK2MAT, INDXG2P, INFOG2L, NUMROC, ILCM ( из библиотеки ScaLAPACK_TOOLS), PDLACPY ( из пакета ScaLAPACK) |
Решение системы A X = B , где A - матрица порядка 9
Пусть матрица А системы имеет вид:
| 19 | 3 | 1 | 12 | 1 | 16 | 1 | 3 | 11 | |
| -19 | 3 | 1 | 12 | 1 | 16 | 1 | 3 | 11 | |
| -19 | -3 | 1 | 12 | 1 | 16 | 1 | 3 | 11 | |
| -19 | -3 | -1 | 12 | 1 | 16 | 1 | 3 | 11 | |
| -19 | -3 | -1 | -12 | 1 | 16 | 1 | 3 | 11 | |
| -19 | -3 | -1 | -12 | -1 | 16 | 1 | 3 | 11 | |
| -19 | -3 | -1 | -12 | -1 | -16 | 1 | 3 | 11 | |
| -19 | 3 | -1 | -12 | -1 | -16 | -1 | 3 | 11 | |
| -19 | -3 | -1 | -12 | -1 | -16 | -1 | -3 | 11 |
Вектор правых частей В имеет вид:
|
0 0 1 0 0 0 0 0 0 |
Матрица разбивается на блоки размеров 2 на 2.
Осуществляется пропуск программы на 6 процессах
и используется двумерная решетка процессов 2 на 3.
Ниже показано, как разбитая на блоки матрица A распределяется
(блочно-циклически) по решетке процессов 2 * 3.
| 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 | |
Распределение матрицы B (блочно-циклическое) по процессам
| 0 | 1 | 2 | |
| 0 |
0 0 | ||
|
0 0 | |||
| 0 | |||
| 1 |
1 0 | ||
|
0 0 |
Фрагмент фортранного текста вызывающей программы
Полный текст теста можно получить в
tdgesv3.zip,
а принятые обозначения - посмотреть в разделе документации
"Обозначения и упрощения в примерах
по использованию подпрограмм комплекса";
в тесте используются: служебная подпрограмма
PDELSET (для распределения элементов
глобальной матрицы/вектора по паралельным процессам) и служебная подпрограмма
PDLAPRNT (для выборки и распечатки элементов
распределенной матрицы/вектора единым массивом)
PROGRAM TDGESV3
include 'mpif.h'
INTEGER DLEN_, IA, JA, IB, JB, N, NB, RSRC, NRHS, NBRHS, NOUT,
$ IAF, JAF, IX, JX, CSRC, MXLLDA, MXLLDAF, MXLLDB, MXLLDX,
$ MXLOCR, MXLOCC, MXRHSC, LWORK, LIWORK
PARAMETER ( DLEN_ = 9, IA = 1, JA = 1, IB = 1, JB = 1,
$ IAF = 1, JAF = 1, IX = 1, JX = 1, N = 9, NB = 2, RSRC = 0,
$ CSRC = 0, MXLLDA = 5, MXLLDB = 5, MXLLDAF = 5,
$ MXLLDX = 5, NRHS = 1, NBRHS = 1, NOUT = 6, MXLOCR = 5,
$ MXLOCC = 4, MXRHSC = 1, LWORK = 26, LIWORK = 5 )
CHARACTER TRANS
PARAMETER ( TRANS = 'N' )
INTEGER ICTXT, INFO, MYCOL, MYROW, NPCOL, NPROW
DOUBLE PRECISION ANORM, BNORM, EPS, RESID, XNORM
INTEGER DESCA( DLEN_ ), DESCB( DLEN_ ), DESCAF( DLEN_ ),
$ DESCX( DLEN_ ), IPIV( MXLOCR+NB ), IWORK( LIWORK )
DOUBLE PRECISION A( MXLLDA, MXLOCC ), B( MXLLDB, MXRHSC ),
$ AF( MXLLDAF, MXLOCC ), X( MXLLDX, MXRHSC ),
$ WORK( LWORK ), FERR( MXRHSC ), BERR( MXRHSC )
EXTERNAL BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
$ DESCINIT, PDMATINIT, PDGESV3, SL_INIT
DATA NPROW / 2 / , NPCOL / 3 /
CALL SL_INIT( ICTXT, NPROW, NPCOL )
CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
IF( MYROW.EQ.-1 ) GO TO 10
CALL DESCINIT( DESCA, N, N, NB, NB, RSRC, CSRC, ICTXT, MXLLDA,
$ INFO )
CALL DESCINIT( DESCAF, N, N, NB, NB, RSRC, CSRC, ICTXT, MXLLDAF,
$ INFO )
CALL DESCINIT( DESCB, N, NRHS, NB, NBRHS, RSRC, CSRC, ICTXT,
$ MXLLDB, INFO )
CALL DESCINIT( DESCX, N, NRHS, NB, NBRHS, RSRC, CSRC, ICTXT,
$ MXLLDX, INFO )
CALL PDMATINIT( N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )
CALL PDGESV3( TRANS, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF, DESCAF,
$ IPIV, B, IB, JB, DESCB, X, IX, JX, DESCX, FERR, BERR,
$ WORK, LWORK, IWORK, LIWORK, INFO )
CALL BLACS_GRIDEXIT( ICTXT )
10 CONTINUE
CALL BLACS_EXIT( 0 )
STOP
END
SUBROUTINE PDMATINIT( N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )
*
* PDMATINIT генерирует и распределяет матрицы A и B
* по решетке процессов 2 x 3
* Параметры
INTEGER IA, IB, INFO, JA, JB, N
INTEGER DESCA( * ), DESCB( * )
DOUBLE PRECISION A( * ), B( * )
*
* Локальные переменные
INTEGER I, J
DOUBLE PRECISION AIJ, BIJ
*
* Внешние подпрограммы
EXTERNAL PDELSET
*
* Выполняемые операторы
INFO = 0
*
IF( IA.NE.1 ) THEN
INFO = -3
ELSE IF(JA .NE.1) THEN
INFO = -4
END IF
*
AIJ = 19.0D0
*
DO 20 J = 1, N
DO 10 I = 1, N
IF( I.LE.J ) THEN
CALL PDELSET( A, I, J, DESCA, AIJ )
ELSE
CALL PDELSET( A, I, J ,DESCA, -AIJ )
END IF
10 CONTINUE
IF( J.EQ.1 .OR. J.EQ.7 ) AIJ = 3.0D0
IF( J.EQ.2 .OR. J.EQ.4 .OR. J.EQ.6 ) AIJ = 1.0D0
IF( J.EQ.3 ) AIJ = 12.0D0
IF( J.EQ.5 ) AIJ = 16.0D0
IF( J.EQ.8 ) AIJ = 11.0D0
20 CONTINUE
CALL PDELSET( A, 8, 2, DESCA, 3.0D0 )
*
BIJ = 0.0D0
*
J = 1
DO 30 I = 1, N
CALL PDELSET( B, I, J, DESCB, BIJ )
30 CONTINUE
CALL PDELSET( B, 3, 1, DESCB, 1.0D0 )
*
RETURN
END
Результаты:
Решение системы
X = ( 0.D0, - 0.16666666666667D0, 0.5D0, 0.D0, 0.D0, 0.D0, - 0.5D0, 0.16666666666667D0, 0.D0 )
Значение INFO = 0
FERR = 0.56066263E-14
BERR = 0.13877788E-16