1.4. Примеры использования пакета ScaLAPACK
В качестве первого примера использования пакета ScaLAPACK рассмотрим уже знакомую нам задачу перемножения матриц. Это позволит сопоставить получаемые результаты и производительность двух программ, созданных с использованием непосредственно среды параллельного программирования MPI [2, раздел 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
PARAMETER ( ONE = 1.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 Если процессор не вошел в сетку, то он ничего не делает;
C такое может случиться, если заказано, например, 5 процессоров
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 Инициализация дескрипторов для 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
C
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, a
C матрица B(I,J) = 1./J
C здесь имеются в виду глобальные индексы.
k = 0
do 250 j = 1,nq
jc = (j-1)/nb
jb = mod(j-1,nb)
jn = mycol*nb + jc*npcol*nb + jb + 1
do 240 i = 1,np
ic = (i-1)/nb
ib = mod(i-1,nb)
in = myrow*nb + ic*nprow*nb + ib + 1
k = k + 1
a(k) = dfloat(in)
b(k) = 1.D+0/dfloat(jn)
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
C
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
Пример 3. Решение системы линейных алгебраических уравнений с ленточной
матрицей
В данном примере решается система линейных алгебраических уравнений
с симметричной положительно определенной ленточной матрицей. Используются
подпрограммы PDPBTRF и PDPBTRS библиотеки ScaLAPACK. Матрица
A формируется следующим образом:
| 6 | -4 | 1 | 0 | 0 | 0 | . . . |
| -4 | 6 | -4 | 1 | 0 | 0 | . . . |
| 1 | -4 | 6 | -4 | 1 | 0 | . . . |
| 0 | 1 | -4 | 6 | -4 | 1 | . . . |
| 0 | 0 | 1 | -4 | 6 | -4 | . . . |
| 0 | 0 | 0 | 1 | -4 | 6 | . . . |
с хранением верхнего треугольника (параметр UPLO ='U')
| * | * | a13 | a24 | a35 | a46 | . . . | an-2,n |
| * | a12 | a23 | a34 | a45 | a56 | . . . | an-1,n |
| a11 | a22 | a33 | a44 | a55 | a66 | . . . | an,n |
Позиции, помеченные *, не используются.
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 Подпрограмма генерации матрицы 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
C
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
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