|
|
|
@ -97,6 +97,7 @@ program spde
|
|
|
|
|
|
|
|
|
|
! solver parameters
|
|
|
|
|
integer :: iter, itmax,itrace, istopc, irst, nlv
|
|
|
|
|
integer(psb_long_int_k_) :: amatsize, precsize, descsize
|
|
|
|
|
real(psb_spk_) :: err, eps
|
|
|
|
|
|
|
|
|
|
type precdata
|
|
|
|
@ -152,7 +153,7 @@ program spde
|
|
|
|
|
|
|
|
|
|
call psb_barrier(ictxt)
|
|
|
|
|
t1 = psb_wtime()
|
|
|
|
|
call create_matrix(idim,a,b,x,desc_a,part_block,ictxt,afmt,info)
|
|
|
|
|
call create_matrix(idim,a,b,x,desc_a,ictxt,afmt,info)
|
|
|
|
|
call psb_barrier(ictxt)
|
|
|
|
|
t2 = psb_wtime() - t1
|
|
|
|
|
if(info /= 0) then
|
|
|
|
@ -233,6 +234,12 @@ program spde
|
|
|
|
|
t2 = psb_wtime() - t1
|
|
|
|
|
call psb_amx(ictxt,t2)
|
|
|
|
|
|
|
|
|
|
amatsize = psb_sizeof(a)
|
|
|
|
|
descsize = psb_sizeof(desc_a)
|
|
|
|
|
precsize = mld_sizeof(prec)
|
|
|
|
|
call psb_sum(ictxt,amatsize)
|
|
|
|
|
call psb_sum(ictxt,descsize)
|
|
|
|
|
call psb_sum(ictxt,precsize)
|
|
|
|
|
if (iam == psb_root_) then
|
|
|
|
|
write(*,'(" ")')
|
|
|
|
|
write(*,'("Time to solve matrix : ",es12.5)')t2
|
|
|
|
@ -240,6 +247,9 @@ program spde
|
|
|
|
|
write(*,'("Number of iterations : ",i0)')iter
|
|
|
|
|
write(*,'("Convergence indicator on exit : ",es12.5)')err
|
|
|
|
|
write(*,'("Info on exit : ",i0)')info
|
|
|
|
|
write(*,'("Total memory occupation for A: ",i12)')amatsize
|
|
|
|
|
write(*,'("Total memory occupation for DESC_A: ",i12)')descsize
|
|
|
|
|
write(*,'("Total memory occupation for PREC: ",i12)')precsize
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
@ -382,7 +392,7 @@ contains
|
|
|
|
|
! subroutine to allocate and fill in the coefficient matrix and
|
|
|
|
|
! the rhs.
|
|
|
|
|
!
|
|
|
|
|
subroutine create_matrix(idim,a,b,xv,desc_a,parts,ictxt,afmt,info)
|
|
|
|
|
subroutine create_matrix(idim,a,b,xv,desc_a,ictxt,afmt,info)
|
|
|
|
|
!
|
|
|
|
|
! discretize the partial diferential equation
|
|
|
|
|
!
|
|
|
|
@ -407,20 +417,11 @@ contains
|
|
|
|
|
type(psb_desc_type) :: desc_a
|
|
|
|
|
integer :: ictxt, info
|
|
|
|
|
character :: afmt*5
|
|
|
|
|
interface
|
|
|
|
|
! .....user passed subroutine.....
|
|
|
|
|
subroutine parts(global_indx,n,np,pv,nv)
|
|
|
|
|
implicit none
|
|
|
|
|
integer, intent(in) :: global_indx, n, np
|
|
|
|
|
integer, intent(out) :: nv
|
|
|
|
|
integer, intent(out) :: pv(*)
|
|
|
|
|
end subroutine parts
|
|
|
|
|
end interface ! local variables
|
|
|
|
|
type(psb_sspmat_type) :: a
|
|
|
|
|
real(psb_spk_) :: zt(nb),glob_x,glob_y,glob_z
|
|
|
|
|
integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k
|
|
|
|
|
integer :: x,y,z,ia,indx_owner
|
|
|
|
|
integer :: np, iam
|
|
|
|
|
integer :: np, iam, nr, nt
|
|
|
|
|
integer :: element
|
|
|
|
|
integer, allocatable :: irow(:),icol(:),myidx(:)
|
|
|
|
|
real(psb_spk_), allocatable :: val(:)
|
|
|
|
@ -432,7 +433,6 @@ contains
|
|
|
|
|
real(psb_spk_) :: a1, a2, a3, a4, b1, b2, b3
|
|
|
|
|
external :: a1, a2, a3, a4, b1, b2, b3
|
|
|
|
|
integer :: err_act
|
|
|
|
|
! common area
|
|
|
|
|
|
|
|
|
|
character(len=20) :: name, ch_err
|
|
|
|
|
|
|
|
|
@ -452,9 +452,18 @@ contains
|
|
|
|
|
nnz = ((n*9)/(np))
|
|
|
|
|
if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0x,")...")')n
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Using a simple BLOCK distribution.
|
|
|
|
|
!
|
|
|
|
|
nt = (m+np-1)/np
|
|
|
|
|
nr = min(nt,m-(iam*nt))
|
|
|
|
|
|
|
|
|
|
nt = nr
|
|
|
|
|
call psb_sum(ictxt,nt)
|
|
|
|
|
if (nt /= m) write(0,*) iam, 'Initialization error ',nr,nt,m
|
|
|
|
|
call psb_barrier(ictxt)
|
|
|
|
|
t0 = psb_wtime()
|
|
|
|
|
call psb_cdall(ictxt,desc_a,info,mg=n,parts=parts)
|
|
|
|
|
call psb_cdall(ictxt,desc_a,info,nl=nr)
|
|
|
|
|
if (info == 0) call psb_spall(a,desc_a,info,nnz=nnz)
|
|
|
|
|
! define rhs from boundary conditions; also build initial guess
|
|
|
|
|
if (info == 0) call psb_geall(b,desc_a,info)
|
|
|
|
@ -607,9 +616,9 @@ contains
|
|
|
|
|
val(element)=-b2(glob_x,glob_y,glob_z)
|
|
|
|
|
val(element) = val(element)/(deltah*&
|
|
|
|
|
& deltah)
|
|
|
|
|
icol(element)=(x-1)*idim*idim+(y)*idim+(z)
|
|
|
|
|
icol(element) = (x-1)*idim*idim+(y)*idim+(z)
|
|
|
|
|
irow(element) = glob_row
|
|
|
|
|
element=element+1
|
|
|
|
|
element = element+1
|
|
|
|
|
endif
|
|
|
|
|
! term depending on (x+1,y,z)
|
|
|
|
|
if (x<idim) then
|
|
|
|
|