c=======================================================================
c ROUTINE: mgv
c
c Solves Poisson equation
c
c Laplacian[u] = rhs
c
c in 1-, 2- and 3-D using fixed V-cycle LCS
c (linear correction scheme) MG.
c
c NOTE: This routine is a real*8 FUNCTION. Be sure to declare
c it as such in the invoking program unit.
c
c TYPICAL USAGE: (I) and (O) denote input and ouptut arguments
c respectively
c
c real*8 mgv
c
c integer rank
c parameter ( rank = 3) ! Problem dimension (I)
c integer nx, ny, nz
c parameter ( nx = 33, ny = 33, nz = 33 )
c real*8 u(nx,ny,nz) ! Soln (O)
c real*8 rhs(nx,ny,nz) ! Rhs (I)
c integer sizeq
c parameter ( sizeq = 3 * nx * ny * nz ) ! Size of work array (I)
c real*8 q(sizeq) ! Work array (I)
c real*8 vh(rank) ! Mesh spacings (I)
c integer shape(rank) ! Mesh dimensions (I)
c integer preswp, pstswp, ncycle ! Multi-grid parameters (I)
c
c
c vh(1) = h ! This example assumes mesh spacings in
c ! all coordinate direction = h
c vh(2) = h
c vh(3) = h
c shape(1) = nx
c shape(2) = ny
c shape(3) = nz
c
c preswp = 3 ! # of pre-CGC sweeps
c pstswp = 3 ! # of post-CGC sweeps
c ncycle = 1 ! # of V-cycles
c
c ! Define rhs array
c .
c .
c .
c
c ! mgv returns norm of residuals after ncycle (preswp,pstswp)
c ! v-cycles.
c
c resid = mgv(rank,shape,vh,u,rhs,q,sizeq,preswp,pstswp,ncycle)
c .
c .
c .
c
c
c IMPLEMENTATION NOTES: This routine is a front end to mg1v0,
c mg2v0 and mg3v0. It invokes those routines with
c
c preswf = presw = preswp
c pstswf = pstsw = pstswp
c=======================================================================
double precision function mgv
* (rank,shape,vh,uf,rhsf,
* q,sizeq,
* preswp,pstswp,ncycle)
implicit none
real*8 mg1v0, mg2v0, mg3v0
integer rank, sizeq
integer shape(rank)
real*8 vh(rank), uf(*), rhsf(*),
* q(sizeq)
integer preswp, pstswp, ncycle
logical ltrace
parameter ( ltrace = .false. )
if( ltrace ) then
write(0,*) 'mgv: Parameter dump'
write(0,*) 'rank, preswp, pstswp, ncycle'
write(0,*) rank, preswp, pstswp, ncycle
call ivdump(shape,rank,'h',0)
call dvdump(vh,rank,'h',0)
end if
if( rank .eq. 1 ) then
mgv = mg1v0(uf,rhsf,shape(1), vh(1),
* q,sizeq,
* preswp,pstswp,preswp,pstswp,ncycle)
else if( rank .eq. 2 ) then
mgv = mg2v0(uf,rhsf,shape(1),shape(2), vh(1),
* q,sizeq,
* preswp,pstswp,preswp,pstswp,ncycle)
else if( rank .eq. 3 ) then
mgv = mg3v0(uf,rhsf,shape(1),shape(2),shape(3),vh(1),
* q,sizeq,
* preswp,pstswp,preswp,pstswp,ncycle)
else
write(0,*) 'mgv: Invalid rank: ', rank
mgv = -1.0d0
end if
return
end
c-----------------------------------------------------------------------
c
c Front end to mg1v0 which invokes mg1v0 with
c
c preswf = presw = preswp
c pstswf = pstsw = pstswp
c
c-----------------------------------------------------------------------
double precision function mg1v
* (uf,rhsf,nxf, hf, q,sizeq,
* preswp,pstswp,ncycle)
implicit none
real*8 mg1v0
integer nxf, sizeq
real*8 uf(nxf), rhsf(nxf),
* q(sizeq)
real*8 hf
integer preswp, pstswp, ncycle
mg1v = mg1v0(uf,rhsf,nxf, hf,q,sizeq,
* preswp,pstswp,preswp,pstswp,ncycle)
return
end
c-----------------------------------------------------------------------
c
c Front end to mg2v0 which invokes mg2v0 with
c
c preswf = presw = preswp
c pstswf = pstsw = pstswp
c
c-----------------------------------------------------------------------
double precision function mg2v
* (uf,rhsf,nxf,nyf, hf, q,sizeq,
* preswp,pstswp,ncycle)
implicit none
real*8 mg2v0
integer nxf, nyf, sizeq
real*8 uf(nxf,nyf), rhsf(nxf,nyf),
* q(sizeq)
real*8 hf
integer preswp, pstswp, ncycle
mg2v = mg2v0(uf,rhsf,nxf,nyf,hf,q,sizeq,
* preswp,pstswp,preswp,pstswp,ncycle)
return
end
c-----------------------------------------------------------------------
c
c Front end to mg3v0 which invokes mg3v0 with
c
c preswf = presw = preswp
c pstswf = pstsw = pstswp
c
c-----------------------------------------------------------------------
double precision function mg3v
* (uf,rhsf,nxf,nyf,nzf, hf, q,sizeq,
* preswp,pstswp,ncycle)
implicit none
real*8 mg3v0
integer nxf, nyf, nzf, sizeq
real*8 uf(nxf,nyf,nzf), rhsf(nxf,nyf,nzf),
* q(sizeq)
real*8 hf
integer preswp, pstswp, ncycle
mg3v = mg3v0(uf,rhsf,nxf,nyf,nzf,hf,q,sizeq,
* preswp,pstswp,preswp,pstswp,ncycle)
return
end
c-----------------------------------------------------------------------
c
c Basic 1d v-cycle multi grid solver for
c
c Laplacian[u(x)] = rhs(x)
c
c This version allows nxf, nyf which are of form
c 2^p * <small number> + 1
c
c where <small number> .le. maxn0 defined in the parameter
c statement below
c
c-----------------------------------------------------------------------
double precision function mg1v0
* (uf,rhsf,nxf, hf, q, sizeq,
* preswf,pstswf,presw,pstsw,ncycle)
implicit none
real*8 relaxc1
integer qalloc, lilog2
logical lp2np1
integer nxf, sizeq
real*8 uf(nxf), rhsf(nxf),
* q(sizeq)
real*8 hf
integer preswf, pstswf, presw, pstsw,
* ncycle
c
c Multi-grid parameters, data structures and pointers
c
real*8 epsic
parameter ( epsic = 1.0d-10 )
integer lmin, lmax
parameter ( lmin = 1, lmax = 20 )
integer nx(lmin:lmax)
real*8 h(lmin:lmax)
integer u(lmin:lmax), rhs(lmin:lmax),
* mres(lmin:lmax), du(lmin:lmax)
integer imposs
parameter ( imposs = -2 000 000 000 )
integer p, pceil
common / comp / p, pceil
c
c For parameter (nxf) checking
c
integer checkfacts
logical badargs
integer nfact_x, nx0
integer maxn0
parameter ( maxn0 = 17 )
c
c Locals ...
c
integer lxf,
* l, lf, lc, nl,
* sweep
real*8 resnrm
logical ltrace
parameter ( ltrace = .false. )
save
mg1v0 = -1.0d0
c
c Parameter check ...
c
if( ltrace ) then
write(0,*) '>>> mg1v0: Parameter dump: '
write(0,*) nxf, hf
write(0,*) preswf, pstswf, presw, pstsw, ncycle
end if
badargs = .false.
if( checkfacts(nxf - 1,maxn0,nx0,nfact_x) .lt. 0 ) then
write(*,*) 'mg1v0: Bad nxf value: ', nxf
badargs = .true.
end if
if( badargs) then
return
end if
nl = nfact_x + 1
lc = 1
lf = lc + nl - 1
nx0 = (nxf - 1) / (2 ** (nl - 1)) + 1
if( ltrace ) then
write(*,*) 'mg1v0: ', nl, ' levels.'
write(*,*) 'mg1v0: Coarse grid ', nx0
end if
if( ncycle .eq. 0 ) then
c
c Set up grid structure, allocate storage ...
c
call qreset(sizeq)
nx(lf) = nxf
h(lf) = hf
u(lf) = imposs
rhs(lf) = imposs
mres(lf) = qalloc(nx(lf))
du(lf) = mres(lf)
do 1000 l = lf - 1 , lc , -1
nx(l) = (nx(l+1) - 1) / 2 + 1
h(l) = 2.0d0 * h(l+1)
u(l) = qalloc(nx(l))
rhs(l) = qalloc(nx(l))
mres(l) = qalloc(nx(l))
du(l) = mres(l)
1000 continue
if( ltrace ) then
call ivdump(nx,nl,'nx',0)
call dvdump(h,nl,'h',0)
call ivdump(u,nl,'u pointers',0)
call ivdump(rhs,nl,'rhs pointers',0)
call ivdump(mres,nl,'mres pointers',0)
call ivdump(du,nl,'du pointers',0)
end if
end if
c
c Cycle from fine to coarse ...
c
do l = lf , lc , -1
if( l .ne. lf .or. ncycle .eq. 0 ) then
if( l .eq. lc ) then
100 continue
resnrm = relaxc1(q(u(l)),q(rhs(l)),
* nx(l),h(l))
if( ltrace ) then
call trace_res('mg1v0: pre sweep',l,resnrm)
call cres1_trace('mg1v0: pre lc',60,
& q(mres(l)),q(u(l)),q(rhs(l)),
& nx(l),h(l))
end if
if( resnrm .gt. epsic ) go to 100
else if( l .eq. lf ) then
do sweep = 1 , preswf
resnrm = relaxc1(uf,rhsf,
* nx(l),h(l))
call cmres1(q(mres(l)),uf,rhsf,
* nx(l),h(l))
if( ltrace ) then
call trace_res('mg1v0: pre sweep',l,resnrm)
call cres1_trace('mg1v0: pre lf',60,
& q(mres(l)),uf,rhsf,
& nx(l),h(l))
end if
end do
else
do sweep = 1 , presw
resnrm = relaxc1(q(u(l)),q(rhs(l)),
* nx(l),h(l))
if( ltrace ) then
call trace_res('mg1v0: pre sweep',l,resnrm)
call cres1_trace('mg1v0: pre inter',60,
& q(mres(l)),q(u(l)),q(rhs(l)),
& nx(l),h(l))
end if
end do
end if
end if
c
c Set up coarse grid problem ...
c
if( l .eq. lf ) then
call cmres1(q(mres(l)),uf,rhsf,nx(l),h(l))
else if( l .ne. lc ) then
call cmres1(q(mres(l)),q(u(l)),q(rhs(l)),
* nx(l),h(l))
end if
if( l .ne. lc ) then
call dvrshw(q(mres(l)),q(rhs(l-1)),
* nx(l-1))
call dvls(q(u(l-1)),0.0d0,
* nx(l-1))
end if
end do
c
c Cycle from coarse to fine ...
c
do l = lc + 1 , lf
c
c Interpolate correction ...
c
call dvlint(q(u(l-1)),q(du(l)),nx(l-1))
if( l .eq. lf ) then
call dvva(uf,q(du(l)),uf,nx(l))
else
call dvva(q(u(l)),q(du(l)),q(u(l)),nx(l))
end if
if( l .eq. lf ) then
do sweep = 1 , pstswf
resnrm = relaxc1(uf,rhsf,nx(l),h(l))
if( ltrace ) then
call trace_res('mg1v0: post sweep',l,resnrm)
call cres1_trace('mg1v0: post lf',60,
& q(mres(l)),uf,rhsf,
& nx(l),h(l))
end if
end do
else
do sweep = 1 , pstsw
resnrm = relaxc1(q(u(l)),q(rhs(l)),
* nx(l),h(l))
if( ltrace ) then
call trace_res('mg1v0: post sweep',l,resnrm)
call cres1_trace('mg1v0: post other',60,
& q(mres(l)),q(u(l)),q(rhs(l)),
& nx(l),h(l))
end if
end do
end if
end do
ncycle = ncycle + 1
mg1v0 = resnrm
return
end
c-----------------------------------------------------------------------
c
c Basic 2d v-cycle multi grid solver for
c
c Laplacian[u(x,y)] = rhs(x,y)
c
c on rectangular domain.
c
c This version allows nxf, nyf which are of form
c 2^p * <small number> + 1
c
c where <small number> .le. maxn0 defined in the parameter
c statement below
c
c-----------------------------------------------------------------------
double precision function mg2v0
* (uf,rhsf,nxf,nyf, hf, q, sizeq,
* preswf,pstswf,presw,pstsw,ncycle)
implicit none
real*8 relaxc2
integer qalloc, lilog2
logical lp2np1
integer nxf, nyf, sizeq
real*8 uf(nxf,nyf), rhsf(nxf,nyf),
* q(sizeq)
real*8 hf
integer preswf, pstswf, presw, pstsw,
* ncycle
c
c Multi-grid parameters, data structures and pointers
c
real*8 epsic
parameter ( epsic = 1.0d-10 )
integer lmin, lmax
parameter ( lmin = 1, lmax = 20 )
integer nx(lmin:lmax), ny(lmin:lmax)
real*8 h(lmin:lmax)
integer u(lmin:lmax), rhs(lmin:lmax),
* mres(lmin:lmax), du(lmin:lmax)
integer imposs
parameter ( imposs = -2 000 000 000 )
integer p, pceil
common / comp / p, pceil
c
c For parameter (nxf, nyf) checking
c
integer checkfacts
logical badargs
integer nfact_x, nx0
integer nfact_y, ny0
integer maxn0
parameter ( maxn0 = 17 )
c
c Locals ...
c
integer lxf, lyf,
* l, lf, lc, nl,
* sweep
real*8 resnrm
logical ltrace
parameter ( ltrace = .false. )
save
mg2v0 = -1.0d0
c
c Parameter check ...
c
if( ltrace ) then
write(0,*) '>>> mg2v0: Parameter dump: '
write(0,*) nxf, nyf, hf
write(0,*) preswf, pstswf, presw, pstsw, ncycle
end if
badargs = .false.
if( checkfacts(nxf - 1,maxn0,nx0,nfact_x) .lt. 0 ) then
write(*,*) 'mg2v0: Bad nxf value: ', nxf
badargs = .true.
end if
if( checkfacts(nyf - 1,maxn0,ny0,nfact_y) .lt. 0 ) then
write(*,*) 'mg2v0: Bad nyf value: ', nyf
badargs = .true.
end if
if( badargs) then
return
end if
nl = min(nfact_x,nfact_y) + 1
lc = 1
lf = lc + nl - 1
nx0 = (nxf - 1) / (2 ** (nl - 1)) + 1
ny0 = (nyf - 1) / (2 ** (nl - 1)) + 1
if( ltrace ) then
write(*,*) 'mg2v0: ', nl, ' levels.'
write(*,*) 'mg2v0: Coarse grid ', nx0, ' x ', ny0
end if
if( ncycle .eq. 0 ) then
c
c Set up grid structure, allocate storage ...
c
call qreset(sizeq)
nx(lf) = nxf
ny(lf) = nyf
h(lf) = hf
u(lf) = imposs
rhs(lf) = imposs
mres(lf) = qalloc(nx(lf)*ny(lf))
du(lf) = mres(lf)
do 1000 l = lf - 1 , lc , -1
nx(l) = (nx(l+1) - 1) / 2 + 1
ny(l) = (ny(l+1) - 1) / 2 + 1
h(l) = 2.0d0 * h(l+1)
u(l) = qalloc(nx(l)*ny(l))
rhs(l) = qalloc(nx(l)*ny(l))
mres(l) = qalloc(nx(l)*ny(l))
du(l) = mres(l)
1000 continue
if( ltrace ) then
call ivdump(nx,nl,'nx',0)
call ivdump(ny,nl,'ny',0)
call dvdump(h,nl,'h',0)
call ivdump(u,nl,'u pointers',0)
call ivdump(rhs,nl,'rhs pointers',0)
call ivdump(mres,nl,'mres pointers',0)
call ivdump(du,nl,'du pointers',0)
end if
end if
c
c Cycle from fine to coarse ...
c
do l = lf , lc , -1
if( l .ne. lf .or. ncycle .eq. 0 ) then
if( l .eq. lc ) then
100 continue
resnrm = relaxc2(q(u(l)),q(rhs(l)),
* nx(l),ny(l),h(l))
if( ltrace ) then
call trace_res('mg2v0: pre sweep',l,resnrm)
end if
if( resnrm .gt. epsic ) go to 100
else if( l .eq. lf ) then
do sweep = 1 , preswf
resnrm = relaxc2(uf,rhsf,
* nx(l),ny(l),h(l))
call cmres2(q(mres(l)),uf,rhsf,
* nx(l),ny(l),h(l))
if( ltrace ) then
call trace_res('mg2v0: pre sweep',l,resnrm)
end if
end do
else
do sweep = 1 , presw
resnrm = relaxc2(q(u(l)),q(rhs(l)),
* nx(l),ny(l),h(l))
if( ltrace ) then
call trace_res('mg2v0: pre sweep',l,resnrm)
end if
end do
end if
end if
c
c Set up coarse grid problem ...
c
if( l .eq. lf ) then
call cmres2(q(mres(l)),uf,rhsf,nx(l),ny(l),h(l))
else if( l .ne. lc ) then
call cmres2(q(mres(l)),q(u(l)),q(rhs(l)),
* nx(l),ny(l),h(l))
end if
if( l .ne. lc ) then
call dmrshw(q(mres(l)),q(rhs(l-1)),
* nx(l-1),ny(l-1))
call dmls(q(u(l-1)),0.0d0,
* nx(l-1),ny(l-1))
end if
end do
c
c Cycle from coarse to fine ...
c
do l = lc + 1 , lf
c
c Interpolate correction ...
c
call dmlint(q(u(l-1)),q(du(l)),nx(l-1),ny(l-1))
if( l .eq. lf ) then
call dmma(uf,q(du(l)),uf,nx(l),ny(l))
else
call dmma(q(u(l)),q(du(l)),q(u(l)),nx(l),ny(l))
end if
if( l .eq. lf ) then
do sweep = 1 , pstswf
resnrm = relaxc2(uf,rhsf,nx(l),ny(l),h(l))
if( ltrace ) then
call trace_res('mg2v0: post sweep',l,resnrm)
end if
end do
else
do sweep = 1 , pstsw
resnrm = relaxc2(q(u(l)),q(rhs(l)),
* nx(l),ny(l),h(l))
if( ltrace ) then
call trace_res('mg2v0: post sweep',l,resnrm)
end if
end do
end if
end do
ncycle = ncycle + 1
mg2v0 = resnrm
return
end
c-----------------------------------------------------------------------
c
c Basic 3d v-cycle multi grid solver for
c
c Laplacian[u(x,y,z)] = rhs(x,y,z)
c
c on rectangular domain.
c
c This version allows nxf, nyf, nzy which are of form
c 2^p * <small number> + 1
c
c where <small number> .le. maxn0 defined in the parameter
c statement below
c
c-----------------------------------------------------------------------
double precision function mg3v0
* (uf,rhsf,nxf,nyf,nzf, hf, q, sizeq,
* preswf,pstswf,presw,pstsw,ncycle)
implicit none
real*8 d3nrm2, relaxc3
integer qalloc, lilog2
logical lp2np1
integer nxf, nyf, nzf, sizeq
real*8 uf(nxf,nyf,nzf), rhsf(nxf,nyf,nzf),
* q(sizeq)
real*8 hf
integer preswf, pstswf, presw, pstsw,
* ncycle
c
c Multi-grid parameters, data structures and pointers
c
real*8 epsic
parameter ( epsic = 1.0d-10 )
integer lmin, lmax
parameter ( lmin = 1, lmax = 20 )
integer nx(lmin:lmax), ny(lmin:lmax),
* nz(lmin:lmax)
real*8 h(lmin:lmax)
integer u(lmin:lmax), rhs(lmin:lmax),
* mres(lmin:lmax), du(lmin:lmax)
integer imposs
parameter ( imposs = -2 000 000 000 )
integer p, pceil
common / comp / p, pceil
c
c For parameter (nxf, nyf, nzf) checking
c
integer checkfacts
logical badargs
integer nfact_x, nx0
integer nfact_y, ny0
integer nfact_z, nz0
integer maxn0
parameter ( maxn0 = 17 )
c
c Locals ...
c
integer lxf, lyf, lzf,
* l, lf, lc, nl,
* sweep
real*8 resnrm
logical ltrace
parameter ( ltrace = .false. )
save
mg3v0 = -1.0d0
c
c Parameter check ...
c
if( ltrace ) then
write(0,*) '>>> mg3v0: Parameter dump: '
write(0,*) nxf, nyf, nzf, hf
write(0,*) preswf, pstswf, presw, pstsw, ncycle
end if
badargs = .false.
if( checkfacts(nxf - 1,maxn0,nx0,nfact_x) .lt. 0 ) then
write(*,*) 'mg3v0: Bad nxf value: ', nxf
badargs = .true.
end if
if( checkfacts(nyf - 1,maxn0,ny0,nfact_y) .lt. 0 ) then
write(*,*) 'mg3v0: Bad nyf value: ', nyf
badargs = .true.
end if
if( checkfacts(nzf - 1,maxn0,nz0,nfact_z) .lt. 0 ) then
write(*,*) 'mg3v0: Bad nzf value: ', nzf
badargs = .true.
end if
if( badargs) then
return
end if
nl = min(nfact_x,nfact_y,nfact_z) + 1
lc = 1
lf = lc + nl - 1
nx0 = (nxf - 1) / (2 ** (nl - 1)) + 1
ny0 = (nyf - 1) / (2 ** (nl - 1)) + 1
nz0 = (nzf - 1) / (2 ** (nl - 1)) + 1
if( ltrace ) then
write(*,*) 'mg3v0: ', nl, ' levels.'
write(*,*) 'mg3v0: Coarse grid ', nx0, ' x ', ny0,
* ' x ', nz0
end if
if( ncycle .eq. 0 ) then
c
c Set up grid structure, allocate storage ...
c
call qreset(sizeq)
nx(lf) = nxf
ny(lf) = nyf
nz(lf) = nzf
h(lf) = hf
u(lf) = imposs
rhs(lf) = imposs
mres(lf) = qalloc(nx(lf)*ny(lf)*nz(lf))
du(lf) = mres(lf)
do 1000 l = lf - 1 , lc , -1
nx(l) = (nx(l+1) - 1) / 2 + 1
ny(l) = (ny(l+1) - 1) / 2 + 1
nz(l) = (nz(l+1) - 1) / 2 + 1
h(l) = 2.0d0 * h(l+1)
u(l) = qalloc(nx(l)*ny(l)*nz(l))
rhs(l) = qalloc(nx(l)*ny(l)*nz(l))
mres(l) = qalloc(nx(l)*ny(l)*nz(l))
du(l) = mres(l)
1000 continue
if( ltrace ) then
call ivdump(nx,nl,'nx',0)
call ivdump(ny,nl,'ny',0)
call ivdump(nz,nl,'nz',0)
call dvdump(h,nl,'h',0)
call ivdump(u,nl,'u pointers',0)
call ivdump(rhs,nl,'rhs pointers',0)
call ivdump(mres,nl,'mres pointers',0)
call ivdump(du,nl,'du pointers',0)
end if
end if
c
c Cycle from fine to coarse ...
c
do l = lf , lc , -1
if( l .ne. lf .or. ncycle .eq. 0 ) then
if( l .eq. lc ) then
100 continue
resnrm = relaxc3(q(u(l)),q(rhs(l)),
* nx(l),ny(l),nz(l),h(l))
if( ltrace ) then
call trace_res('mg3v0: pre sweep',l,resnrm)
end if
if( resnrm .gt. epsic ) go to 100
else if( l .eq. lf ) then
do sweep = 1 , preswf
resnrm = relaxc3(uf,rhsf,
* nx(l),ny(l),nz(l),h(l))
call cmres3(q(mres(l)),uf,rhsf,
* nx(l),ny(l),nz(l),h(l))
if( ltrace ) then
call trace_res('mg3v0: pre sweep',l,resnrm)
end if
end do
else
do sweep = 1 , presw
resnrm = relaxc3(q(u(l)),q(rhs(l)),
* nx(l),ny(l),nz(l),h(l))
if( ltrace ) then
call trace_res('mg3v0: pre sweep',l,resnrm)
end if
end do
end if
end if
c
c Set up coarse grid problem ...
c
if( l .eq. lf ) then
call cmres3(q(mres(l)),uf,rhsf,nx(l),ny(l),nz(l),h(l))
else if( l .ne. lc ) then
call cmres3(q(mres(l)),q(u(l)),q(rhs(l)),
* nx(l),ny(l),nz(l),h(l))
end if
if( l .ne. lc ) then
call d3rshw(q(mres(l)),q(rhs(l-1)),
* nx(l-1),ny(l-1),nz(l-1))
call d3ls(q(u(l-1)),0.0d0,
* nx(l-1),ny(l-1),nz(l-1))
end if
end do
c
c Cycle from coarse to fine ...
c
do l = lc + 1 , lf
c
c Interpolate correction ...
c
call d3lint(q(u(l-1)),q(du(l)),nx(l-1),ny(l-1),nz(l-1))
if( l .eq. lf ) then
call d33a(uf,q(du(l)),uf,nx(l),ny(l),nz(l))
else
call d33a(q(u(l)),q(du(l)),q(u(l)),nx(l),ny(l),nz(l))
end if
if( l .eq. lf ) then
do sweep = 1 , pstswf
resnrm = relaxc3(uf,rhsf,nx(l),ny(l),nz(l),h(l))
if( ltrace ) then
call trace_res('mg3v0: post sweep',l,resnrm)
end if
end do
else
do sweep = 1 , pstsw
resnrm = relaxc3(q(u(l)),q(rhs(l)),
* nx(l),ny(l),nz(l),h(l))
if( ltrace ) then
call trace_res('mg3v0: post sweep',l,resnrm)
end if
end do
end if
end do
ncycle = ncycle + 1
mg3v0 = resnrm
return
end
c-----------------------------------------------------------------------
c
c Resets pointer into main storage array.
c
c-----------------------------------------------------------------------
subroutine qreset(ceil)
implicit none
integer ceil
integer p, pceil
common / comp / p, pceil
pceil = ceil
p = 1
return
end
c-----------------------------------------------------------------------
c
c Returns next available location in main storage array and updates
c pointer.
c
c-----------------------------------------------------------------------
integer function qalloc(size)
integer size
integer p, pceil
common / comp / p, pceil
if( size .gt. 0 ) then
if( p + size .le. pceil ) then
qalloc = p
p = p + size
else
write(0,1000) int(size / 1.0d6 + 0.5d0),
& int(pceil / 1.0d6 + 0.5d0)
1000 format(' >>> qalloc: Out of memory'/
& ' >>> qalloc: Requested block size:', i5,
& ' Mbytes'/
& ' >>> qalloc: Total available: ', i5,
& ' Mbytes')
stop
end if
else
write(0,*) '>>> qalloc: size ', size, ' <= 0.'
qalloc = -1
end if
return
end
c-----------------------------------------------------------------------
c
c Returns number of words of memory used.
c
c-----------------------------------------------------------------------
integer function qused()
integer p, pceil
common / comp / p, pceil
qused = p
return
end
c-----------------------------------------------------------------------
c
c Returns number of words of available memory.
c
c-----------------------------------------------------------------------
integer function qleft()
integer p, pceil
common / comp / p, pceil
qleft = pceil - p + 1
return
end
c-----------------------------------------------------------------------
c
c Tracing routine
c
c-----------------------------------------------------------------------
subroutine trace_res(string,l,resnrm)
implicit none
character*(*) string
integer l
real*8 resnrm
write(0,1000) string, l, resnrm
1000 format(a,' Level: ',i2,' |res.|: ',1pe12.3)
return
end
c-----------------------------------------------------------------------
c
c Does one-step 1d Point-GS relaxation and returns norm of
c running ('instantaneous') residuals---checkerboard ordering.
c
c-----------------------------------------------------------------------
double precision function relaxc1(u,rhs,nx,h)
implicit none
integer nx
real*8 u(nx), rhs(nx)
real*8 h
real*8 hm2, jacel, rresl
integer i, pass
relaxc1 = 0.0d0
hm2 = 1.0d0 / (h * h)
jacel = -2.0d0 * hm2
do pass = 1 , 2
do i = pass + 1 , nx - 1 , 2
rresl = hm2 * (u(i+1) + u(i-1) - 2.0d0*u(i)) -
* rhs(i)
u(i) = u(i) - rresl / jacel
relaxc1 = relaxc1 + rresl * rresl
end do
end do
relaxc1 = sqrt(relaxc1 / nx)
return
end
c-----------------------------------------------------------------------
c
c Computes negative of 1d residuals ...
c
c-----------------------------------------------------------------------
subroutine cmres1(mres,u,rhs,nx,h)
implicit none
integer nx
real*8 mres(nx), u(nx),
* rhs(nx)
real*8 h
integer i
real*8 hm2
hm2 = 1.0d0 / (h * h)
do i = 2 , nx - 1
mres(i) = rhs(i) -
* (hm2 * (u(i-1) - 2.0d0 * u(i) + u(i+1)))
end do
call dvlbs(mres,0.0d0,nx)
return
end
c-----------------------------------------------------------------------
c
c Does one-step 2d Point-GS relaxation and returns norm of
c running ('instantaneous') residuals---checkerboard ordering.
c
c-----------------------------------------------------------------------
double precision function relaxc2(u,rhs,nx,ny,h)
implicit none
integer nx, ny
real*8 u(nx,ny), rhs(nx,ny)
real*8 h
real*8 hm2, mjelm1, rresl
integer i, j,
* isw, jsw, pass
relaxc2 = 0.0d0
hm2 = 1.0d0 / (h * h)
mjelm1 = h * h / 4.0d0
jsw = 1
do pass = 1 , 2
isw = jsw
do j = 2 , ny - 1
do i = isw + 1 , nx - 1, 2
rresl = hm2 * (u(i+1,j) + u(i-1,j) +
* u(i,j+1) + u(i,j-1) -
* 4.0d0 * u(i,j)) -
* rhs(i,j)
u(i,j) = u(i,j) + mjelm1 * rresl
relaxc2 = relaxc2 + rresl * rresl
end do
isw = 3 - isw
end do
jsw = 3 - jsw
end do
relaxc2 = sqrt(relaxc2 / (nx * ny))
return
end
c-----------------------------------------------------------------------
c
c Computes negative of 2d residuals ...
c
c-----------------------------------------------------------------------
subroutine cmres2(mres,u,rhs,nx,ny,h)
implicit none
integer nx, ny
real*8 mres(nx,ny), u(nx,ny),
* rhs(nx,ny)
real*8 h
integer i, j
real*8 hm2
hm2 = 1.0d0 / (h * h)
do j = 2 , ny - 1
do i = 2 , nx - 1
mres(i,j) = rhs(i,j) -
* hm2 * (u(i+1,j) + u(i-1,j) +
* u(i,j+1) + u(i,j-1) -
* 4.0d0*u(i,j))
end do
end do
call dmlbs(mres,0.0d0,nx,ny)
return
end
c-----------------------------------------------------------------------
c
c Does one-step 3d Point-GS relaxation and returns norm of
c running ('instantaneous') residuals---checkerboard ordering.
c
c-----------------------------------------------------------------------
double precision function relaxc3(u,rhs,nx,ny,nz,h)
implicit none
integer nx, ny, nz
real*8 u(nx,ny,nz), rhs(nx,ny,nz)
real*8 h
real*8 hm2, mjelm1, rresl
integer i, j, k,
* isw, ksw, pass
relaxc3 = 0.0d0
hm2 = 1.0d0 / (h * h)
mjelm1 = h * h / 6.0d0
ksw = 2
do pass = 1 , 2
isw = ksw
do k = 2 , nz - 1
do j = 2 , ny - 1
do i = isw + 1 , nx - 1, 2
rresl = hm2 * (u(i+1,j,k) + u(i-1,j,k) +
* u(i,j+1,k) + u(i,j-1,k) +
* u(i,j,k+1) + u(i,j,k-1) -
* 6.0d0 * u(i,j,k)) -
* rhs(i,j,k)
u(i,j,k) = u(i,j,k) + mjelm1 * rresl
relaxc3 = relaxc3 + rresl * rresl
end do
isw = 3 - isw
end do
end do
ksw = 3 - ksw
end do
relaxc3 = sqrt(relaxc3 / (nx * ny * nz))
return
end
c-----------------------------------------------------------------------
c
c Computes negative of 3d residuals ...
c
c-----------------------------------------------------------------------
subroutine cmres3(mres,u,rhs,nx,ny,nz,h)
implicit none
integer nx, ny, nz
real*8 mres(nx,ny,nz), u(nx,ny,nz),
* rhs(nx,ny,nz)
real*8 h
integer i, j, k
real*8 hm2
hm2 = 1.0d0 / (h * h)
do k = 2 , nz - 1
do j = 2 , ny - 1
do i = 2 , nx - 1
mres(i,j,k) = rhs(i,j,k) -
* hm2 * (u(i+1,j,k) + u(i-1,j,k) +
* u(i,j+1,k) + u(i,j-1,k) +
* u(i,j,k+1) + u(i,j,k-1) -
* 6.0d0*u(i,j,k))
end do
end do
end do
call d3lbs(mres,0.0d0,nx,ny,nz)
return
end
c-----------------------------------------------------------------------
c
c Helper routine for argument checking ...
c
c This version doesn't demand that
c
c num = prime * 2^n
c
c-----------------------------------------------------------------------
integer function checkfacts(num,maxbase,base,npow2)
implicit none
integer ivprimedecomp
integer num, maxbase, base, npow2
integer maxnfact
parameter ( maxnfact = 100 )
integer vfact(maxnfact)
integer nfact, i
logical ltrace
parameter ( ltrace = .false. )
if( ivprimedecomp(num,vfact,nfact) .gt. 0 ) then
npow2 = 0
do i = 1 , nfact - 1
if( vfact(i) .ne. 2 ) goto 100
npow2 = npow2 + 1
end do
100 continue
base = 1
do i = npow2 + 1 , nfact
base = base * vfact(i)
end do
if( base .le. maxbase ) then
checkfacts = npow2
else
checkfacts = -1
end if
else
checkfacts = -1
end if
if( ltrace ) then
write(*,*) 'checkfacts: Number being factored ', num
call ivdump(vfact,nfact,'factors',6)
end if
return
end
c-----------------------------------------------------------------------
c Computes and dumps true residuals
c-----------------------------------------------------------------------
subroutine cres1_trace(label,unit,res,u,rhs,nx,h)
implicit none
character*(*) label
integer nx, unit
real*8 res(*), u(*), rhs(*)
real*8 h
integer shape(1)
real*8 vh(1)
shape(1) = nx
vh(1) = h
call cres(res,u,rhs,1,shape,vh)
write(unit,1000) label
1000 format('cres1_trace: ',a)
call dvdump(u,nx,'u',unit)
call dvdump(res,nx,'residual',unit)
return
end
c-----------------------------------------------------------------------
c Computes true (non-negated) residual.
c-----------------------------------------------------------------------
subroutine cres(res,u,rhs,rank,shape,vh)
implicit none
integer rank
real*8 res(*), u(*), rhs(*),
* vh(rank)
integer shape(rank)
call cmres(res,u,rhs,rank,shape,vh)
if( rank .eq. 1 ) then
call dvsm(res,-1.0d0,res,shape(1))
else if( rank .eq. 2 ) then
call dvsm(res,-1.0d0,res,shape(1)*shape(2))
else if( rank .eq. 3 ) then
call dvsm(res,-1.0d0,res,shape(1)*shape(2)*shape(3))
else
write(0,*) 'cmres: rank = ', rank,' not implemented'
end if
return
end
c-----------------------------------------------------------------------
c Front end to cmres1, cmres2, cmres3 ...
c-----------------------------------------------------------------------
subroutine cmres(mres,u,rhs,rank,shape,vh)
implicit none
integer rank
real*8 mres(*), u(*), rhs(*),
* vh(rank)
integer shape(rank)
if( rank .eq. 1 ) then
call cmres1(mres,u,rhs,shape(1), vh(1))
else if( rank .eq. 2 ) then
call cmres2(mres,u,rhs,shape(1),shape(2), vh(1))
else if( rank .eq. 3 ) then
call cmres3(mres,u,rhs,shape(1),shape(2),shape(3),vh(1))
else
write(0,*) 'cmres: rank = ', rank,' not implemented'
end if
return
end
c-----------------------------------------------------------------------
c Eventually should be incorporated into dmatlib.f
c
c Multiply matrix with scalar ...
c-----------------------------------------------------------------------
subroutine dmsm(a1,s1,a2,d1,d2)
implicit none
integer d1, d2
real*8 a1(d1,d2), a2(d1,d2)
real*8 s1
call dvsm(a1,s1,a2,d1*d2)
return
end