Below is a complete nested iteration V cycle multigrid code.
It appears in S. F. McCormick, Multigrid Methods, Frontiers in Applied
Mathematics, vol. 3, SIAM Books, Philadelphia, 1987.
Below is a complete nested iteration V cycle multigrid code.
It appears in S. F. McCormick, Multigrid Methods, Frontiers in Applied
Mathematics, vol. 3, SIAM Books, Philadelphia, 1987.
program fmv
implicit real*8 (a-h,o-z)
dimension u(8000), f(8000), ks(10), ex(8000)
do 500 i = 1,8000
u(i) = 0.
f(i) = 0.
ex(i) = 0.
500 continue
h = .5
kst = 1
nlev = 6
n = 2 ** nlev
write (*,5) n
5 format (// ' One FMV cycle for n = ', i3 // )
do 10 k = 1,nlev
ks(k) = kst
n = 2 ** k + 1
call fset ( f(kst), ex(kst), n, h )
kst = kst + n*n
h = h * .5
10 continue
do 300 n = 1,nlev
if ( n.gt.1 ) then
do 100 k = n,2,-1
kst = ks(k)
kstc = ks(k-1)
nf = 2 ** k + 1
nc = 2 ** (k-1) + 1
call relax ( u(kst), f(kst), nf )
call relax ( u(kst), f(kst), nf )
call resinj ( u(kst), f(kst), f(kstc), nf, nc )
100 continue
endif
do 200 k = 1,n
kst = ks(k)
kstf = ks(k+1)
nc = 2 ** k + 1
nf = 2 ** (k+1) + 1
c call relax ( u(kst), f(kst), nc )
call relax ( u(kst), f(kst), nc )
if ( k.eq.n ) call err ( u(kst), ex(kst), f(kst), nc )
if ( k.lt.nlev ) call intadd ( u(kst), u(kstf), nc, nf )
200 continue
300 continue
end
subroutine err ( u, ex, f, n )
implicit real*8 (a-h,o-z)
dimension u(n,n), f(n,n), ex(n,n)
data once / 0 /
h = 1. / (n-1)
emax = 0.
esum = 0.
rmax = 0.
rsum = 0.
n1 = n - 1
do 10 i = 2,n1
do 10 j = 2,n1
er = abs( u(i,j) - ex(i,j) )
if ( er.gt.emax ) emax = er
r = abs ( f(i,j) + u(i+1,j) + u(i-1,j) + u(i,j+1) +
* u(i,j-1) - 4. * u(i,j) )
if ( r.gt.rmax ) rmax = r
esum = esum + er*er
rsum = rsum + r*r
10 continue
eh = (h * esum) ** .5
rh = (h * rsum) ** .5
if ( once.eq.0 ) write (*,19) 'n', 'rmax', 'rh', 'emax', 'eh'
once = 1
write (*,20) n, rmax, rh, emax, eh
19 format ( 1x, a3, 4a12 /)
20 format ( 1x, i3, 4e12.3 )
return
end
subroutine fset ( g, ex, n, h )
implicit real*8 (a-h,o-z)
dimension g(n,n), ex(n,n)
h2 = h * h
n1 = n - 1
do 10 i = 2,n1
x = h * (i - 1)
x2 = x * x
xx = x2 * (x2 - 1.)
x6 = 6.*x2 - 1.
do 10 j = 2,n1
y = h * (j - 1)
y2 = y * y
yy = y2 * (y2 - 1.)
ex(i,j) = xx * yy
g(i,j) = -2. * h2 * (x6 * yy + (6. * y2 - 1.) * xx)
10 continue
return
end
subroutine relax ( v, g, n )
implicit real*8 (a-h,o-z)
dimension v(n,n), g(n,n)
n1 = n - 1
irb = 1
do 10 krb = 1,2
do 20 i = 2,n1
irb = 1 - irb
do 20 j = irb+2,n1,2
v(i,j) = (g(i,j) + v(i+1,j) + v(i-1,j) + v(i,j+1) +
* v(i,j-1)) * .25
20 continue
10 continue
return
end
subroutine resinj ( v, f, g, n, nc )
implicit real*8 (a-h,o-z)
dimension f(n,n), v(n,n), g(nc,nc)
n1 = nc - 1
do 10 i = 2,n1
i2 = 2*i - 1
do 10 j = 2,n1
j2 = 2*j - 1
t = f(i2,j2) + v(i2+1,j2) + v(i2-1,j2) + v(i2,j2+1) +
* v(i2,j2-1)
g(i,j) = 2.0 * t - 8.0 * v(i2,j2)
10 continue
return
end
subroutine intadd ( v, w, n, nf )
implicit real*8 (a-h,o-z)
dimension v(n,n), w(nf,nf)
do 10 i = 2,n
i2 = 2*i - 1
do 10 j = 2,n
j2 = 2*j - 1
w(i2,j2-1) = w(i2,j2-1) + (v(i,j)+v(i,j1-1)) * .5
w(i2-1,j2) = w(i2-1,j2) + (v(i,j)+v(i-1,j1)) * .5
10 continue
do 20 i = 1,n
do 20 j = 1,n
v(i,j) = 0.
20 continue
return
end