I would like to use the PDD algorithm developed by Sun to solve tridiagonal matrices in parallel for the following compact finite difference scheme:
$ \begin{align} \dfrac{1}{4}f^{'}_{i-1} + f^{'}_i + \dfrac{1}{4}f^{'}_{i+1} &= \dfrac{3}{4} \dfrac{-f_{i-1} + f_{i+1}}{dx} \; \; \; \; \;\;\;\;\;\;\;\;\;\; \text{(Interior Points)} \\ f^{'}_0 + 3f^{'}_1 &= \dfrac{-17 f_0 + 9 f_1 + 9 f_2 - f_3}{6 \Delta x} \; \; \; \text{(Boundary points)} \end{align} $
The PDD algorithm consists in splitting the diagonal matrix into a block diagonal, $\tilde{A}$, and offset matrix, $\Delta A$. The following example is for a 8 grid points domain divided into 2 sub-domains (processors):
$ \begin{align} A &= \begin{bmatrix} 1 & 3 & & & \\ \alpha & 1 & \alpha \\ &\alpha & 1 & \alpha \\ &&\alpha & 1 & \alpha \\ &&&\alpha & 1 & \alpha \\ &&&&\alpha & 1 & \alpha \\ &&&&& \alpha & 1 & \alpha \\ &&&&&& 3 & 1 \end{bmatrix} \\ &= \begin{bmatrix} 1 & 3 & & & \\ \alpha & 1 & \alpha \\ &\alpha & 1 & \alpha \\ &&\alpha & 1 & \\ &&& & 1 & \alpha \\ &&&&\alpha & 1 & \alpha \\ &&&&& \alpha & 1 & \alpha \\ &&&&&& 3 & 1 \end{bmatrix} + \begin{bmatrix} 0& & & & \\ & 0 & & \\ && 0 & & \\ &&& 0 & \alpha \\ &&& \alpha & 0 & \\ &&&&& 0 & \\ &&&&&& 0 & \\ &&&&&&& 0 \end{bmatrix} \\ \end{align} \\ \Delta A = [\alpha \, e_4, \alpha \, e_3][e^T_3, e^T_4] = VE^T $
Where $e_i$ represents a column vector with its ith elements being 1 and zero otherwise for $ 0 < i < n-1$. With this splitting, we have:
$ x^{-1} = A^{-1} d = (\tilde{A} + VE^T)^{-1} d = \tilde{A}^{-1} d - \tilde{A}^{-1}(1+E^T\tilde{A}^{-1}V)^{-1}E^T\tilde{A}^{-1}d $
This generally leads to have a pentadiagonal reduced matrix solved in parallel. However, the author shows that for a sufficiently large number of grid points per processors, the block-diagonal terms are below machine precision and the matrix reduces into 2x2 independent block problems. This is even more true for diagonally dominant matrices.
The algorithm consists of the following steps:
- Solve independently on each $i$ processor: $A_i[\tilde{x}^i, v^i, w^i] = [d^i, \alpha e_0, \alpha e_{m-1}]$
- Send $\tilde{x}^i_0, v_0^i$ from the $i^{\text{th}}$ node to the $(i-1)^{\text{th}})$ node.
- Solve: $ \begin{bmatrix} 1 & w_{m-1}^i \\ v_0^{i+1} & 1 \end{bmatrix} \begin{bmatrix} y_{2i} \\ y_{2i+1} \end{bmatrix} = \begin{bmatrix} \tilde{x}^i_{m-1} \\ \tilde{x}_0^{i+1}\end{bmatrix} $
- Send $y_{2i+1}$ from the $i^{\text{th}}$ to the $(i+1)^{\text{th}}$ node.
- $x^i = \tilde{x}^i - \begin{bmatrix}v^i, w^i\end{bmatrix} \begin{bmatrix} y_{2i-1} \\ y_{2i} \end{bmatrix}$
I wrote a MPI fortran code which follows the algorithm. I tested the code with explicit finite difference and also in serial to make sure the parallelism and compact scheme were implemented correctly. In parallel, I am getting large errors right at the nodes boundary. I have verified that my implementation follows the algorithm perfectly. Anyone know where the problem may lie?
Here is the standalone PDD algorithm for the compact scheme:
subroutine tridiag_parallel(nPoints, dx, phi, dphi, src, des)
integer, intent(in) :: nPoints, src, des
real(dp), intent(in) :: dx
real(dp), dimension(-buf:), intent(in) :: phi
real(dp), dimension( 0:), intent(out) :: dphi
integer :: ii, min_x, max_x
integer :: recv_rqst, send_rqst
real(dp) :: a_im, c_ip1m, x_recv, v_recv
real(dp) :: y_2im1, y_2i, y_2ip1
real(dp), dimension(0:nPoints-1) :: Ai_low, Ai_mid, Ai_upp, e_m
real(dp), dimension(0:nPoints-1) :: x_tilde, v, w, rhs, rhs_v, rhs_w
Ai_low = alpha; Ai_mid = 1._dp; Ai_upp = alpha
y_2i = 0._dp; y_2ip1 = 0._dp; y_2im1 = 0._dp; x_recv = 0._dp; v_recv = 0._dp
! Build RHS
min_x = 0; max_x = nPoints-1
if(src == mpi_proc_null) then
min_x = 1; Ai_upp(0) = 3._dp
rhs(0) = 1._dp/(6._dp*dx) * (-17._dp*phi(0) + 9._dp*(phi(1)+phi(2)) - phi(3))
end if
if(des == mpi_proc_null) then
max_x = nPoints-2; Ai_low(nPoints-1) = 3._dp
rhs(nPoints-1) = 1._dp/(6._dp*dx) * (phi(nPoints-4) - 9._dp*(phi(nPoints-3)+phi(nPoints-2)) &
+ 17._dp*phi(nPoints-1))
end if
do ii=min_x,max_x
rhs(ii) = 3._dp/(4._dp*dx) * (-phi(ii-1) + phi(ii+1))
end do
! Build a_im, c_im, e_m, rhs_v, rhs_w
a_im = alpha; c_ip1m = alpha; e_m = 0._dp
if(src == mpi_proc_null) a_im = 0._dp
if(des == mpi_proc_null) c_ip1m = 0._dp
rhs_v = 0._dp; rhs_v(0) = a_im
rhs_w = 0._dp; rhs_w(nPoints-1) = c_ip1m
! Solve for x_tilde, v and w
call tridiag(nPoints, Ai_low, Ai_mid, Ai_upp, rhs, x_tilde)
call tridiag(nPoints, Ai_low, Ai_mid, Ai_upp, rhs_v, v)
call tridiag(nPoints, Ai_low, Ai_mid, Ai_upp, rhs_w, w)
! First communication: send tilde{x_0}, v_0 to previous processor
call mpi_isend(x_tilde(0), 1, mpi_dp, src, 200, mpi_comm_world, send_rqst, ierror)
call mpi_irecv(x_recv, 1, mpi_dp, des, 200, mpi_comm_world, recv_rqst, ierror)
call mpi_wait(send_rqst, stat, ierror); call mpi_wait(recv_rqst, stat, ierror)
call mpi_isend(v(0), 1, mpi_dp, src, 300, mpi_comm_world, send_rqst, ierror)
call mpi_irecv(v_recv, 1, mpi_dp, des, 300, mpi_comm_world, recv_rqst, ierror)
call mpi_wait(send_rqst, stat, ierror); call mpi_wait(recv_rqst, stat, ierror)
! Solve for y_2i, y_2ip1
y_2i = (x_tilde(nPoints-1) - w(nPoints-1)*x_recv) / (1._dp - w(nPoints-1)*v_recv)
y_2ip1 = x_recv - v_recv * y_2i
! Second communication: send y_2i to next processor
call mpi_isend(y_2i , 1, mpi_dp, des, 400, mpi_comm_world, send_rqst, ierror)
call mpi_irecv(y_2im1, 1, mpi_dp, src, 400, mpi_comm_world, recv_rqst, ierror)
call mpi_wait(send_rqst, stat, ierror); call mpi_wait(recv_rqst, stat, ierror)
! Derivative
do ii=0,nPoints-1
dphi(ii) = x_tilde(ii) - (v(ii)*y_2im1 + w(ii)*y_2i)
end do
end subroutine tridiag_parallel
Below is the actual program (not sure if necessary but here it is anyway):
program main
use, intrinsic :: iso_fortran_env, dp => real64
use mpi, mpi_dp => mpi_double_precision, stat => mpi_status_ignore
implicit none
! Grid
integer :: nx, grid_nx
integer, parameter :: buf = 2
real(dp) :: dx, proc_xmin
real(dp), parameter :: grid_xmin = -1._dp, grid_xmax = 1._dp
real(dp), allocatable, dimension(:) :: x, x_grid
! MPI
integer :: nProcs, myrank, ierror, Xsrc, Xdes
logical :: master_proc = .false.
! Compact
real(dp), parameter :: alpha = 1._dp/4._dp
! Subscripts => s: serial, p: parallel, e: exact
real(dp), allocatable, dimension(:) :: u_s, u_p, dudx_s, dudx_p, dudx_e
integer :: ii, nn, send_rqst, recv_rqst
! Init mpi
call mpi_init(ierror)
call mpi_comm_size(mpi_comm_world, nProcs, ierror)
call mpi_comm_rank(mpi_comm_world, myrank, ierror)
if(myrank == 0) master_proc = .true.
call init_grid
! Serial
if(master_proc) then
Xsrc = mpi_proc_null; Xdes = mpi_proc_null
call setup_case(grid_nx, x_grid, u_s, dudx_s)
call tridiag_parallel(grid_nx, dx, u_s, dudx_s, Xsrc, Xdes)
!call FD2_x(grid_nx, dx, u_s, dudx_s, Xsrc, Xdes)
end if
! Parallel
call init_blocks
call setup_case(nx, x, u_p, dudx_p)
call tridiag_parallel(nx, dx, u_p, dudx_p, Xsrc, Xdes)
!call FD2_x(nx, dx, u_p, dudx_p, Xsrc, Xdes)
call compute_exact
! Compare
call compare(dudx_s, dudx_p)
call compute_errors(dudx_e, dudx_p)
! Finalize
call mpi_finalize(ierror)
contains
subroutine init_grid
! Read in grid size
if(master_proc) then
print'(A)', "Enter the total number of grid points"
read(*,*) grid_nx
do nn=1,nProcs-1
call mpi_isend(grid_nx, 1, mpi_int, nn, 100, mpi_comm_world, send_rqst, ierror)
call mpi_wait(send_rqst, stat, ierror)
end do
else
call mpi_irecv(grid_nx, 1, mpi_int, 0, 100, mpi_comm_world, recv_rqst, ierror)
call mpi_wait(recv_rqst, stat, ierror)
end if
if(mod(grid_nx, nProcs) /= 0) then
print'(A)', "The number of points is not divisible by the number of processors."
print'(A)', "Stopping the simulation"
call mpi_abort(mpi_comm_world, 10, ierror)
end if
dx = (grid_xmax - grid_xmin) / real(grid_nx - 1, dp)
nx = grid_nx / nProcs
proc_xmin = grid_xmin + mod(myrank, nProcs) * dx * nx
allocate(x(0:nx-1), x_grid(0:grid_nx-1))
x = [ (proc_xmin + ii*dx, ii=0,nx-1) ]
x_grid = [ (grid_xmin + ii*dx, ii=0,grid_nx-1) ]
end subroutine init_grid
subroutine init_blocks
Xsrc = myrank-1; Xdes = myrank+1
if(mod(myrank, nProcs) == 0) Xsrc = mpi_proc_null
if(mod(myrank+1, nProcs) == 0) Xdes = mpi_proc_null
end subroutine
subroutine apply_bc(u)
real(dp), dimension(-2:), intent(inout) :: u
call mpi_sendrecv(u(0), 2, mpi_dp, Xsrc, 100, &
u(nx), 2, mpi_dp, Xdes, 100, mpi_comm_world, stat, ierror)
call mpi_sendrecv(u(nx-2), 2, mpi_dp, Xdes, 200, &
u(-2), 2, mpi_dp, Xsrc, 200, mpi_comm_world, stat, ierror)
end subroutine apply_bc
subroutine setup_case(nx, x, u, dudx)
integer, intent(in) :: nx
real(dp), dimension(0:), intent(in) :: x
real(dp), allocatable, dimension(:), intent(out) :: u, dudx
allocate(u(-2:nx+1), dudx(0:nx-1)); u = 0._dp
u(0:nx-1) = [(exp(-10._dp*x(ii)**2), ii=0,nx-1)]
call apply_bc(u)
end subroutine setup_case
subroutine compute_exact
allocate(dudx_e(0:nx-1))
dudx_e = [(-20._dp*x(ii)*exp(-10._dp*x(ii)**2), ii=0,nx-1)]
end subroutine compute_exact
! 1st order derivative, 2nd order explicit central difference
subroutine FD2_x(nx, dx, phi, dphi, Xsrc, Xdes)
integer, intent(in) :: nx, Xsrc, Xdes
real(dp), intent(in) :: dx
real(dp), dimension(-2:), intent(in) :: phi
real(dp), dimension( 0:), intent(inout) :: dphi
integer :: nn, min_x, max_x
min_x = 0; max_x = nx-1
! BC
if(Xsrc == mpi_proc_null) then
min_x = 1
dphi(0) = (phi(1)-phi(0)) / dx
end if
if(Xdes == mpi_proc_null) then
max_x = nx-2
dphi(nx-1) = (-phi(nx-1)+phi(nx-2)) / dx
end if
! Interior points
do nn=min_x,max_x
dphi(nn) = (phi(nn+1)-phi(nn-1)) / (2._dp*dx)
end do
end subroutine FD2_x
subroutine tridiag(n, a, b, c, r, u)
integer, intent(in) :: n
real(dp), dimension(0:), intent(in) :: a, b, c, r
real(dp), dimension(0:), intent(out) :: u
real(dp), dimension(0:n-1) :: gam
integer :: j
real(dp) :: bet
bet = b(0)
u(0) = r(0)/bet
do j = 1,n-1
gam(j) = c(j-1) / bet
bet = b(j) - a(j) * gam(j)
u(j) = (r(j) - a(j) * u(j-1))/bet
end do
do j =n-2,0,-1
u(j) = u(j) - gam(j+1) * u(j+1)
end do
end subroutine tridiag
subroutine compare(serial, parallel)
real(dp), dimension(0:), intent(in) :: serial, parallel
integer :: proc_id
real(dp) :: Ltwo, Linf
real(dp), dimension(0:nx-1,0:nProcs-1) :: U_all
real(dp), dimension(0:grid_nx-1) :: grid_parallel
if(master_proc) then
U_all(:,0) = parallel(:)
do proc_id=1,nProcs-1
call mpi_recv(U_all(:,proc_id), nx, mpi_dp, proc_id, proc_id+100, mpi_comm_world, &
stat, ierror)
end do
else
call mpi_send(parallel(0:nx-1), nx, mpi_dp, 0, myrank+100, mpi_comm_world, ierror)
end if
call mpi_barrier(mpi_comm_world, ierror)
if(master_proc) then
grid_parallel = [(U_all(mod(ii,nx), int(ii/nx)), ii=0,grid_nx-1)]
Ltwo = sqrt(sum((serial-grid_parallel)**2)/real(grid_nx,dp))
Linf = maxval(abs(serial-grid_parallel))
print'(A)',"Serial vs Parallel:"
print'(2(A,ES15.6))', "Ltwo =", Ltwo, "; Linf =", Linf
open(unit=100,file='postProcessing/compare.dat')
do ii=0,grid_nx-1
write(100,*) x_grid(ii), serial(ii), grid_parallel(ii), abs(grid_parallel(ii)-serial(ii))
end do
close(100)
end if
end subroutine compare
subroutine compute_errors(exact, finite)
real(dp), dimension(0:), intent(in) :: exact, finite
real(dp) :: m_sum, m_max, Ltwo, Linf
m_sum = sum((exact-finite)**2)
call mpi_reduce(m_sum, Ltwo, 1, mpi_dp, mpi_sum, 0, mpi_comm_world, ierror)
Ltwo = sqrt(Ltwo/real(nx,dp))
m_max = maxval(abs(exact-finite))
call mpi_reduce(m_max, Linf, 1, mpi_dp, mpi_max, 0, mpi_comm_world, ierror)
if(master_proc) then
print'(A)', "Finite difference vs Exact:"
print'(2(A,ES15.6))', "Ltwo =", Ltwo, "; Linf =", Linf
end if
end subroutine compute_errors
end program
