Содержимое файлов с исходными матрицами приведено следом за головной программой.
PROGRAM TDGESV25
*
* Тест к подпрограмме PDGESV
*
include 'mpif.h'
* Именованные константы - параметры задачи ..
INTEGER NOUT, MEMSIZE, NIN, NPR
PARAMETER ( MEMSIZE = 20000,
$ NOUT =13, NIN = 11, NPR = 6 )
*
* Локальные переменные ..
INTEGER CTXT, I, IAM, INFO, IPA, IPB, IPW,
$ IPPIV, MYCOL, MYROW, NP, NQ, NQRHS, NRHS,
$ N, NB, NPCOL, NPROCS, NPROW, IA, JA, IB, JB,
$ WORKSIZ
*
* Локальные массивы ..
INTEGER DESCA( 9 ), DESCB( 9 ), DESCX( 9)
DOUBLE PRECISION MEM( MEMSIZE )
*
* Внешние функции ..
INTEGER ICEIL, NUMROC
DOUBLE PRECISION PDLAMCH, PDLANGE
EXTERNAL ICEIL, NUMROC, PDLAMCH, PDLANGE
*
* Внешние подпрограммы ..
EXTERNAL BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
$ BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
$ BLACS_SETUP, DESCINIT, PDLAPRNT, PDLAREAD,
$ PDLAWRITE, PDGESV
*
* Выполняемые операторы ..
*
* Постановка задачи
*
N = 9
NB = 2
NPROW = 2
NPCOL = 2
IA = 1
JA = 1
IB = 1
JB = 1
NRHS = 1
*
* Инициализация пакета BLACS
*
CALL BLACS_PINFO( IAM, NPROCS )
*
IF( NPROCS.LT.1 ) THEN
CALL BLACS_SETUP( IAM, NPROW*NPCOL )
END IF
*
* Инициализация контекста BLACS'а и решетки процессов
*
CALL BLACS_GET( -1, 0, CTXT )
CALL BLACS_GRIDINIT( CTXT, 'R', NPROW, NPCOL )
CALL BLACS_GRIDINFO( CTXT, NPROW, NPCOL, MYROW, MYCOL )
*
* Если процесс не является частью данного контекста,
* переход на конец программы
*
IF( MYROW.GE.NPROW .OR. MYCOL.GE.NPCOL ) GO TO 20
*
* Вычисление размеров локальных частей матриц A и B на всех процессах
*
NP = NUMROC( N, NB, MYROW, 0, NPROW )
NQ = NUMROC( N, NB, MYCOL, 0, NPCOL )
NQRHS = NUMROC( NRHS, NB, MYCOL, 0, NPCOL )
*
* Инициализация дескрипторов для матриц A, B и X
*
CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, CTXT, MAX(1,NP), INFO)
CALL DESCINIT(DESCB, N, NRHS, NB,NB, 0,0, CTXT, MAX(1,NP), INFO)
CALL DESCINIT(DESCX, N, NRHS, NB,NB, 0,0, CTXT, MAX(1,NP), INFO)
*
* Распределение памяти MEM
*
* Указатель на первый элемент локальной части матрицы A
IPA = 1
* Указатель на первый элемент локальной части матрицы(вектора) B
IPB = IPA + DESCA( 9 )*NQ
* Указатель на первый элемент массива перестановок строк матрицы A
IPPIV = IPB + DESCB( 9 )*NQRHS
LIPIV = ICEIL( 4*(NP+NB), 8 )
* Указатель на рабочий массив
IPW = IPPIV + MAX( NP, LIPIV )
*
* Проверка величины рабочей памяти
*
WORKSIZ = MAX( NB, NP )
*
INFO = 0
IF(IPW + WORKSIZ .GT. MEMSIZE ) THEN
IF( IAM .EQ. 0 ) WRITE( NOUT, FMT = * )' test'
INFO = 1
END IF
*
CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, INFO, 1, -1, 0 )
IF( INFO .GT. 0 ) THEN
IF( IAM .EQ. 0 )
$ WRITE( NOUT, FMT = * ) 'BAD MEMORY'
GO TO 10
END IF
*
* Чтение из файлов и распределение матриц A и B
*
CALL PDLAREAD('tdgesv2A.dat', MEM(IPA), DESCA, 0, 0, MEM(IPW))
CALL PDLAREAD('tdgesv2B.dat', MEM(IPB), DESCB, 0, 0, MEM(IPW))
*
* Печать входных данных
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
WRITE( NPR, FMT = 92 )
WRITE( NPR, FMT = 98 )N, N, NB, NB
WRITE( NPR, FMT = 97 )NPROW*NPCOL, NPROW, NPCOL
WRITE( NPR, FMT = * ) ' Исходные данные '
WRITE( NPR, FMT = * )' IPA, IPB, IPPIV, IPW'
WRITE( NPR, FMT = 94 )IPA, IPB, IPPIV, IPW
WRITE( NPR, FMT = * ) ' DESCA'
WRITE( NPR, FMT = 85 )DESCA
WRITE( NPR, FMT = * ) ' DESCB'
WRITE( NPR, FMT = 85 )DESCB
WRITE( NPR, FMT = * ) ' DESCX'
WRITE( NPR, FMT = 85 )DESCX
WRITE( NPR, FMT = * ) ' Матрица A'
END IF
94 FORMAT( / 4I7 /)
85 FORMAT( / 9I5 /)
CALL PDLAPRNT( N, N, MEM(IPA), IA, JA, DESCA, 0, 0, 'A',
$ NPR, MEM( IPW ) )
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
WRITE( NPR, FMT = * ) ' Матрица B'
END IF
*
CALL PDLAPRNT( N, NRHS, MEM(IPB), IB, JB, DESCB, 0, 0, 'B',
$ NPR, MEM( IPW ) )
*
* Решение линейной системы A * X = B
*
CALL PDGESV( N, NRHS, MEM(IPA), IA, JA, DESCA, MEM(IPPIV),
$ MEM(IPB), IB, JB, DESCB, INFO )
*
* Печать результатов X
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
WRITE( NPR, FMT = * ) ' Результаты '
WRITE( NPR, FMT = 96 )INFO
WRITE( NPR, FMT = * ) ' Bектор X'
END IF
*
CALL PDLAPRNT( N, NRHS, MEM(IPB), IB, JB, DESCB, 0, 0, 'X', NPR,
$ MEM( IPW ) )
*
* Запись в файл результатов работы подпрограммы PDGESV.f
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
OPEN( NOUT, FILE = 'tdgesv2X.res', STATUS = 'UNKNOWN' )
WRITE( NOUT, FMT = * ) ' Решение линейной системы A * X = B'
END IF
*
CALL PDLAWRITE( 'tdgesv2X.res', N, NRHS, MEM(IPB), 1, 1, DESCB,
$ 0, 0, MEM(IPW) )
10 CONTINUE
CALL BLACS_GRIDEXIT( CTXT )
20 CONTINUE
IF( IAM .EQ. 0 ) THEN
IF( NOUT .NE. 6 .AND. NOUT .NE. 0 ) CLOSE( NOUT )
END IF
*
CALL BLACS_EXIT( 0 )
*
92 FORMAT( / 'Тест к подпрограмме PDGESV' )
98 FORMAT(/' Решение A X = B ,'//' где A - матрица ',I2,' на',I2,',',
$ /' разбитая на блоки', I2, ' на ', I2, /)
97 FORMAT( / ' Запуск на ', I2, ' процессах,',
$ ' образующих решетку размером ', I2, ' на ', I2 /)
96 FORMAT( / ' Значение INFO = ', I6 /)
STOP
END
Файл с исходной матрицей A:
9 9
0.19000D+02
-0.19000D+02
-0.19000D+02
-0.19000D+02
-0.19000D+02
-0.19000D+02
-0.19000D+02
-0.19000D+02
-0.19000D+02
0.30000D+01
0.30000D+01
-0.30000D+01
-0.30000D+01
-0.30000D+01
-0.30000D+01
-0.30000D+01
0.30000D+01
-0.30000D+01
0.10000D+01
0.10000D+01
0.10000D+01
-0.10000D+01
-0.10000D+01
-0.10000D+01
-0.10000D+01
-0.10000D+01
-0.10000D+01
0.12000D+02
0.12000D+02
0.12000D+02
0.12000D+02
-0.12000D+02
-0.12000D+02
-0.12000D+02
-0.12000D+02
-0.12000D+02
0.10000D+01
0.10000D+01
0.10000D+01
0.10000D+01
0.10000D+01
-0.10000D+01
-0.10000D+01
-0.10000D+01
-0.10000D+01
0.16000D+02
0.16000D+02
0.16000D+02
0.16000D+02
0.16000D+02
0.16000D+02
-0.16000D+02
-0.16000D+02
-0.16000D+02
0.10000D+01
0.10000D+01
0.10000D+01
0.10000D+01
0.10000D+01
0.10000D+01
0.10000D+01
-0.10000D+01
-0.10000D+01
0.30000D+01
0.30000D+01
0.30000D+01
0.30000D+01
0.30000D+01
0.30000D+01
0.30000D+01
0.30000D+01
-0.30000D+01
0.11000D+02
0.11000D+02
0.11000D+02
0.11000D+02
0.11000D+02
0.11000D+02
0.11000D+02
0.11000D+02
0.11000D+02
Файл с исходным вектором B
9 1
0.00000D+00
0.00000D+00
0.10000D+01
0.00000D+00
0.00000D+00
0.00000D+00
0.00000D+00
0.00000D+00
0.00000D+00
Ниже приводится содержимое файла с результатами работы подпрограммы
PDGESV.
Решение линейной системы A * X = B
9 1
-0.146081976924362695D-17
-0.166666666666666657D+00
0.500000000000000000D+00
0.000000000000000000D+00
0.000000000000000000D+00
0.173472347597680709D-17
-0.500000000000000000D+00
0.166666666666666657D+00
0.000000000000000000D+00