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

Пример составления головной программы
при чтении исходных матриц из файла и записи результатов в файл.

Содержимое файлов с исходными матрицами приведено следом за головной программой.

      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