1.4. Примеры использования пакета ScaLAPACK

В качестве первого примера использования пакета ScaLAPACK рассмотрим уже знакомую нам задачу перемножения матриц. Это позволит сопоставить получаемые результаты и производительность двух программ, созданных с использованием непосредственно среды параллельного программирования MPI (часть II, раздел 8.2) и с помощью пакета ScaLAPACK.

Пример 1. Перемножение матриц

Используется подпрограмма PDGEMM из PBLAS, которая выполняет матричную операцию C = aA*B + bC, где A, В и С - матрицы, a и b - константы. В нашем случае мы полагаем a = 1, b = 0.

    program abcsl
    include 'mpif.h'

c Параметр nm определяет максимальную размерность блока матрицы
c на одном процессоре, массивы описаны как одномерные

    parameter (nm = 1000, nxn = nm*nm)
    double precision a(nxn), b(nxn), c(nxn), mem(nm)
    double precision time(6), ops, total, t1

c Параметр NOUT - номер выходного устройства (терминал)
    PARAMETER ( NOUT = 6 )
    DOUBLE PRECISION ONE, ZERO
    PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0)
    INTEGER DESCA(9), DESCB(9), DESCC(9)

c Инициализация BLACS
    CALL BLACS_PINFO( IAM, NPROCS )
c Вычисление формата сетки процессоров,
c наиболее приближенного к квадратному
    NPROW = INT(SQRT(REAL(NPROCS)))
    NPCOL = NPROCS/NPROW
c Считывание параметров решаемой задачи ( N - размер матриц и
c NB - размер блоков ) 0-м процессором и печать этой информации
    IF( IAM.EQ.0 ) THEN
    WRITE(*,* ) ' Input N and NB: '
    READ( *, * ) N, NB
    WRITE( NOUT, FMT = * )
    WRITE( NOUT, FMT = 9999 )
   $      'The following parameter values will be used:'
    WRITE( NOUT, FMT = 9998 ) 'N ', N
    WRITE( NOUT, FMT = 9998 ) 'NB ', NB
    WRITE( NOUT, FMT = 9998 ) 'P ', NPROW
    WRITE( NOUT, FMT = 9998 ) 'Q ', NPCOL
    WRITE( NOUT, FMT = * )
    END IF

c Рассылка считанной информации всем процессорам
    call MPI_BCAST(N, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
    call MPI_BCAST(NB,1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
c Теоретическое количество операций при перемножении
c двух квадратных матриц
    ops = (2.0d0*dfloat(n)-1)*dfloat(n)*dfloat(n)

c Инициализация сетки процессоров
    CALL BLACS_GET( -1, 0, ICTXT )
    CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL )
    CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )

c Если процессор не вошел в сетку, то он ничего не делает;
c такое может случиться, если заказано, например, 5 процессоров
    IF( MYROW.GE.NPROW .OR. MYCOL.GE.NPCOL )
   $ GO TO 500

c Вычисление реальных размеров матриц на процессоре
    NP = NUMROC( N, NB, MYROW, 0, NPROW )
    NQ = NUMROC( N, NB, MYCOL, 0, NPCOL )
c Инициализация дескрипторов для 3-х матриц
    CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, ICTXT, MAX(1,NP ), INFO )
    CALL DESCINIT( DESCB, N, N, NB, NB, 0, 0, ICTXT, MAX(1,NP ), INFO )
    CALL DESCINIT( DESCC, N, N, NB, NB, 0, 0, ICTXT, MAX(1,NP ), INFO )
*
    lda = DESCA(9)
c Вызов процедуры генерации матриц А и В
    call pmatgen(a, DESCA, np, nq, b, DESCB, nprow, npcol, myrow, mycol)
    
    t1 = MPI_Wtime()
*
* Вызов процедуры перемножения матриц
    CALL PDGEMM('N','N', N, N, N, ONE, A, 1, 1, DESCA,
   $             B, 1, 1, DESCB, 0. 0, C, 1, 1, DESCC)
*
    time(2) = MPI_Wtime() - t1
c Печать угловых элементов матрицы C
c с помощью служебной подпрограммы
    if (IAM.EQ.0) write(*,*) 'Matrix C...'
    CALL PDLAPRNT( 1, 1, C, 1, 1, DESCC, 0, 0, 'C', 6, MEM )
    CALL PDLAPRNT( 1, 1, C, 1, N, DESCC, 0, 0, 'C', 6, MEM )
    CALL PDLAPRNT( 1, 1, C, N, 1, DESCC, 0, 0, 'C', 6, MEM )
    CALL PDLAPRNT( 1, 1, C, N, N, DESCC, 0, 0, 'C', 6, MEM )
c Вычисление времени, затраченного на перемножение,
c и оценка производительности в Mflops.
    total = time(2)
    time(4) = ops/(1.0d6*total)
    if (IAM.EQ.0) then
    write(6,80) lda
 80 format(' times for array with leading dimension of',i4)
    write(6,110) time(2), time(4)
 110 format(2x,'Time calculation: ',f12.4, ' sec.',
  $        ' Mflops = ',f12.4)
    end if
c Закрытие BLACS-процессов
    CALL BLACS_GRIDEXIT( ICTXT )
    CALL BLACS_EXIT(0)
9998 FORMAT( 2X, A5, '   :   ', I6 )
9999 FORMAT( 2X, 60A )
 500 continue
    stop
    end

    subroutine pmatgen(a, DESCA, np, nq, b, DESCB, nprow, npcol, myrow, mycol)
    integer n, i, j, DESCA(*), DESCB(*), nprow, npcol, myrow, mycol
    double precision a(*), b(*)
c
    nb = DESCA(5)
c Заполнение локальных частей матриц A и B,
c матрица A формируется по алгоритму A(I,J) = I,
c матрица B(I,J) = 1./J
c Здесь имеются в виду глобальные индексы.
    k = 0
    do 250 j = 1,nq
    jc = (j-1)/nb
    jb = mod(j-1,nb)
    jg = mycol*nb + jc*npcol*nb + jb + 1
    do 240 i = 1,np
    ic = (i-1)/nb
    ib = mod(i-1,nb)
    ig = myrow*nb + ic*nprow*nb + ib + 1
    k = k + 1
    a(k) = dfloat(ig)
    b(k) = 1.D+0/dfloat(jg)
 240 continue
 250 continue
    return
    end

Пример 2. Решение системы линейных уравнений с матрицей общего вида

В данном примере решается система линейных уравнений с матрицей общего вида, которая формируется с помощью генератора случайных чисел. Правая часть системы формируется таким образом, чтобы получить единичное решение. Для решения системы используются две вычислительных подпрограммы: PDGETRF (для факторизации матрицы) и PDGETRS (для нахождения решения). Общий шаблон вызовов функций мало отличается от предыдущего примера. Отличие, главным образом, состоит в том, что в этом примере все коммуникационные операции выполняются с помощью подпрограмм из библиотеки BLACS (чисто из иллюстративных соображений), хотя многие операци можно компактнее записать на MPI.

    program pdlusl
    include 'mpif.h'
    parameter (nsize = 3000, nxn = nsize*nsize)
    double precision a(nxn), b(nsize), x(nsize)
    double precision time(6), cray, ops, total, norma, normx, t1
    double precision resid, residn, eps, epslon, rab, RANN
    integer ipvt(nsize), iwork(5), init(4)
    PARAMETER ( NOUT = 6 )
    DOUBLE PRECISION ONE
    PARAMETER ( ONE = 1.0D+0 )
    INTEGER DESCA( 9 ), DESCB( 9 ), DESCX( 9 )
c Параметр NRHS - количество правых частей
    NRHS = 1
    CALL BLACS_PINFO( IAM, NPROCS )
*
    NPROW = INT(SQRT(REAL(NPROCS)))
    NPCOL = NPROCS/NPROW
c Теоретическое число операций, которое необходимо выполнить для
c решения системы
    ops = (2.0d0*dfloat(n)**3)/3.0d0 + 2.0d0*dfloat(n)**2

c Формирование одномерной сетки процессоров
c (только для рассылки параметров)
    CALL BLACS_GET( -1, 0, ICTXT )
    CALL BLACS_GRIDINIT( ICTXT, 'Row-major', 1, NPROCS )
c На 0-м процессоре считываем параметры, упаковываем в массив
c и рассылаем всем остальным с помощью процедуры IGEBS2D.
    IF( IAM.EQ.0 ) THEN
    WRITE( *,* ) ' Input N and NB: '
    READ ( *,* ) N, NB
    IWORK( 1 ) = N
    IWORK( 2 ) = NB
    CALL IGEBS2D( ICTXT, 'All', ' ', 2, 1, IWORK, 2 )

    WRITE( NOUT, FMT = 9999 )
  $         'The following parameter values will be used:'
    WRITE( NOUT, FMT = 9998 ) 'N ', N
    WRITE( NOUT, FMT = 9998 ) 'NB ', NB
    WRITE( NOUT, FMT = 9998 ) 'P ', NPROW
    WRITE( NOUT, FMT = 9998 ) 'Q ', NPCOL
    WRITE( NOUT, FMT = * )
*
    ELSE
c На не 0-м процессоре получаем массив с процессора (0,0) и
c распаковываем его.
      CALL IGEBR2D( ICTXT, 'All', ' ', 2, 1, IWORK, 2, 0, 0 )
      N  = IWORK( 1 )
      NB = IWORK( 2 )
    END IF

c Уничтожаем временную сетку процессоров
    CALL BLACS_GRIDEXIT( ICTXT )
c Формируем рабочую сетку процессоров
    CALL BLACS_GET( -1, 0, ICTXT )
    CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL )
    CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
c
c Проверка процессоров, не вошедших в сетку
    IF( MYROW.GE.NPROW .OR. MYCOL.GE.NPCOL )
  $    GO TO 500
c Определяем точное число строк и столбцов распределенной матрицы
c в процессоре
    NP = NUMROC( N, NB, MYROW, 0, NPROW )
    NQ = NUMROC( N, NB, MYCOL, 0, NPCOL )
c
c Формируем дескрипторы
*
    CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, ICTXT, MAX(1,NP ), INFO )
    CALL DESCINIT( DESCB, N, NRHS, NB, NB, 0, 0, ICTXT, MAX(1,NP), INFO )
    CALL DESCINIT( DESCX, N, NRHS, NB, NB, 0, 0, ICTXT, MAX(1,NP), INFO )
*
    lda = DESCA(9)
c Вызов процедуры генерации матрицы A и вектора B
    call pmatgenl(a, DESCA, NP, NQ, b, DESCB, nprow, npcol, myrow, mycol)
    t1 = MPI_Wtime()
*
* Обращение к процедере факторизации матрицы A
    CALL PDGETRF(N, N, A, 1, 1, DESCA, ipvt, INFO )
*
    time(1) = MPI_Wtime() - t1
    t1 = MPI_Wtime()
c
c Обращение к процедуре решения системы уравнений с факторизованной
c матрицей
    CALL PDGETRS('No', N, NRHS, A, 1, 1, DESCA, ipvt, B, 1, 1, DESCB, INFO)
*
    time(2) = MPI_Wtime() - t1
    total = time(1) + time(2)
c
c На этом собственно решение задачи заканчивается, далее идет
c подготовка печати и сама печать
    if (iam.eq.0) then
    write(6,40)
 40 format( '      x(1)         x(nр)')
    write(6,50) x(1),x(np)
 50 format(1p5e16.8)
c
    write(6,60) n
 60 format(//'   times are reported for matrices of order ',i5)
    write(6,70)
 70 format(6x,'factor',5x,'solve',6x,'total',5x,'mflops',7x,'unit',
 $         6x,'ratio')
c
    time(3) = total
    time(4) = ops/(1.0d6*total)
    time(5) = 2.0d0/time(4)
    time(6) = total/cray
    write(6,80) lda
 80 format(' times for array with leading dimension of',i4)
    write(6,110) (time(i),i=1,6)
 110 format(6(1pe11.3))
    write(6,*)' end of test'
    end if
c
    CALL BLACS_GRIDEXIT( ICTXT )
    CALL BLACS_EXIT(0)
 9998 FORMAT( 2X, A5, 'аа :аа', I6 )
 9999 FORMAT(2X, 60A )
 500  continue
    stop
    end

c Процедура генерации матрицы A с помощью генератора случайных
c чисел RANN.
c Последовательности случайных чисел на процессорах должны быть
c независимые, чтобы матрица не оказалась вырожденной.
    subroutine pmatgenl(a, DESCA, NP, NQ, b, DESCB, nprow, npcol, myrow, mycol)
    integer n, init(4), i, j, DESCA(*), DESCB(*), nprow, npcol, myrow, mycol
    double precision a(*),b(*),rann
c
    nb = DESCA(5)
    ICTXT = DESCA(2)
c Инициализация генератора случайных чисел
    init(1) = 1
    init(2) = myrow + 2
    init(3) = mycol + 3
    init(4) = 1325
c Заполнение матрицы A
    k = 0
    do 250 j = 1,nq
    do 240 i = 1,np
    k = k + 1
    a(k) = rann(init) - 0.5
 240 continue
 250 continue
    na = k
c Вычисление вектора B такого, чтобы получить единичное решение,
c сначала вычисляются локальные суммы по строке на каждом процессоре,
c а затем выполняется суммирование по всем процессорам.
    do 350 i = 1,np
    k = i
    b(i) = 0
    do 351 j = 1,nq
    b(i) = b(i) + a(k)а
    k = k + np
 351 continue
    CALL BLACS_BARRIER( ICTXT, 'All' )аа
    CALL DGSUM2D( ICTXT, 'Row', ' ', 1, 1, b(i), 1, -1, 0)
 350 continue
    return
    end

В данном примере решается система линейных алгебраических уравнений с симметричной положительно определенной ленточной матрицей. Используются подпрограммы PDPBTRF и PDPBTRS библиотеки ScaLAPACK. Матрица A формируется следующим образом:

с хранением верхнего треугольника (параметр UPLO ='U'), где звездочкой помечены неиспользуемые позиции:

    program bandu
c nsize - максимальное число столбцов матрицы в одном процессоре
    parameter (nsize = 30000)
    double precision a(3*nsize), b(nsize), x(nsize), bg(nsize)
    double precision AF(3*nsize), WORK(10)
    integer ipvt(nsize), BW
    PARAMETER ( NOUT=6 )
    INTEGER DESCA(7), DESCB(7), DESCX(7)
    CALL BLACS_PINFO( IAM, NPROCS )

c Задаем размеры матрицы и сетки процессоров
    N = 9000
    NRHS = 1
    NPROW = 1
    NPCOL = 4
c BW - ширина ленты над диагональю
    BW=2
c Вычисляем NB - длину блока матрицы A
    NDD = mod(N,NPCOL)
    IF(NDD.EQ.0) THEN
    NB = N/NPCOL
    ELSE
    NB = N/NPCOL + 1
    END IF
    NB = MAX(NB,2*BW)
    NB = MIN(N,NB)
    IF(IAM.EQ.0) THEN
    WRITE( 6, FMT = 9998 ) 'N ', N
    WRITE( 6, FMT = 9998 ) 'NRHS ', NRHS
    WRITE( 6, FMT = 9998 ) 'BW   ', BW
    WRITE( 6, FMT = 9998 ) 'P ', NPROW
    WRITE( 6, FMT = 9998 ) 'Q ', NPCOL
    WRITE( 6, FMT = 9998 ) 'NB ', NB
    END IF
c Вычисляем размеры рабочих массивов
    LWORK = 2*BW*BW
    LAF = (NB+2*BW)*BW
c
c Инициализируем сетку процессоров
    CALL BLACS_GET( -1, 0, ICTXT )
    CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL )
    CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
    IF( MYROW.GE.NPROW .OR. MYCOL.GE.NPCOL )
 $  GO TO 500
c
    NP = NUMROC( (BW+1), (BW+1), MYROW, 0, NPROW )
    NQ = NUMROC( N, NB, MYCOL, 0, NPCOL )
c Формируем дескрипторы для левых и правых частей уравнения
    DESCA(1) = 501
    DESCA(2) = ICTXT
    DESCA(3) = N
    DESCA(4) = NB
    DESCA(5) = 0
    DESCA(6) = BW+1
    DESCA(7) = 0
c
    DESCB(1) = 502
    DESCB(2) = ICTXT
    DESCB(3) = N
    DESCB(4) = NB
    DESCB(5) = 0
    DESCB(6) = NB
    DESCB(7) = 0
c
    lda = NB
c Вызов подпрограммы генерации ленточной матрицы и правой части
    call pmatgenb(a, DESCA, bw, b, DESCB, nprow, npcol,
 $                myrow, mycol, n, bg)
c Факторизация матрицы
    CALL PDPBTRF('U', N, BW, A, 1, DESCA, AF, LAF, WORK, LWORK, INFO3 )
c Решение системы
    CALL PDPBTRS('U', N, BW, NRHS, A, 1, DESCA, B, 1, DESCB,
 $                AF, LAF, WORK, LWORK, INFO)
    if (iam.eq.0) then
    write(6,40)
 40 format('x(1),...,x(4)')
    write(6,50) (b(i),i=1,4)
 50 format(4d16.8)
    end if
    CALL BLACS_GRIDEXIT( ICTXT )
    CALL BLACS_GRIDEXIT( ICTXTB
    CALL BLACS_EXIT (0)
500 continue
    stop
    end


c Подпрограмма-функция генерации (I,J)-го матричного элемента
    double precision function matij(i,j)
    double precision rab
    rab = 0.0d0
    if(i.eq.j)               rab =  6.0d0
    if(i.eq.j+1.or.j.eq.i+1) rab = -4.0d0
    if(i.eq.j+2.or.j.eq.i+2) rab =  1.0d0
    matij = rab
    return
    end

c Подпрограмма генерации матрицы A и вектора В
    subroutine pmatgenb(a, DESCA, bw, b, DESCB, nprow,npcol,
 $                      myrow, mycol, n, bg)
    integer i, j, DESCA(*), DESCB(*), nprow, npcol, bw, bw1, myrow, mycol
    double precision a(bw+1,*), b(*), bg(*), matij

    nb = DESCA(4)
    ICTXT = DESCA(2)
    n = DESCA(3)
    BW1 = BW + 1
c Генерация всех компонент вектора B таким образом,
c чтобы решение X(I) = I
    do 231 i = 1,n
    bg(i) = 0.0
    n1 = max(1,i-bw)
    n2 = min(n,i+bw)
    do 231 j = n1,n2
    bg(i) = bg(i) + matij(i,j)*j
 231 continue
c
c Вычисление локальной части матрицы A
    jcs = MYCOL*NB
    NC  = MIN(NB,N-jcs)
    do 250 j = 1,NC
    jc = jcs + j
    do 240 i = 1,BW1
    ic = jc - BW1 + i
    if (ic.ge.1 )   a(i,j) = matij(ic,jc)
 240 continue
 250 continue
c Заполнение локальной части вектора B
    do 350 i = 1,NC
    b(i) = bg(jcs+i)
 351 continue
 350 continue
    return
    end