В заключение настоящего параграфа рассмотрим пример программы, в которой Aztec используется для решения системы уравнений, возникающей при конечно-разностной аппроксимации 3-х мерного оператора Лапласа. Для простоты мы не будем использовать систему уравнений, связанную с решением какой-либо реальной физической задачи, а вместо этого будем формировать правую часть таким образом, чтобы получить заранее известное решение, например, x(i) = i + 1.
Программа azt.f
program azt
include 'mpif.h'
include 'az_aztecf.h'
integer NULL
parameter (NULL = 0)
integer n
common /global/ n
integer i,j,IAM,NPROCS,ierr,nz,K,in,ik,nz1,nzz
integer proc_config(0:AZ_PROC_SIZE), options(0:AZ_OPTIONS_SIZE)
double precision params(0:AZ_PARAMS_SIZE)
double precision status(0:AZ_STATUS_SIZE)
double precision s,time,total
integer N_update,ierror,recb(0:63),disp(0:63)
CHARACTER*8 met(0:4)
integer, dimension(:), allocatable:: data_org, update
integer, dimension(:), allocatable:: external, update_index
integer, dimension(:), allocatable:: extern_index, bindx
real(KIND=8), dimension(:), allocatable:: b,x,tmp,val
CALL MPI_Init(ierror)
met(0) = 'cg'
met(1) = 'gmres'
met(2) = 'cgs'
met(3) = 'tfqmr'
met(4) = 'bicgstab'
C Инициализация коммуникационного окружения
call AZ_processor_info(proc_config)
IAM = proc_config(AZ_node)
NPROCS = proc_config (AZ_N_procs)
if (IAM.eq.0) then
write(6,1)
1 format(' '//
$ ' Victor N. Datsyuk'/
$ ' High Performance Computers Department'/
$ ' Rostov State University'/
$ ' Rostov-on-Don, phone 219-97-13'//
$ ' email: root@rsusu1.rnd.runnet.ru'/
$ ' Internet: http://rsusu1.rnd.runnet.ru'/)
write(6,2)
2 format(' '/
$ ' Aztec example driver'//
$ ' This program create an MSR matrix corresponding to'/
$ ' a 7pt discrete approximation to the 3D Poisson operator'/
$ ' on an n x n x n square and solve it various methods.'//)
open(10, FILE='size.dat', FORM='FORMATTED')
read(10,*) n
endif
с Рассылка параметров задачи на все процессоры. Можно использовать
с либо процедуру из Aztec, либо из MPI
c call AZ_broadcast(n,4,proc_config,AZ_PACK)
c call AZ_broadcast(NULL,0,proc_config,AZ_SEND)
c
call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
nz = n*n*n
nzz = nz/NPROCS + 1
c Резервирование памяти под локальные массивы в зависимости от
с размера решаемой задачи и числа процессоров.
allocate( data_org(0:nzz) )
allocate( update(0:nzz) )
allocate( external(0:nzz) )
allocate( update_index(0:nzz) )
allocate( extern_index(0:nzz) )
allocate( bindx(0:nzz*7) )
allocate( b(0:nzz) )
allocate( x(0:nz) )
allocate( tmp(0:nzz) )
allocate( val(0:nzz*7) )
total = 0.0
с Цикл по методам решения системы уравнений. Метод GMRES(номер 1)
с пропускается из-за очевидных проблем с его реализацией
DO 777 K = 0,4
IF (K .EQ.1) GOTO 777
time = MPI_Wtime()
с Определение параметров разбиения системы уравнений по
с процессорам
call AZ_read_update(N_update,update,proc_config,nz,1,0)
с Формирование распределенной матрицы системы уравннений.
с Каждый процессор заполняет только свою часть матрицы
bindx(0) = N_update+1
do 250 i = 0, N_update-1
call add_row_7pt(update(i),i,val,bindx)
250 continue
с Формирование правой части системы уравнений по алгоритму:
c b(i) = sum(A(i,j)*j), в следствии чего решение будет иметь вид:
с x(i) = i ;
do 341 i = 0, N_update-1
in = bindx(i)
ik = bindx(i+1) - 1
s = val(i)*(update(i)+1)
do 342 j = in,ik
s = s + val(j)*(bindx(j)+1)
342 continue
tmp(i) = s
341 continue
c Трансформация глобальной индексации массивов в локальную
call AZ_transform(proc_config,external,bindx,val,update,
$ update_index,extern_index,data_org,
$ N_update,NULL,NULL,NULL,NULL,AZ_MSR_MATRIX)
с Инициализация параметров решателя
call AZ_defaults(options, params)
options(AZ_solver) = K
с options(AZ_precond) = AZ_dom_decomp
с options(AZ_subdomain_solve) = AZ_ilut
options(4) = 1
options(AZ_output) = AZ_none
options(AZ_max_iter) = 2000
params(AZ_tol) = 1.d-8
if (IAM.eq.0) then
write(6,780)
780 format (16x,11(2H**))
write (6,*) ' Method solve is =',met(k)
write(6,780)
write (6,*) ' Dimension matrices =',nz,
+' NPROCS = ', NPROCS
endif
с Задание начального приближения и переиндексация массивов
do 350 i = 0, N_update-1
x(update_index(i)) = 0.0
b(update_index(i)) = tmp(i)
350 continue
с Решение системы уравнений
call AZ_solve(x,b, options, params, NULL,bindx,NULL,NULL,
$ NULL,val, data_org, status, proc_config)
c Обратная переиндексация массивов
DO 415 I = 0,N_update - 1
415 tmp(i) = x(update_index(i))
c Сборка частей решения со всех процессоров
CALL MPI_ALLgather(N_update,1,MPI_INTEGER,recb,1,
*MPI_INTEGER,MPI_COMM_WORLD,ierror)
CALL MPI_BARRIER(MPI_COMM_WORLD,ierror)
DISP(0) = 0
DO 404 I = 1,NPROCS-1
404 DISP(I) = DISP(I-1) + RECB(I-1)
CALL MPI_ALLgatherv(tmp,N_update,MPI_DOUBLE_PRECISION,x,recb,
*disp,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,ierror)
CALL MPI_BARRIER(MPI_COMM_WORLD,ierror)
c Печать контрольных компонент решения
nz1 = nz-1
if (IAM.eq.0) then
do 348 i = 0,nz1,nz1
s = x(i)
write (6,956) i,s
348 continue
956 format (2x,'i=',I12,2x,'x(i)=',F18.8)
endif
time = MPI_Wtime() - time
total = total + time
if(iam.eq.0) write(*,435) time
435 FORMAT(2x,'TIME CALCULATION (SEC.) = ',F12.4)
CALL FLUSH(6)
777 continue
if(iam.eq.0) write(*,436) total
436 FORMAT(/,2x,' TOTAL TIME CALCULATION (SEC.) = ',F12.4)
CALL MPI_FINALIZE(IERROR)
stop
end
C****************************************************************
C****************************************************************
с Подпрограмма добавления очередной строки матрицы
subroutine add_row_7pt(row,location,val,bindx)
integer row,location,bindx(0:*)
double precision val(0:*)
integer n,n3
common /global/ n
C
C Параметры:
C row == глобальный номер строки, добавляемый к матрице
C location == локальный номер строки
C val,bindx == массивы для хранения ненулевых матричных
C элементов и номеров их столбцов
C
integer k
с проверка наличия соседей в расчетной сетке и добавление
с ненулевых элементов, если таковые существуют
n3 = n*n*n - 1
k = bindx(location)
bindx(k) = row + 1
if (bindx(k).le.n3) then
val(k) = -1.
k = k + 1
endif
bindx(k) = row - 1
if (bindx(k) .ge. 0) then
val(k) = -1.
k = k + 1
endif
bindx(k) = row + n
if (bindx(k).le.n3) then
val(k) = -1.
k = k + 1
endif
bindx(k) = row - n
if (bindx(k).ge. 0) then
val(k) = -1.
k = k + 1
endif
bindx(k) = row + n*n
if (bindx(k) .le. n3) then
val(k) = -1.
k = k + 1
endif
bindx(k) = row - n*n
if (bindx(k) .ge. 0) then
val(k) = -1.
k = k + 1
endif
bindx(location+1) = k
val(location) = 6.0
return
end
Для компиляции программы использовался Makefile:
PROG = azt
OBJ = $(PROG).o
FC = mpif90
CC = mpicc
FFLAGS = -O
LDFLAGS = -O
LIBS = -laztec -lmkl -lguide
all: exe
exe: $(OBJ)
${FC} ${LDFLAGS} -o $(PROG) ${OBJ} ${LIBS}
clean :
rm -f *.o
.f.o : ; $(FC) -c ${FFLAGS} $*.f
.c.o : ; $(CC) -c $(CCFLAGS) $(CDEFS) $*.c
Следует обратить внимание на то, что при компиляции программы помимо библиотеки Aztec необходимо подключать оптимизированную библиотеку blas. В данном случае используется библиотека Intel MKL. Текстовая версия библиотеки blas включена в дистрибутив библиотеки Aztec и может быть собрана вместе с библиотекой aztec, однако производительность при этом будет значительно ниже, чем при использовании бинарной оптимизированной библиотеки blas.