Приложение 2_2

Пример составления головной программы
для использования целевой подпрограммы из Приложения 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