2
$\begingroup$

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:

  1. Solve independently on each $i$ processor: $A_i[\tilde{x}^i, v^i, w^i] = [d^i, \alpha e_0, \alpha e_{m-1}]$
  2. Send $\tilde{x}^i_0, v_0^i$ from the $i^{\text{th}}$ node to the $(i-1)^{\text{th}})$ node.
  3. 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} $
  4. Send $y_{2i+1}$ from the $i^{\text{th}}$ to the $(i+1)^{\text{th}}$ node.
  5. $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
$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.