Пример составления головной программы
для использования целевой подпрограммы из Приложения 2_1.
PROGRAM TPDGESV
*
* Тест к подпрограмме PDGESV
*
include 'mpif.h'
* Именованные константы - параметры задачи
INTEGER DLEN_, IA, JA, IB, JB, M, N, MB, NB, RSRC,
$ CSRC, MXLLDA, MXLLDB, NRHS, NBRHS, NOUT,
$ MXLOCR, MXLOCC, MXRHSC
PARAMETER ( DLEN_ = 9, IA = 1, JA = 1, IB = 1, JB = 1,
$ M = 9, N = 9, MB = 2, NB = 2, RSRC = 0,
$ CSRC = 0, MXLLDA = 5, MXLLDB = 5, NRHS = 1,
$ NBRHS = 1, NOUT = 6, MXLOCR = 5, MXLOCC = 4,
$ MXRHSC = 1 )
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
*
* Локальные переменные
INTEGER ICTXT, INFO, MYCOL, MYROW, NPCOL, NPROW
DOUBLE PRECISION ANORM, BNORM, EPS, RESID, XNORM
*
* Локальные массивы
INTEGER DESCA( DLEN_ ), DESCB( DLEN_ ),
$ IPIV( MXLOCR+NB )
DOUBLE PRECISION A( MXLLDA, MXLOCC ), A0( MXLLDA, MXLOCC ),
$ B( MXLLDB, MXRHSC ), B0( MXLLDB, MXRHSC ),
$ WORK( MXLOCR ), WORK1( MB )
*
* Внешние функции
DOUBLE PRECISION PDLAMCH, PDLANGE
EXTERNAL PDLAMCH, PDLANGE
*
* Внешние подпрограммы
EXTERNAL BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
$ DESCINIT, PDMATINIT, PDGEMM, PDGESV, PDLACPY,
$ PDLAPRNT, SL_INIT
*
* Внутренние функции
INTRINSIC DBLE
*
* Операторы DATA
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
*
* Распределение матриц по решетке процессов
*
* Инициализация дескрипторов для матриц A и B
*
CALL DESCINIT( DESCA, M, N, MB, NB, RSRC, CSRC, ICTXT, MXLLDA,
$ INFO )
CALL DESCINIT( DESCB, N, NRHS, NB, NBRHS, RSRC, CSRC, ICTXT,
$ MXLLDB, INFO )
*
* Генерация матриц A и B и распределение их по решетке процессов
*
CALL PDMATINIT( N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )
*
* Копирование A и B для проверки невязки
*
CALL PDLACPY( 'All', N, N, A, 1, 1, DESCA, A0, 1, 1, DESCA )
CALL PDLACPY( 'All', N, NRHS, B, 1, 1, DESCB, B0, 1, 1, DESCB )
*
* Печать входных данных
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
WRITE( NOUT, FMT = 99 )
WRITE( NOUT, FMT = 98 )M, N, MB, NB
WRITE( NOUT, FMT = 97 )NPROW*NPCOL, NPROW, NPCOL
WRITE( NOUT, FMT = * ) ' Исходные данные '
WRITE( NOUT, FMT = * ) ' DESCA'
WRITE( NOUT, FMT = 85 )DESCA
WRITE( NOUT, FMT = * ) ' DESCB'
WRITE( NOUT, FMT = 85 )DESCB
WRITE( NOUT, FMT = * ) 'Матрица A'
END IF
85 FORMAT( / 9I5 /)
*
CALL PDLAPRNT( M, N, A, IA, JA, DESCA, 0, 0, 'A', NOUT, WORK1 )
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
WRITE( NOUT, FMT = * ) 'Матрица B'
END IF
CALL PDLAPRNT( N, NRHS, B, IB, JB, DESCB, 0, 0, 'B', NOUT, WORK1 )
*
* Вызов подпрограмм комплекса
*
* Решение линейной системы A * X = B
*
CALL PDGESV( N, NRHS, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB,
$ INFO )
*
* Печать результатов X
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
WRITE( NOUT, FMT = * ) ' Результаты '
WRITE( NOUT, FMT = 96 )INFO
WRITE( NOUT, FMT = * ) 'Решение X'
END IF
*
CALL PDLAPRNT( N, NRHS, B, IB, JB, DESCB, 0, 0, 'X', NOUT, WORK1 )
*
* Вычисление невязки ||A * X - B|| / ( ||X|| * ||A|| * eps * N )
*
EPS = PDLAMCH( ICTXT, 'Epsilon' )
ANORM = PDLANGE( 'I', N, N, A, 1, 1, DESCA, WORK )
BNORM = PDLANGE( 'I', N, NRHS, B, 1, 1, DESCB, WORK )
CALL PDGEMM( 'N', 'N', N, NRHS, N, ONE, A0, 1, 1, DESCA, B, 1, 1,
$ DESCB, -ONE, B0, 1, 1, DESCB )
XNORM = PDLANGE( 'I', N, NRHS, B0, 1, 1, DESCB, WORK )
RESID = XNORM / ( ANORM*BNORM*EPS*DBLE( N ) )
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
IF( RESID.LT.10.0D+0 ) THEN
WRITE( NOUT, FMT = 95 )
WRITE( NOUT, FMT = 93 )RESID
ELSE
WRITE( NOUT, FMT = 94 )
WRITE( NOUT, FMT = 93 )RESID
END IF
END IF
*
* Освобождение решетки процессов
*
* Освобождение контекста BLACS'а
*
CALL BLACS_GRIDEXIT( ICTXT )
10 CONTINUE
*
* Выход из BLACS
*
CALL BLACS_EXIT( 0 )
*
99 FORMAT( / ' Тест к подпрограмме PDGESV' )
98 FORMAT(/' Решение A X = B ,'//' где A - матрица ',I2,' на',I2,',',
$ /' разбитая на блоки', I2, ' на ', I2, /)
97 FORMAT( ' Пропуск на ', I2, ' процессах,',
$ ' образующих решетку размером ', I2, ' на ', I2 /)
96 FORMAT( / ' Значение INFO = ', I6 /)
95 FORMAT( /' Невязка мала')
94 FORMAT( /' Невязка велика')
93 FORMAT( / '||A*x - b|| / ( ||x||*||A||*eps*N ) = ', 1P, E16.8 )
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 CTXT_
PARAMETER ( CTXT_ = 2 )
*
INTEGER ICTXT, I, J, MYCOL, MYROW, NPCOL, NPROW
DOUBLE PRECISION AIJ
EXTERNAL BLACS_GRIDINFO, PDELSET
*
INFO = 0
ICTXT = DESCA( CTXT_ )
CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
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