!===========================================
! A subroutine for assembling Laplace matrix
! in the compressed row storage format
! (sequential version)
!===========================================
use header
implicit none
! Variables as per assignment
integer, intent(in) :: m
type(Matrix), intent(inout) :: Delta
! Loop counters: total, horizontal, vertical, nonzeros, neighbour
integer :: irow, i, j, inz, next
! Total problem size
integer :: n
inz = 1 ! index of current nonzero element
do irow=1,n
! Calculate Cartesian index splitting
j = (irow-1)/m + 1
i = irow - (j-1)*m
! Init current row position
Delta%ii(irow) = inz
! Diagonal element
Delta%aa(inz) = -4.0_rk * m**2 ! -4/h^2
Delta%jj(inz) = irow ! column index = row index on diagonal
! Done with diagonal element, shift the counter of elements
inz = inz + 1
! Off-diagonal elements, including periodically looped ones
! left
Delta%aa(inz) = m**2 ! 1/h^2
next = i-1
if (next<1) next = m ! loop over if we are on a boundary
Delta%jj(inz) = next + (j-1)*m ! column index back in global range
! shift the nnz counter
inz = inz + 1
! right
Delta%aa(inz) = m**2
next = i+1
if (next>m) next = 1
Delta%jj(inz) = next + (j-1)*m
inz = inz + 1
! bottom
Delta%aa(inz) = m**2
next = j-1
if (next<1) next = m
Delta%jj(inz) = i + (next-1)*m
inz = inz + 1
! top
Delta%aa(inz) = m**2
next = j+1
if (next>m) next = 1
Delta%jj(inz) = i + (next-1)*m
inz = inz + 1
end do