! Full Multigrid V-cycle for a one dimensional pde, u''(x) = g(x),
! u(0) = 1, u(1) = e on [0,1].
program FMV
! u contains the approximate solution to the current problem on all grids.
! f contains the right hand side for the current problem on all grids.
! gridLevelInfo contains the starting index (into u and
! f) of the current grid as well as the dimension of the current
! grid.
! fullSize is the total static storage necessary to maintain
! results for all grids during execution.
! nSweepsBefore is the number of relaxation sweeps to perform
! before injection of the residual to the next coarser level.
! nSweepsAfter is the number of relaxation sweeps to perform
! after correction to the solution on the current grid.
! g is the right hand side of the de.
implicit none
integer, parameter :: K = 15 ! The number of grid levels
integer, parameter :: fullSize = 2**(K+1) + K - 2
real u(fullSize), f(fullSize)
integer gridLevelInfo(K,2)
integer nSweepsBefore, nSweepsAfter
parameter (nSweepsBefore = 1)
parameter (nSweepsAfter = 1)
integer i, n ! Index variables
integer index ! Current displacement into u and f
integer fineIndex, coarseIndex ! Displacements into u and f
integer nFine, nCoarse ! Size of fine and coarse grids
integer level ! Current grid level
real h ! Current step size
interface
function g( x )
real g
real, intent( in ) :: x
end function g
end interface
! Initialization
h = 0.5 ! Begin on coarsest grid
index = 1
! Set up right hand side and displacements into the solution
! and right hand side arrays, u and f.
do i = 1, K
N = 2**i + 1 ! Size of current grid
gridLevelInfo(i,1) = index
gridLevelInfo(i,2) = N
! Initialize right hand side for all grids for the FMV
! cycle. Note that this initialization is used to solve
! the problem Au = f on coarser and coarser grids. It is
! not for residual correction.
call setF( f( index:index+N-1 ), h, g ) ! Note f is a vector in R**N
index = index + N
h = h/2
enddo
! Set up the boundary conditions for the coarsest grid.
! These are copied by interpolate() into finer grids.
index = gridLevelInfo(1,1)
N = gridLevelInfo(1,2)
call setCoarseBC( u( index:index+N-1 ) )
! Begin FMV-cycle
! Relax on coarsest grid.
! Note that this code is duplicated within the general loop. Here,
! it sets up the coarsest grid for the following loop.
! Notice that relax with one unknown is a direct solve on the
! coarsest grid.
call relax( u( index:index+N-1 ), f( index:index+N-1 ), 1 )
do level = 2, K
! Retrieve intial guess from previous coarser grid using
! interpolation. Note that this information has been
! calculated either by the direct solve above or the
! ascent phase of the V-cycle in the general case.
fineIndex = gridLevelInfo(level,1)
nFine = gridLevelInfo(level,2)
coarseIndex = gridLevelInfo(level-1,1)
nCoarse = gridLevelInfo(level-1,2)
call interpolate( u( fineIndex:fineIndex+nFine-1 ), &
u( coarseIndex:coarseIndex+nCoarse-1 ) )
! Begin the descent phase of a single V-cycle.
! Descend to the coarsest grid for the current problem.
do i = level, 2, -1
fineIndex = gridLevelInfo(i,1)
coarseIndex = gridLevelInfo(i-1,1)
nFine = gridLevelInfo(i,2) ! Dimension of fine grid
nCoarse = gridLevelInfo(i-1,2) ! Dimension of coarse grid
! Relax on current grid using the initial guess just
! interpolated from the previous coarser grid.
call relax( u( fineIndex:fineIndex+nFine-1 ), &
f( fineIndex:fineIndex+nFine-1 ), &
nSweepsBefore)
! Inject residual of current equation into next coarser
! level. Note that this overwrites the previous right
! hand side. After the call to the subroutine
! injectResidual(), the right hand side contains
! the residual, r. The solution to Ae = r, on this coarser
! grid yields the correction for the current problem at
! level i. The correction is added in the ascent phase of
! the V-cycle.
call injectResidual( u( fineIndex:fineIndex+nFine-1 ), &
f( fineIndex:fineIndex+nFine-1 ), &
f( coarseIndex:coarseIndex+nCoarse-1 ) )
! Force a zero initial guess for the residual equation,
! Ae = r for coarse level.
call fillZeros( u( coarseIndex:coarseIndex+nCoarse-1 ) )
enddo ! End downward phase of V-cycle
! Now at coarsest grid, N = 3.
! Solve the current problem exactly.
! Residual equations are solved yielding corrections to the
! current equation on the next finer grid.
! Notice that relax with one unknown is a direct solve on the
! coarsest grid.
index = gridLevelInfo(1,1)
N = gridLevelInfo(1,2)
call relax( u( index:index+N-1 ), f( index:index+N-1 ), 1 )
! Begin ascent phase of V-cycle.
! Ascend to the finest grid (level) for the current
! problem.
do i = 2, level
coarseIndex = gridLevelInfo(i-1,1)
fineIndex = gridLevelInfo(i,1)
nCoarse = gridLevelInfo(i-1,2)
nFine = gridLevelInfo(i,2)
! Get correction for current problem from previous
! level using interpolation.
! Add to the current solution of the current problem.
call correct( u( fineIndex:fineIndex+nFine-1 ), &
u( coarseIndex:coarseIndex+nCoarse-1 ) )
! Relax on the fine grid problem after the correction.
call relax( u( fineIndex:fineIndex+nFine-1 ), &
f( fineIndex:fineIndex+nFine-1 ), &
nSweepsAfter)
enddo ! End upward phase of V-cycle
enddo ! End nested iteration
! Out a summary of results.
index = gridLevelInfo( K,1 )
N = gridLevelInfo( K,2 )
call writeResults( u( index:index+N-1 ), K )
contains
! Set up the boundary conditions for the coarsest grid.
! These are copied by interpolate() into finer grids.
subroutine setCoarseBC( u )
implicit none
real u(0:)
integer N
N = size( u )
u( 0 ) = 1.0
u( N-1 ) = exp( 1.0 )
end subroutine setCoarseBC
! Set the right hand size, f, using the step size h and the right
! hand side of the pde, g(x).
subroutine setF( f, h, g )
implicit none
real f(0:)
integer N
real h
interface
function g( x )
real g
real, intent( in ) :: x
end function g
end interface
integer i
real hSquared
hSquared = h * h
N = size( f )
do i = 1, N-2
f(i) = hSquared * exp(i * h) ! Note, x = i * h
end do
end subroutine setF
! Set the solution vector, u, to zero so that the residual equation
! is solved with initial guess 0.
subroutine fillZeros( u )
implicit none
real u(0:)
u = 0.0
end subroutine fillZeros
! Relaxation method.
! The relaxation method currently used is Jacobi.
! Note that u(0) and u(N) are both 0.0.
subroutine relax( u, f, nSweeps )
implicit none
integer nSweeps
real u(0:), f(0:)
integer i, j
integer N
N = size( u )
do j = 1, nSweeps
do i = 1, N-2
u(i) = 0.5 * ( u(i-1) + u(i+1) - f(i) )
end do
end do
end subroutine relax
subroutine interpolate( uFine, uCoarse )
implicit none
! This subroutine interpolates uCoarse to uFine.
! At even-numbered fine grid points, the values are transferred
! directly from the coarse grid.
! At odd-numbered fine grid points, value is the average of the
! two adjacent coarse grid points.
! uFine(0) and uFine(nFine-1) are boundary points.
! Note that 0 and nFine-1 are even numbered fine grid points.
real uFine(0:), uCoarse(0:)
integer i
integer nFine, nCoarse
nFine = size( uFine ); nCoarse = size( uCoarse )
do i = 0, nCoarse-1
uFine(2*i) = uCoarse(i)
end do
do i = 1, nFine-2, 2
uFine(i) = 0.5 * ( uFine(i-1) + uFine(i+1) )
end do
end subroutine interpolate
subroutine correct( uFine, uCoarse )
implicit none
! Correct the solution uFine by adding the interpolated correction
! given by uCoarse.
! Note that correct() uses the same interpolation algorithm as
! interpolate. Code was duplicated to avoid extra storage and the
! overhead of another subroutine call.
! Note that this subroutine and interpolate() could be combined if
! uFine were initialized to zero before the call to interpolate
! in the descent phase of each V-cycle.
real uFine(0:), uCoarse(0:)
integer i
integer nFine, nCoarse
nFine = size( uFine ); nCoarse = size( uCoarse )
do i = 0, nCoarse-2
uFine(2*i) = uFine(2*i) + uCoarse(i)
uFine(2*i+1) = uFine(2*i+1) + 0.5 * &
( uCoarse(i) + uCoarse(i+1) )
end do
uFine(2*(nCoarse-1)) = uFine(2*(nCoarse-1)) + uCoarse(nCoarse-1)
end subroutine correct
subroutine injectResidual( uFine, fFine, fCoarse )
implicit none
! Calculate the residual on the fine grid and inject it down to the
! coarse grid using a full weighting restriction operator.
real uFine(0:), fFine(0:)
real fCoarse(0:)
integer i ! Index variable
integer nFine, nCoarse
nFine = size( uFine ); nCoarse = size( fCoarse )
! Calculate the (2*i+1)st, (2*i+2)nd, and (2*i+3)rd components of
! the residual on the fine grid.
! Inject this component of the residual into the right hand side
! for the residual equation on the coarse grid.
do i = 1, nCoarse-2
fCoarse(i) = 0.25 * ( fFine(2*i-1) - &
( uFine(2*i-2) - 2 * uFine(2*i-1) + uFine(2*i) ) + 2 * &
( fFine(2*i) - ( uFine(2*i-1) - 2 * uFine(2*i) + &
uFine(2*i+1) ) ) + fFine(2*i+1) - &
( uFine(2*i) - 2 * uFine(2*i+1) + uFine(2*i+2) ) )
end do
end subroutine injectResidual
subroutine writeResults( u, K )
implicit none
real u(0:)
integer K
integer N
N = size( u )
print *, '-----------------------------------------'
print *, 'Summary of Results'
print *, '------------------'
print *, 'K =', K
print *, 'N =', N
print *, 'h =', 1.0/(N-1)
print *, 'u(0) =', u(0)
print *, 'u(0.5) =', u( (N-1)/2 )
print *, 'u(1) =', u(N-1)
print *, '-----------------------------------------'
end subroutine writeResults
end ! End program FMV
! g(x) is the right hand side of the pde being solved where x is
! the independent variable.
function g(x)
real g
real, intent( in ) :: x
g = exp(x)
end function g