hypre_obj_BuildMatrixIJ Subroutine

    impure subroutine hypre_obj_BuildMatrixIJ(this)
      !> Sets the coefficients of the matrix.
      implicit none
      class(hypre_obj), intent(inout) :: this                                  !! Hypre machinery
      ! Work variables
      integer  :: ierr
      integer  :: n,i,j,k
      real(wp) :: coef
      integer  :: Ng(3)
      integer  :: st

      associate (x=>this%block%x, y=>this%block%y, z=>this%block%z,   &
        xm=>this%block%xm, ym=>this%block%ym, zm=>this%block%zm,      &
        dxm=>this%block%dxm, dym=>this%block%dym, dzm=>this%block%dzm,&
        lo=>this%block%lo, hi =>this%block%hi, irow=>this%irow)

        ! Get number of grid points in each direction
        call this%parallel%Max(hi(1),Ng(1))
        call this%parallel%Max(Ng(1)-lo(1)+1,Ng(1))
        call this%parallel%Max(hi(2),Ng(2))
        call this%parallel%Max(Ng(2)-lo(2)+1,Ng(2))
        call this%parallel%Max(hi(3),Ng(3))
        call this%parallel%Max(Ng(3)-lo(3)+1,Ng(3))

        ! Setup the values of the matrix
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)

              ! Set column ID and operator value
              n = 0

              ! Set off-diagonal coefficients in the x3-dir
              if (Ng(3).gt.1) then
                do st=-this%stm,this%stm
                  ! Treat diagonals seperately
                  if (st.eq.0) cycle
                  if (irow%cell(i,j,k+st).ne.-1) then
                    n = n+1
                    coef = this%mat(i,j,k,3,st)
                    this%cols(n) = irow%cell(i,j,k+st); this%values(n) = coef*this%h2
                  end if
                end do
              end if
              ! Set off-diagonal coefficients in the x2-dir
              if (Ng(2).gt.1) then
                do st=-this%stm,this%stm
                  ! Treat diagonals seperately
                  if (st.eq.0) cycle
                  if (irow%cell(i,j+st,k).ne.-1) then
                    n = n+1
                    coef = this%mat(i,j,k,2,st)
                    this%cols(n) = irow%cell(i,j+st,k); this%values(n) = coef*this%h2
                  end if
                end do
              end if
              ! Set off-diagonal coefficients in the x1-dir
              if (Ng(1).gt.1) then
                do st=-this%stm,this%stm
                  ! Treat diagonals seperately
                  if (st.eq.0) cycle
                  if (irow%cell(i+st,j,k).ne.-1) then
                    n = n+1
                    coef = this%mat(i,j,k,1,st)
                    this%cols(n) = irow%cell(i+st,j,k); this%values(n) = coef*this%h2
                  end if
                end do
              end if

              ! - Diagonal
              n = n+1
              coef = 0.0_wp
              if (Ng(1).gt.1) coef = coef + this%mat(i,j,k,1,0)
              if (Ng(2).gt.1) coef = coef + this%mat(i,j,k,2,0)
              if (Ng(3).gt.1) coef = coef + this%mat(i,j,k,3,0)
              this%cols(n) = irow%cell(i,j,k)  ; this%values(n) = coef*this%h2

              ! number of columns (stencil size) and ID of row to set
              this%tmpi(1) = n
              this%tmpi(2) = irow%cell(i,j,k)

              ! Set values one row at a time
              call HYPRE_IJMatrixAddToValues(this%matrix,1,this%tmpi(1),this%tmpi(2),this%cols(1:n),this%values(1:n),ierr)
            end do
          end do
        end do
      end associate

      return
    end subroutine hypre_obj_BuildMatrixIJ