*******************************************************************************
* Sample driver for AZTEC package. The software is tested by setting up a *
* system of equations and a right hand side and then solving the system of *
* equations using AZTECs iterative solvers. *
* *
* 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. *
* n must be <= 50. *
*******************************************************************************
program main
include 'mpif.h'
include 'az_aztecf.h'
integer n
common /global/ n
C
C
integer i,j,IAM,NPROCS,ierr,nz,K,in,ik,nz1,nzz,no
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)
C
met(0) = 'cg'
met(1) = 'gmres'
met(2) = 'cgs'
met(3) = 'tfqmr'
met(4) = 'bicgstab'
C
C get number of processors and the name of this processor
C
call AZ_processor_info(proc_config)
IAM = proc_config(AZ_node)
NPROCS = proc_config (AZ_N_procs)
C
if (IAM.eq.0) then
write(6,1)
1 format(' '//
$ ' Victor N. Datsyuk'/
$ ' High Performance Computers Department'/
$ ' Rostov State University'/
$ ' Rostov-on-Don, phone 28-54-66'//
$ ' 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.'/
$ ' n must be <= 50.'//)
open(10, FILE='size.dat', FORM='FORMATTED')
read(10,*) n
endif
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
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
DO 777 K = 0,4
if(k.eq.1) goto 777
time = MPI_Wtime()
C
C Define paritioning:matrix rows in ascending order assigned
C to this node
C
call AZ_read_update(N_update,update,proc_config,nz,1,0)
C
C create the matrix: each processor creates only rows
C appearing in update[] (using global col. numbers).
C
bindx(0) = N_update+1
do 250 i = 0, N_update-1
call add_row_7pt(update(i),i,val,bindx)
250 continue
c
c create the right hand side of the system under rule:
c b(i) = sum(A(i,j)*j), asuming that x(i) = i ;
c
do 341 i = 0, N_update-1
in = bindx(i)
ik = bindx(i+1) - 1
s = val(i)*(update(i)+1)
c s = val(i)
do 342 j = in,ik
s = s + val(j)*(bindx(j)+1)
c s = s + val(j)
342 continue
tmp(i) = s
341 continue
c
C convert matrix to a local distributed matrix */
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)
C
C initialize AZTEC options
C
call AZ_defaults(options, params)
options(AZ_solver) = K
c options(AZ_precond) = AZ_dom_decomp
c options(AZ_subdomain_solve) = AZ_ilut
options(4) = 1
options(AZ_output) = AZ_none
options(AZ_max_iter) = 2000
params(AZ_tol) = 1.d-12
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
C
C Set rhs (delta function at grid center) and initialize guess
C
do 350 i = 0, N_update-1
x(update_index(i)) = 0.0
b(update_index(i)) = tmp(i)
350 continue
C
C solve the system of equations using b as the right hand side
C
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
c gather all parts of solve on all processors
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)
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*********************************************************************
C
subroutine add_row_7pt(row,location,val,bindx)
C
integer row,location,bindx(0:*)
double precision val(0:*)
integer n,n3
common /global/ n
C
C
C Parameters:
C row == global row number of the new row to be added.
C location == local row where diagonal of the new row will be stored.
C val,bindx == (see user's guide). On output, val[] and bindx[]
C are appended such that the new row has been added.
C
integer k
C
C check neighbors in each direction and add nonzero if neighbor exits
C
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