Решение системы линейных уравнений с ленточной матрицей общего вида методом Гаусса
Подпрограмма вычисляет решение вещественной системы линейных уравнений
A(1 : N, JA : JA+N-1) * X = B( IB : IB+N-1, 1 : NRHS) ,
где A(1 : N, JA : JA+N-1) - квадратная вещественная ленточная
распределенная матрица порядка N
с шириной ленты BWL, BWU.
Используется метод исключения Гаусса с выбором ведущего элемента для разложения матрицы на множители L и U (нижняя и верхняя треугольные матрицы).
Матрицы A и B должны быть заданы в виде массивов двойной точности. Все вычисления проводятся в режиме DOUBLE PRECISION.
Литература:
http://www.netlib.org/scalapack/slug/index.html
A. Cleary and J. Dongarra, Implementation in scalapack of divide - and - conquer algorithms for banded and tridiagonal linear systems, Computer Science Dept. Technical Report CS - 97 - 358, University of Tennessee, Knoxville, TN, April 1997. (Also LAPACK Working Note #125).
CALL PDGBSV ( N, BWL, BWU, NRHS, A, JA, DESCA, IPIV,
B, IB, DESCB, WORK, LWORK, INFO)
Параметры
| N - | порядок распределенной подматрицы A (1 : N, JA : JA + N - 1), N ≥ 0 (глобальный входной параметр, тип целый) |
| BWL - | число поддиагоналей, 0 ≤ BWL ≤ N - 1 (глобальный входной параметр, тип целый) |
| BWU - | число наддиагоналей, 0 ≤ BWU ≤ N - 1 (глобальный входной параметр, тип целый) |
| NRHS - |
число правых частей, т.е. число столбцов распределенной подматрицы B ( IB : IB + N - 1, 1 : NRHS), NRHS ≥ 0 (глобальный входной параметр, тип целый) |
| A - |
указатель на локальную память, занимаемую массивом
двойной точности с первой размерностью LLD_A ≥ (2 * bwl + 2 * bwu + 1) (хранится в DESCA); на входе - этот массив содержит локальные части факторизуемой распределенной матрицы A (1 : N, JA : JA + N - 1). Эта локальная порция хранится в запакованном ленточном формате, используемом в пакете LAPACK. Более подробно см. в разделе документации "Схема размещения в локальной памяти и блочно - столбцовое разбиение ленточных матриц"; на выходе - массив А содержит информацию о деталях разложения матрицы; заметим, что перестановки выполняются над матрицей так, что выдаваемые множители отличаются от множителей, выдаваемых LAPACK'ом (локальный входной и локальный выходной параметр) |
| JA - | номер столбца в глобальном массиве A, указывающий на начало матрицы, которая должна быть обработана (которая может быть либо всей матрицей A, либо подматрицей матрицы A) (глобальный входной параметр, тип целый) |
| DESCA - |
дескриптор распределенной матрицы A (одномерный массив длины DLEN); если тип дескриптора одномерный ( DTYPE_A = 501), DLEN ≥ 7, если тип дескриптора двумерный ( DTYPE_A = 1), DLEN ≥ 9; Содержит информацию о размещении A в памяти. Полное описание DESCA см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр тип целый) |
| IPIV - | одномерный массив длины ≥ DESCA (NB). Содержит индексы выбираемых ведущих элементов в процессе локального разложения. Пользователи не должны (не имеют права) изменять его содержимое между процессом разложения матрицы и процессом решения системы |
| B - |
указатель на локальную память, занимаемую массивом
двойной точности с локальной ведущей размерностью
lld_b ≥ NB; на входе - этот массив содержит локальные куски правой части B ( IB : IB + N - 1, 1 : NRHS); на выходе - массив содержит локальный кусок решений распределенной матрицы X, столбцы которой образуют NRHS векторов решений исходной системы (локальный входной и локальный выходной параметр) |
| IB - | номер строки в глобальном массиве B, указывающий на первую строку матрицы, которая должна быть обработана (которая может быть либо всей матрицей B, либо подматрицей матрицы B) (глобальный входной параметр, тип целый) |
| DESCB - |
дескриптор распределенной матрицы B (одномерный массив длины DLEN); если тип дескриптора одномерный ( DTYPE_B = 502), DLEN ≥ 7, если тип дескриптора двумерный ( DTYPE_B = 1), DLEN ≥ 9; Содержит информацию о размещении B в памяти. Полное описание DESCB см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый) |
| WORK - | временное рабочее пространство двойной точности. Это пространство может изменять свое содержимое между вызовами подпрограмм. WORK должно иметь длину, задаваемую параметром LWORK; на выходе - WORK (1) содержит минимальное значение LWORK (локальное рабочее пространство, локальный выходной параметр) |
| LWORK - |
размер входного рабочего пространства пользователя с
именем WORK. Если LWORK слишком мало, минимальный
приемлемый (допустимый) размер будет возвращен в
элементе WORK (1) и будет выдан код ошибки. LWORK ≥ (NB + bwu) * (bwl + bwu) + (bwl + bwu) * (bwl + 2 * bwu) + max (NRHS * (NB + 2 * bwl + 4 * bwu),1) (локальный или глобальный входной параметр, тип целый) |
| INFO - | целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр) |
| = 0 - | успешное завершение работы; |
| < 0 - |
если i - ый фактический параметр является массивом
и его j - ый элемент имеет недопустимое значение, то INFO = - ( i * 100 + j ), если i - ый фактический параметр является скаляром и имеет недопустимое значение, то INFO = - i; |
| > 0 - |
если INFO = K ≤ NPROCS,
подматрица, хранящаяся
на процессоре INFO и подвергающаяся разложению
локально, не является подматрицей с диагональным
преобладанием, и процесс разложения не был завершен; если INFO = K > NPROCS, то подматрица, хранящаяся на процессоре INFO-NPROCS, представляющая взаимодействие с другими процессорами, не была разложена, и процесс разложения не был завершен. |
Версии
| PSGBSV - PCGBSV PZGBSV | решение системы с ленточной матрицей общего вида методом Гаусса с выбором ведущих элементов для вещественных данных одинарной точности, комплексных данных одинарной точности и комплексных данных двойной точности соответственно |
Вызываемые подпрограммы
Здесь указаны только базовые подпрограммы (2 - ого уровня), которые вызываются из целевой подпрограммы (1 - ого уровня).
| PDGBTRF - PZGBTRF | LU - разложение ленточной матрицы общего вида методом Гаусса с выбором ведущего элемента по столбцу |
| PDGBTRS - PZGBTRS | Решение системы с ленточной матрицей общего вида на основе LU - разложения, полученного подпрограммой PDGBTRF(PZGBTRF) |
Замечания по использованию
| В подпрограммах PSGBSV, PCGBSV, PZGBSV параметры A, B и WORK имеют тип REAL, COMPLEX и DOUBLE COMPLEX, соответственно | |
| Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS), PXERBLA ( из пакета PBLAS) | |
| Для того чтобы избежать дублирования, большинство проверок ошибок в этой подпрограмме не выполняется, эти проверки выполняются в подпрограммах PDGBTRF(PZGBTRF) и PDGBTRS(PZGBTRS) | |
| Имеются следующие ограничения на входные параметры. |
| 1. |
Не циклическое ограничение P * NB ≥ mod ( JA - 1, NB) + N |
| 2. | Ни один из процессоров не может содержать больше одного блока глобальной матрицы. |
| 3. |
Размер блока (NB) не может быть слишком маленьким: если матрица обрабатывается более чем на одном процессоре, то должно быть выполнено следующее неравенство: NB ≥ 2*BW Это объясняется тем, что на каждом процессоре действия выполняются над блоками матрицы размера O (NB). Если этот размер слишком мал, то алгоритм "разделяй и властвуй" является неэффективным. |
| 4. | Чтобы избежать лишних передач сообщений между процессорами, при работе с подматрицами следует полагать JA = IB. |
Решается система уравнений A X = B , где A - квадратная ленточная матрица порядка 9
Пусть матрица А системы имеет вид:
| 2 | -1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
| -1 | 2 | -1 | 0 | 0 | 0 | 0 | 0 | 0 | |
| 0 | -1 | 2 | -1 | 0 | 0 | 0 | 0 | 0 | |
| 0 | 0 | -1 | 2 | -1 | 0 | 0 | 0 | 0 | |
| 0 | 0 | 0 | -1 | 2 | -1 | 0 | 0 | 0 | |
| 0 | 0 | 0 | 0 | -1 | 2 | -1 | 0 | 0 | |
| 0 | 0 | 0 | 0 | 0 | -1 | 2 | -1 | 0 | |
| 0 | 0 | 0 | 0 | 0 | 0 | -1 | 2 | -1 | |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 1 |
Вектор правых частей В имеет вид:
| 0 |
| 0 |
| 0 |
| 0 |
| 0 |
| 0 |
| 0 |
| 0 |
| 1 |
Матрица разбивается по столбцам на 3 блока с размером блоков NB = 3.
Осуществляется пропуск программы на 3 процессах
и используется одномерная решетка процессов 1 на 3.
Фрагмент фортранного текста вызывающей программы
Полный текст теста можно получить в tpdgbsv.zip, а принятые обозначения - посмотреть в разделе документации "Обозначения и упрощения в примерах по использованию подпрограмм комплекса"
Другой (более удобный) вариант теста к этой подпрограмме с использованием служебной подпрограммы PDELSET для распределения элементов исходной глобальной матрицы по параллельным процессам (как это описано в разделах методического пособия в пп.12, 15, 16) и с использованием служебной подпрограммы PDLAPRNT (для выборки всех элементов распределенной матрицы/вектора со всех процессов и последующей распечатки их единым массивом) можно получить в tdgbsv1.zip.
PROGRAM TPDGBSV
include 'mpif.h'
INTEGER DLEN_, JA, IB, N, NB, RSRC, CSRC, MXLLDA,
$ MXLLDB, NRHS, NBRHS, NOUT, MXLOCR, MXLOCC,
$ MXRHSC, BWL, BWU, LWORK
PARAMETER ( DLEN_ = 7, JA = 1, IB = 1,
$ N = 9, NB = 3, RSRC = 0, CSRC = 0,
$ MXLLDA = 5, MXLLDB = 3, NRHS = 1, NBRHS = 1,
$ NOUT = 6, MXLOCR = 5, MXLOCC = 3, MXRHSC = 1,
$ BWU = 1, BWL = 1, LWORK = 53 )
INTEGER ICTXT, INFO, MYCOL, MYROW, NPCOL, NPROW
DOUBLE PRECISION ANORM, BNORM, EPS, RESID, XNORM
INTEGER DESCA( DLEN_ ), DESCB( DLEN_ ),
$ IPIV( MXLOCR )
DOUBLE PRECISION A( MXLLDA, MXLOCC ), B( MXLLDB, MXRHSC ),
$ WORK( LWORK )
EXTERNAL BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
$ MATINIT, PDGBSV, SL_INIT
DATA NPROW / 1 / , NPCOL / 3 /
CALL SL_INIT( ICTXT, NPROW, NPCOL )
CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
IF( MYROW.EQ.-1 ) GO TO 10
DESCA(1) = 501
DESCA(2) = ICTXT
DESCA(3) = N
DESCA(4) = NB
DESCA(5) = CSRC
DESCA(6) = 5
DESCA(7) = 0
DESCB(1) = 502
DESCB(2) = ICTXT
DESCB(3) = N
DESCB(4) = NB
DESCB(5) = RSRC
DESCB(6) = 3
DESCB(7) = 0
CALL MATINIT( A, DESCA, B, DESCB )
CALL PDGBSV( N, BWL, BWU, NRHS, A, JA, DESCA, IPIV,
$ B, IB, DESCB, WORK, LWORK, INFO)
CALL BLACS_GRIDEXIT( ICTXT )
10 CONTINUE
CALL BLACS_EXIT( 0 )
STOP
END
SUBROUTINE MATINIT( A, DESCA, B, DESCB )
* Подпрограмма генерации и распределения по процессам матриц A и B
INTEGER DESCA( * ), DESCB( * )
DOUBLE PRECISION A( * ), B( * )
INTEGER CTXT_
PARAMETER ( CTXT_ = 2 )
INTEGER ICTXT, MYCOL, MYROW, NPCOL, NPROW
DOUBLE PRECISION MO, T, Z
EXTERNAL BLACS_GRIDINFO
ICTXT = DESCA( CTXT_ )
CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
T = 2.0D0
Z = 0.0D0
MO = -1.0D0
* Локальные части матриц A и B для процесса (0, 0)
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
A( 1 ) = Z
A( 2 ) = Z
A( 3 ) = Z
A( 4 ) = T
A( 5 ) = MO
A( 6 ) = Z
A( 7 ) = Z
A( 8 ) = MO
A( 9 ) = T
A( 10 ) = MO
A( 11 ) = Z
A( 12 ) = Z
A( 13 ) = MO
A( 14 ) = T
A( 15 ) = MO
B( 1 ) = Z
B( 2 ) = Z
B( 3 ) = Z
* Локальные части матриц A и B для процесса (0, 1)
ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.1 ) THEN
A( 1 ) = Z
A( 2 ) = Z
A( 3 ) = MO
A( 4 ) = T
A( 5 ) = MO
A( 6 ) = Z
A( 7 ) = Z
A( 8 ) = MO
A( 9 ) = T
A( 10 ) = MO
A( 11 ) = Z
A( 12 ) = Z
A( 13 ) = MO
A( 14 ) = T
A( 15 ) = MO
B( 1 ) = Z
B( 2 ) = Z
B( 3 ) = Z
* Локальные части матриц A и B для процесса (0, 2)
ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.2 ) THEN
A( 1 ) = Z
A( 2 ) = Z
A( 3 ) = MO
A( 4 ) = T
A( 5 ) = MO
A( 6 ) = Z
A( 7 ) = Z
A( 8 ) = MO
A( 9 ) = T
A( 10 ) = MO
A( 11 ) = Z
A( 12 ) = Z
A( 13 ) = MO
A( 14 ) = 1.0D0
A( 15 ) = Z
B( 1 ) = Z
B( 2 ) = Z
B( 3 ) = 1.0D0
END IF
RETURN
END
Результаты:
Решение системы
X = ( 1.D0, 2.D0, 3.D0, 4.D0, 5.D0, 6.D0, 7.D0, 8.D0, 9.D0 )
Значение INFO = 0