nl_operator.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module nl_operator_oct_m
  use accel_oct_m
  use batch_oct_m
  use boundaries_oct_m
  use debug_oct_m
  use global_oct_m
  use index_oct_m
  use iso_c_binding
  use math_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use operate_f_oct_m
  use par_vec_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use sort_oct_m
  use stencil_oct_m
  use types_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::                     &
    nl_operator_t,              &
    nl_operator_index_t,        &
    nl_operator_global_init,    &
    nl_operator_global_end,     &
    nl_operator_init,           &
    nl_operator_build,          &
    dnl_operator_operate,       &
    znl_operator_operate,       &
    dnl_operator_operate_batch, &
    znl_operator_operate_batch, &
    dnl_operator_operate_diag,  &
    znl_operator_operate_diag,  &
    nl_operator_end,            &
    nl_operator_get_index,      &
    nl_operator_output_weights, &
    nl_operator_np_zero_bc,     &
    nl_operator_allocate_gpu_buffers, &
    nl_operator_update_gpu_buffers, &
    nl_operator_remove_zero_weight_points, &
    nl_operator_build_symmetric_weights
 
  type nl_operator_index_t
    private
    integer              :: nri = 0
    integer, allocatable :: imin(:)
    integer, allocatable :: imax(:)
    integer, allocatable :: ri(:, :)
    integer, allocatable :: ri_pos(:,:,:)
    integer, allocatable :: ri_neg(:,:,:)
  end type nl_operator_index_t
 
  integer, public, parameter ::  &
    OP_GENERAL       = 1,   &
    op_symmetric     = 2,   &
    op_antisymmetric = 3
 
  type nl_operator_t
    private
    type(stencil_t),   public :: stencil
    type(mesh_t), pointer     :: mesh => null()  
    integer,      allocatable :: nn(:)
    integer,           public :: np = 0
    ! When running in parallel mode, the next three arrays are unique on each node.
    real(real64),   allocatable, public :: w(:,:)
 
    logical,          public :: const_w = .true.   
 
    type(accel_mem_t), public :: buff_weights
    type(accel_mem_t), public :: buff_half_weights
 
    integer                   :: symmetry = op_general 
 
    character(len=40) :: label
 
    integer, public :: nri = 0
    integer, allocatable, public :: ri(:,:)
    integer, allocatable, public :: rimap(:)
    integer, allocatable, public :: rimap_inv(:)
 
    integer :: npairs = 0
    real(real64),allocatable :: wpair(:)
    real(real64)             :: wcenter
    integer, allocatable :: ri_pos(:,:,:)
    integer, allocatable :: ri_neg(:,:,:)
    integer :: max_allocated_ri_pair = 0
 
    integer                   :: ninner = 0
    integer                   :: nouter = 0
 
    type(nl_operator_index_t) :: inner
    type(nl_operator_index_t) :: outer
 
    type(accel_kernel_t) :: kernel
    type(accel_mem_t) :: buff_imin
    type(accel_mem_t) :: buff_imax
    type(accel_mem_t) :: buff_ri
    type(accel_mem_t) :: buff_map
    type(accel_mem_t) :: buff_all
    type(accel_mem_t) :: buff_inner
    type(accel_mem_t) :: buff_outer
    type(accel_mem_t) :: buff_stencil
    type(accel_mem_t) :: buff_ip_to_xyz
    type(accel_mem_t) :: buff_xyz_to_ip
 
    ! For multigrid solvers
    type(nl_operator_t), public, pointer :: coarser => null()
 
  end type nl_operator_t
 
  integer, parameter :: &
    OP_FORTRAN = 0,  &
    op_vec     = 1,  &
    op_min     = op_fortran, &
    op_max     = op_vec
 
  integer, parameter ::     &
    OP_INVMAP    = 1,       &
    op_map       = 2,       &
    op_nomap     = 3
 
  integer, public, parameter :: OP_ALL = 3, op_inner = 1, op_outer = 2
 
  interface
    integer function op_is_available(opid, type)
      implicit none
      integer, intent(in) :: opid, type
    end function op_is_available
  end interface
 
  integer :: dfunction_global = -1
  integer :: zfunction_global = -1
  integer :: function_accel
  logical :: use_symmetries = .false.
 
contains
 
  ! ---------------------------------------------------------
  subroutine nl_operator_global_init(namespace)
    type(namespace_t),         intent(in)    :: namespace
 
    integer :: default
 
    push_sub(nl_operator_global_init)
 
    !%Variable OperateDouble
    !%Type integer
    !%Section Execution::Optimization
    !%Default optimized
    !%Description
    !% This variable selects the subroutine used to apply non-local
    !% operators over the grid for real functions.
    !%Option fortran 0
    !% The standard Fortran function.
    !%Option optimized 1
    !% This version is optimized using vector primitives (if available).
    !%End
 
    !%Variable OperateComplex
    !%Type integer
    !%Section Execution::Optimization
    !%Default optimized
    !%Description
    !% This variable selects the subroutine used to apply non-local
    !% operators over the grid for complex functions.
    !%Option fortran 0
    !% The standard Fortran function.
    !%Option optimized 1
    !% This version is optimized using vector primitives (if available).
    !%End
 
    default = op_vec
 
    call parse_variable(namespace, 'OperateDouble', default, dfunction_global)
    if (.not. varinfo_valid_option('OperateDouble', dfunction_global)) call messages_input_error(namespace, 'OperateDouble')
 
    call parse_variable(namespace, 'OperateComplex', default, zfunction_global)
    if (.not. varinfo_valid_option('OperateComplex', zfunction_global)) call messages_input_error(namespace, 'OperateComplex')
 
    !%Variable OperateUseSymmetries
    !%Type logical
    !%Section Execution::Optimization
    !%Default no
    !%Description
    !% This variable selects if the operators are built using symmetries or not.
    !%End
    call parse_variable(namespace, 'OperateUseSymmetries', .false., use_symmetries)
    if (accel_is_enabled() .and. use_symmetries) then
      call messages_not_implemented("OperateUseSymmetries=yes with GPUs")
    end if
 
    if (accel_is_enabled()) then
 
      !%Variable OperateAccel
      !%Type integer
      !%Default map
      !%Section Execution::Optimization
      !%Description
      !% This variable selects the subroutine used to apply non-local
      !% operators over the grid when an accelerator device is used.
      !%Option invmap 1
      !% The standard implementation ported to GPUs.
      !%Option map 2
      !% A different version, more suitable for GPUs.
      !%End
      call parse_variable(namespace, 'OperateAccel',  op_map, function_accel)
 
      call messages_obsolete_variable(namespace, 'OperateOpenCL', 'OperateAccel')
 
    end if
 
    pop_sub(nl_operator_global_init)
  end subroutine nl_operator_global_init
 
  ! ---------------------------------------------------------
 
  subroutine nl_operator_global_end()
    push_sub(nl_operator_global_end)
 
    pop_sub(nl_operator_global_end)
  end subroutine nl_operator_global_end
 
  ! ---------------------------------------------------------
  subroutine nl_operator_init(op, label, symm)
    type(nl_operator_t), intent(inout) :: op
    character(len=*),    intent(in)    :: label
    integer, optional,   intent(in)    :: symm
 
    push_sub(nl_operator_init)
 
    op%label = label
    op%symmetry = op_general
    if (use_symmetries) then
      op%symmetry = optional_default(symm, op_general)
    end if
 
    pop_sub(nl_operator_init)
  end subroutine nl_operator_init
 
  ! ---------------------------------------------------------
  subroutine nl_operator_build(space, mesh, op, np, const_w, regenerate)
    class(space_t),       intent(in)    :: space
    type(mesh_t), target, intent(in)    :: mesh
    type(nl_operator_t),  intent(inout) :: op
    integer,              intent(in)    :: np
    logical, optional,    intent(in)    :: const_w
    logical, optional,    intent(in)    :: regenerate
 
    integer :: ii, jj, p1(space%dim), time, current
    integer, allocatable :: st1(:), st2(:), st1r(:)
    integer :: ir, maxp, iinner, iouter
    logical :: change, force_change
    character(len=200) :: flags
    integer, allocatable :: inner_points(:), outer_points(:), all_points(:)
 
    push_sub(nl_operator_build)
 
    if (mesh%parallel_in_domains .and. .not. const_w) then
      call messages_experimental('Domain parallelization with curvilinear coordinates')
    end if
 
    assert(np > 0)
 
    ! store values in structure
    op%np       = np
    op%mesh     => mesh
    if (.not. optional_default(regenerate, .false.)) then
      op%const_w  = optional_default(const_w, .false.)
    end if
 
    ! allocate weights op%w
    if (op%const_w) then
      safe_allocate(op%w(1:op%stencil%size, 1))
      message(1) = 'Debug: nl_operator_build: working with constant weights.'
      call messages_info(1, debug_only=.true.)
    else
      safe_allocate(op%w(1:op%stencil%size, 1:op%np))
      message(1) = 'Debug: nl_operator_build: working with non-constant weights.'
      call messages_info(1, debug_only=.true.)
    end if
 
    ! set initially to zero
    op%w = m_zero
 
    ! Build lookup table
    safe_allocate(st1(1:op%stencil%size))
    safe_allocate(st1r(1:op%stencil%size))
    safe_allocate(st2(1:op%stencil%size))
 
    op%nri = 0
    do time = 1, 2
      st2 = 0
      do ii = 1, np
        p1 = 0
        call mesh_local_index_to_coords(mesh, ii, p1)
 
        do jj = 1, op%stencil%size
          ! Get local index of p1 plus current stencil point.
          st1(jj) = mesh_local_index_from_coords(mesh, p1 + op%stencil%points(:, jj))
 
          assert(st1(jj) > 0)
        end do
 
        st1(1:op%stencil%size) = st1(1:op%stencil%size) - ii
 
        change = any(st1 /= st2)
 
        !the next is to detect when we move from a point that does not
        !have boundary points as neighbours to one that has
        force_change = any(st1 + ii > mesh%np) .and. all(st2 + ii - 1 <= mesh%np)
 
        ! if the stencil changes
        if (change .or. force_change) then
          !store it
          st2(:) = st1(:)
 
          !first time, just count
          if (time == 1) op%nri = op%nri + 1
 
          !second time, store
          if (time == 2) then
            current = current + 1
            op%ri(1:op%stencil%size, current) = st1(1:op%stencil%size)
          end if
        end if
 
        if (time == 2) op%rimap(ii) = current
 
      end do
 
      !after counting, allocate
      if (time == 1) then
        safe_deallocate_a(op%ri)
        safe_deallocate_a(op%rimap)
        safe_deallocate_a(op%rimap_inv)
 
        safe_allocate(op%ri(1:op%stencil%size, 1:op%nri))
        safe_allocate(op%rimap(1:op%np))
        safe_allocate(op%rimap_inv(1:op%nri + 1))
        op%ri        = 0
        op%rimap     = 0
        op%rimap_inv = 0
        current      = 0
 
        ! the sizes
        if (mesh%use_curvilinear) then
          safe_allocate(op%nn(1:op%nri))
          ! for the moment all the sizes are the same
          op%nn = op%stencil%size
        end if
      end if
 
    end do
 
    !the inverse mapping
    op%rimap_inv(1) = 0
    do jj = 1, op%np
      op%rimap_inv(op%rimap(jj) + 1) = jj
    end do
    op%rimap_inv(op%nri + 1) = op%np
 
    safe_deallocate_a(st1)
    safe_deallocate_a(st1r)
    safe_deallocate_a(st2)
 
    if (op%mesh%parallel_in_domains) then
      !now build the arrays required to apply the nl_operator by parts
 
      !count points
      op%inner%nri = 0
      op%outer%nri = 0
      do ir = 1, op%nri
        maxp = op%rimap_inv(ir + 1) + maxval(op%ri(1:op%stencil%size, ir))
        if (maxp <= np) then
          !inner point
          op%inner%nri = op%inner%nri + 1
          assert(op%inner%nri <= op%nri)
        else
          !outer point
          op%outer%nri = op%outer%nri + 1
          assert(op%outer%nri <= op%nri)
        end if
      end do
 
      assert(op%inner%nri + op%outer%nri == op%nri)
 
      if (optional_default(regenerate, .false.)) then
        safe_deallocate_a(op%inner%imin)
        safe_deallocate_a(op%inner%imax)
        safe_deallocate_a(op%inner%ri)
        safe_deallocate_a(op%outer%imin)
        safe_deallocate_a(op%outer%imax)
        safe_deallocate_a(op%outer%ri)
      end if
      safe_allocate(op%inner%imin(1:op%inner%nri + 1))
      safe_allocate(op%inner%imax(1:op%inner%nri))
      safe_allocate(op%inner%ri(1:op%stencil%size, 1:op%inner%nri))
 
      safe_allocate(op%outer%imin(1:op%outer%nri + 1))
      safe_allocate(op%outer%imax(1:op%outer%nri))
      safe_allocate(op%outer%ri(1:op%stencil%size, 1:op%outer%nri))
 
      !now populate the arrays
      iinner = 0
      iouter = 0
      do ir = 1, op%nri
        maxp = op%rimap_inv(ir + 1) + maxval(op%ri(1:op%stencil%size, ir))
        if (maxp <= np) then
          !inner point
          iinner = iinner + 1
          op%inner%imin(iinner) = op%rimap_inv(ir)
          op%inner%imax(iinner) = op%rimap_inv(ir + 1)
          op%inner%ri(1:op%stencil%size, iinner) = op%ri(1:op%stencil%size, ir)
        else
          !outer point
          iouter = iouter + 1
          op%outer%imin(iouter) = op%rimap_inv(ir)
          op%outer%imax(iouter) = op%rimap_inv(ir + 1)
          op%outer%ri(1:op%stencil%size, iouter) = op%ri(1:op%stencil%size, ir)
        end if
      end do
 
      !verify that all points in the inner operator are actually inner
      do ir = 1, op%inner%nri
        do ii = op%inner%imin(ir) + 1, op%inner%imax(ir)
          assert(all(ii + op%inner%ri(1:op%stencil%size, ir) <= mesh%np))
        end do
      end do
 
    end if
 
    if (accel_is_enabled() .and. op%const_w) then
 
      write(flags, '(i5)') op%stencil%size
      flags='-DNDIM=3 -DSTENCIL_SIZE='//trim(adjustl(flags))
 
      if (op%mesh%parallel_in_domains) flags = '-DINDIRECT '//trim(flags)
 
      select case (function_accel)
      case (op_invmap)
        call accel_kernel_build(op%kernel, 'operate.cu', 'operate', flags)
      case (op_map)
        call accel_kernel_build(op%kernel, 'operate.cu', 'operate_map', flags)
      end select
 
      ! conversion to i8 needed to avoid integer overflow
      call accel_create_buffer(op%buff_ri, accel_mem_read_only, type_integer, int(op%nri, int64)*op%stencil%size)
      call accel_write_buffer(op%buff_ri, op%stencil%size, op%nri, op%ri)
 
      select case (function_accel)
      case (op_invmap)
        call accel_create_buffer(op%buff_imin, accel_mem_read_only, type_integer, op%nri)
        call accel_write_buffer(op%buff_imin, op%nri, op%rimap_inv(1:))
        call accel_create_buffer(op%buff_imax, accel_mem_read_only, type_integer, op%nri)
        call accel_write_buffer(op%buff_imax, op%nri, op%rimap_inv(2:))
 
      case (op_map)
 
        call accel_create_buffer(op%buff_map, accel_mem_read_only, type_integer, pad(op%mesh%np, accel_max_block_size()))
        call accel_write_buffer(op%buff_map, op%mesh%np, (op%rimap - 1)*op%stencil%size)
 
        if (op%mesh%parallel_in_domains) then
 
          safe_allocate(inner_points(1:op%mesh%np))
          safe_allocate(outer_points(1:op%mesh%np))
          safe_allocate(all_points(1:op%mesh%np))
 
          op%ninner = 0
          op%nouter = 0
 
          do ii = 1, op%mesh%np
            all_points(ii) = ii - 1
            maxp = ii + maxval(op%ri(1:op%stencil%size, op%rimap(ii)))
            if (maxp <= op%mesh%np) then
              op%ninner = op%ninner + 1
              inner_points(op%ninner) = ii - 1
            else
              op%nouter = op%nouter + 1
              outer_points(op%nouter) = ii - 1
            end if
          end do
 
          call accel_create_buffer(op%buff_all, accel_mem_read_only, type_integer, pad(op%mesh%np, accel_max_block_size()))
          call accel_write_buffer(op%buff_all, op%mesh%np, all_points)
 
          call accel_create_buffer(op%buff_inner, accel_mem_read_only, type_integer, pad(op%ninner, accel_max_block_size()))
          call accel_write_buffer(op%buff_inner, op%ninner, inner_points)
 
          call accel_create_buffer(op%buff_outer, accel_mem_read_only, type_integer, pad(op%nouter, accel_max_block_size()))
          call accel_write_buffer(op%buff_outer, op%nouter, outer_points)
 
          safe_deallocate_a(inner_points)
          safe_deallocate_a(outer_points)
          safe_deallocate_a(all_points)
 
        end if
      end select
    end if
 
    pop_sub(nl_operator_build)
 
  end subroutine nl_operator_build
 
  ! ---------------------------------------------------------
  subroutine nl_operator_output_weights(this)
    type(nl_operator_t), intent(inout)  :: this
 
    integer :: istencil, idir
 
    push_sub(nl_operator_output_weights)
 
    write(message(1), '(3a)') 'Debug info: Finite difference weights for ', trim(this%label), '.'
    write(message(2), '(a)')  '            Spacing:'
    do idir = 1, this%mesh%box%dim
      write(message(2), '(a,f16.8)') trim(message(2)), this%mesh%spacing(idir)
    end do
    call messages_info(2, debug_only=.true.)
 
    do istencil = 1, this%stencil%size
      select case(this%mesh%box%dim)
      case(1)
        write(message(1), '(a,i3,1i4,f25.10)') '      ', istencil, this%stencil%points(1:1, istencil), this%w(istencil, 1)
      case(2)
        write(message(1), '(a,i3,2i4,f25.10)') '      ', istencil, this%stencil%points(1:2, istencil), this%w(istencil, 1)
      case(3)
        write(message(1), '(a,i3,3i4,f25.10)') '      ', istencil, this%stencil%points(1:3, istencil), this%w(istencil, 1)
      end select
      call messages_info(1, debug_only=.true.)
    end do
 
    pop_sub(nl_operator_output_weights)
 
  end subroutine nl_operator_output_weights
 
  ! ---------------------------------------------------------
  subroutine nl_operator_end(op)
    type(nl_operator_t), intent(inout) :: op
 
    push_sub(nl_operator_end)
 
    if (accel_is_enabled() .and. op%const_w) then
      call nl_operator_clear_gpu_buffers(op)
    end if
 
    safe_deallocate_a(op%inner%imin)
    safe_deallocate_a(op%inner%imax)
    safe_deallocate_a(op%inner%ri)
    safe_deallocate_a(op%outer%imin)
    safe_deallocate_a(op%outer%imax)
    safe_deallocate_a(op%outer%ri)
 
    safe_deallocate_a(op%w)
 
    safe_deallocate_a(op%ri)
    safe_deallocate_a(op%rimap)
    safe_deallocate_a(op%rimap_inv)
    safe_deallocate_a(op%nn)
 
    safe_deallocate_a(op%wpair)
    safe_deallocate_a(op%ri_pos)
    safe_deallocate_a(op%ri_neg)
 
    safe_deallocate_a(op%inner%ri_pos)
    safe_deallocate_a(op%inner%ri_neg)
    safe_deallocate_a(op%outer%ri_pos)
    safe_deallocate_a(op%outer%ri_neg)
    call stencil_end(op%stencil)
 
    pop_sub(nl_operator_end)
  end subroutine nl_operator_end
 
 
  subroutine nl_operator_clear_gpu_buffers(op)
    type(nl_operator_t), intent(inout) :: op
 
    push_sub(nl_operator_clear_gpu_buffers)
 
    call accel_free_buffer(op%buff_ri)
    select case (function_accel)
    case (op_invmap)
      call accel_free_buffer(op%buff_imin)
      call accel_free_buffer(op%buff_imax)
 
    case (op_map)
      call accel_free_buffer(op%buff_map)
      if (op%mesh%parallel_in_domains) then
        call accel_free_buffer(op%buff_all)
        call accel_free_buffer(op%buff_inner)
        call accel_free_buffer(op%buff_outer)
      end if
 
    case (op_nomap)
      call accel_free_buffer(op%buff_map)
      call accel_free_buffer(op%buff_stencil)
      call accel_free_buffer(op%buff_xyz_to_ip)
      call accel_free_buffer(op%buff_ip_to_xyz)
    end select
 
    call accel_free_buffer(op%buff_weights)
    call accel_free_buffer(op%buff_half_weights)
 
    pop_sub(nl_operator_clear_gpu_buffers)
  end subroutine nl_operator_clear_gpu_buffers
 
  ! ---------------------------------------------------------
  integer pure function nl_operator_get_index(op, is, ip) result(res)
    type(nl_operator_t), intent(in)   :: op
    integer,             intent(in)   :: is
    integer,             intent(in)   :: ip
 
    res = ip + op%ri(is, op%rimap(ip))
  end function nl_operator_get_index
 
  ! ---------------------------------------------------------
 
  subroutine nl_operator_allocate_gpu_buffers(op)
    type(nl_operator_t), intent(inout)   :: op
 
    push_sub(nl_operator_allocate_gpu_buffers)
 
    ! Update the GPU weights
    if (accel_is_enabled() .and. op%const_w) then
      call accel_create_buffer(op%buff_weights, accel_mem_read_only, type_float, op%stencil%size)
      call accel_create_buffer(op%buff_half_weights, accel_mem_read_only, type_float, op%stencil%size)
    end if
 
    pop_sub(nl_operator_allocate_gpu_buffers)
  end subroutine nl_operator_allocate_gpu_buffers
 
 
  ! ---------------------------------------------------------
 
  subroutine nl_operator_update_gpu_buffers(op)
    type(nl_operator_t), intent(inout)   :: op
 
    push_sub(nl_operator_update_gpu_buffers)
 
    ! Update the GPU weights
    if (accel_is_enabled() .and. op%const_w) then
      call accel_write_buffer(op%buff_weights, op%stencil%size, op%w(:, 1))
      call accel_write_buffer(op%buff_half_weights, op%stencil%size, -m_half*op%w(:, 1))
    end if
 
    pop_sub(nl_operator_update_gpu_buffers)
  end subroutine nl_operator_update_gpu_buffers
 
  ! ---------------------------------------------------------
 
  integer pure function nl_operator_np_zero_bc(op) result(np_bc)
    type(nl_operator_t), intent(in)   :: op
 
    integer :: jj, ii
 
    np_bc = 0
    do jj = 1, op%nri
      ii = op%rimap_inv(jj + 1) + maxval(op%ri(1:op%stencil%size, jj))
      np_bc = max(np_bc, ii)
    end do
 
  end function nl_operator_np_zero_bc
 
 
  ! ---------------------------------------------------------
  subroutine nl_operator_remove_zero_weight_points(op, space, mesh)
    type(nl_operator_t), intent(inout)   :: op
    type(space_t),       intent(in)      :: space
    class(mesh_t),       intent(in)      :: mesh
 
    integer :: ip, size
    real(real64), parameter :: tol = 1.0e-14_real64
    real(real64) :: max_weight, new_w(op%stencil%size)
    integer :: new_points(space%dim, op%stencil%size)
 
    if (.not. op%const_w) return
 
    push_sub(nl_operator_remove_zero_weight_points)
 
    max_weight = maxval(abs(op%w(:, 1)))
    size = 0
    do ip = 1, op%stencil%size
      if (abs(op%w(ip, 1)) > tol * max_weight) then
        size = size +1
        new_w(size) = op%w(ip, 1)
        new_points(:,size) = op%stencil%points(:, ip)
      end if
    end do
 
    ! We regenerate the stencil without the zero-weight points
    op%stencil%size = size
    safe_deallocate_a(op%stencil%points)
    safe_allocate(op%stencil%points(space%dim, op%stencil%size))
    op%stencil%points(:, :) = new_points(:, 1:size)
    safe_deallocate_a(op%w)
    call nl_operator_clear_gpu_buffers(op)
    call nl_operator_build(space, mesh, op, mesh%np, const_w=op%const_w, regenerate=.true.)
    op%w(1:size, 1) = new_w(1:size)
 
    !Update Stencil%center
    call stencil_init_center(op%stencil)
 
    pop_sub(nl_operator_remove_zero_weight_points)
  end subroutine nl_operator_remove_zero_weight_points
 
  subroutine group_by_pairs_sym(size, ldf, offsets, wre, ri, nri, npairs, wpair, pair_pos, pair_neg, wcenter)
    integer,                    intent(in)  :: size
    integer,                    intent(in)  :: ldf
    integer,                    intent(in)  :: offsets(:, :)
    real(real64),               intent(in)  :: wre(:)
    integer,                    intent(in)  :: ri(:, :)
    integer,                    intent(in)  :: nri
    integer,                    intent(out) :: npairs
    real(real64),             intent(inout) :: wpair(:)
    integer,                  intent(inout) :: pair_pos(:,:), pair_neg(:,:)
    real(real64),               intent(out) :: wcenter
 
    logical, allocatable :: used(:)
    integer :: i, j, ndim, s
    logical :: same
    integer, allocatable :: idx(:)
 
    real(real64), parameter :: tol = 1.0e-12_real64
 
    push_sub(group_by_pairs_sym)
 
    assert(mod(size,2) == 1)
 
    safe_allocate(used(1:size))
    used = .false.
    npairs = 0
 
    ndim = ubound(offsets, dim=1)
 
    safe_allocate(idx(1:size))
    call robust_sort_by_abs(wre, offsets, idx)
 
    do i = 1, size
      if (used(i)) cycle
 
      if (all(offsets(:, idx(i))==0)) then
        wcenter = wre(idx(i))
        used(i) = .true.
        cycle
      end if
 
      ! Try to find symmetric partner j
      do j = i+1, size
        if (used(j)) cycle
 
        ! Weight equality
        same = abs(wre(idx(i)) - wre(idx(j))) <= tol
        if (.not. same) cycle
 
        ! Offsets equal and opposite
        if (any(offsets(:,idx(j))+offsets(:, idx(i)) /= 0)) cycle
 
        npairs = npairs + 1
        do s = 1, nri
          pair_pos(npairs, s) = ri(idx(i), s) * 2**ldf
          pair_neg(npairs, s) = ri(idx(j), s) * 2**ldf
        end do
        wpair(npairs)  = m_half*(wre(idx(i)) + wre(idx(j)))
 
        used(i) = .true.
        used(j) = .true.
        exit
      end do
    end do
 
    assert(npairs == size/2)
 
    safe_deallocate_a(idx)
 
    pop_sub(group_by_pairs_sym)
  end subroutine group_by_pairs_sym
 
  subroutine group_by_pairs_antisym(size, ldf, offsets, wre, ri, nri, npairs, wpair, pair_pos, pair_neg)
    integer,                    intent(in)  :: size
    integer,                    intent(in)  :: ldf
    integer,                    intent(in)  :: offsets(:, :)
    real(real64),               intent(in)  :: wre(:)
    integer,                    intent(in)  :: ri(:, :)
    integer,                    intent(in)  :: nri
    integer,                   intent(out)  :: npairs
    real(real64),             intent(inout) :: wpair(:)
    integer,                  intent(inout) :: pair_pos(:,:), pair_neg(:,:)
 
    logical, allocatable :: used(:)
    integer :: i, j, ndim, s
    logical :: same
    integer, allocatable :: idx(:)
 
    real(real64), parameter :: tol = 1.0e-12_real64
 
    push_sub(group_by_pairs_antisym)
 
    assert(mod(size,2) == 0)
 
    safe_allocate(used(1:size))
    used = .false.
    npairs = 0
 
    ndim = ubound(offsets, dim=1)
 
    safe_allocate(idx(1:size))
    call robust_sort_by_abs(wre, offsets, idx)
 
    ! Max pairs = n/2
    do i = 1, size
      if (used(i)) cycle
 
      ! Try to find symmetric partner j
      do j = i+1, size
        if (used(j)) cycle
 
        ! Weight equality
        same = abs(wre(idx(i)) + wre(idx(j))) <= tol
        if (.not. same) cycle
 
        ! Offsets equal and opposite
        if (any(offsets(:,idx(j))+offsets(:, idx(i)) /= 0)) cycle
 
        npairs = npairs + 1
        do s = 1, nri
          pair_pos(npairs, s) = ri(idx(i), s) * 2**ldf
          pair_neg(npairs, s) = ri(idx(j), s) * 2**ldf
        end do
        wpair(npairs)  = m_half*(wre(idx(i)) - wre(idx(j)))
 
        used(i) = .true.
        used(j) = .true.
        exit
      end do
    end do
 
    assert(npairs == size/2)
 
    safe_deallocate_a(idx)
    pop_sub(group_by_pairs_antisym)
  end subroutine group_by_pairs_antisym
 
  subroutine nl_operator_build_symmetric_weights(op, max_size)
    type(nl_operator_t), intent(inout)   :: op
    integer, optional,   intent(in)      :: max_size
 
    integer :: ldf, start, end, ipair
 
    if (op%symmetry == op_general) return
 
    push_sub(nl_operator_build_symmetric_weights)
 
    assert(op%const_w)
    assert(conf%target_states_block_size > 0)
 
    if(present(max_size)) then
      start = op%max_allocated_ri_pair + 1
      end = max_size
      call reallocate_array(op%ri_pos, op%stencil%size/2, op%nri, op%max_allocated_ri_pair, end)
      call reallocate_array(op%ri_neg, op%stencil%size/2, op%nri, op%max_allocated_ri_pair, end)
      if (op%mesh%parallel_in_domains) then
        call reallocate_array(op%inner%ri_pos, op%stencil%size/2, op%inner%nri, op%max_allocated_ri_pair, end)
        call reallocate_array(op%inner%ri_neg, op%stencil%size/2, op%inner%nri, op%max_allocated_ri_pair, end)
        call reallocate_array(op%outer%ri_pos, op%stencil%size/2, op%outer%nri, op%max_allocated_ri_pair, end)
        call reallocate_array(op%outer%ri_neg, op%stencil%size/2, op%outer%nri, op%max_allocated_ri_pair, end)
      end if
    else
      ! Max pairs = n/2
      safe_allocate(op%wpair(1:op%stencil%size/2))
      start = 1
      end = log2(conf%target_states_block_size)+1
 
      safe_allocate(op%ri_pos(1:op%stencil%size/2, 1:op%nri, 1:end))
      safe_allocate(op%ri_neg(1:op%stencil%size/2, 1:op%nri, 1:end))
      if (op%mesh%parallel_in_domains) then
        safe_allocate(op%inner%ri_pos(1:op%stencil%size/2, 1:op%inner%nri, 1:end))
        safe_allocate(op%inner%ri_neg(1:op%stencil%size/2, 1:op%inner%nri, 1:end))
        safe_allocate(op%outer%ri_pos(1:op%stencil%size/2, 1:op%outer%nri, 1:end))
        safe_allocate(op%outer%ri_neg(1:op%stencil%size/2, 1:op%outer%nri, 1:end))
      end if
    end if
    op%max_allocated_ri_pair = end
 
    do ldf = start-1, end-1
      select case(op%symmetry)
      case(op_symmetric)
        call group_by_pairs_sym(op%stencil%size, ldf, op%stencil%points, op%w(:,1), op%ri, op%nri, &
          op%npairs, op%wpair, op%ri_pos(:,:,ldf+1), op%ri_neg(:,:,ldf+1), op%wcenter)
        if (op%mesh%parallel_in_domains) then
          call group_by_pairs_sym(op%stencil%size, ldf, op%stencil%points, op%w(:,1), op%inner%ri, op%inner%nri, &
            op%npairs, op%wpair, op%inner%ri_pos(:,:,ldf+1), op%inner%ri_neg(:,:,ldf+1), op%wcenter)
          call group_by_pairs_sym(op%stencil%size, ldf, op%stencil%points, op%w(:,1), op%outer%ri, op%outer%nri, &
            op%npairs, op%wpair, op%outer%ri_pos(:,:,ldf+1), op%outer%ri_neg(:,:,ldf+1), op%wcenter)
        end if
      case(op_antisymmetric)
        call group_by_pairs_antisym(op%stencil%size, ldf, op%stencil%points, op%w(:,1), op%ri, op%nri, &
          op%npairs, op%wpair, op%ri_pos(:,:,ldf+1), op%ri_neg(:,:,ldf+1))
        if (op%mesh%parallel_in_domains) then
          call group_by_pairs_antisym(op%stencil%size, ldf, op%stencil%points, op%w(:,1), op%inner%ri, op%inner%nri, &
            op%npairs, op%wpair, op%inner%ri_pos(:,:,ldf+1), op%inner%ri_neg(:,:,ldf+1))
          call group_by_pairs_antisym(op%stencil%size, ldf, op%stencil%points, op%w(:,1), op%outer%ri, op%outer%nri, &
            op%npairs, op%wpair, op%outer%ri_pos(:,:,ldf+1), op%outer%ri_neg(:,:,ldf+1))
        end if
      end select
    end do
 
    if (.not. present(max_size)) then
      write(message(1), '(3a)') 'Debug info: Sorted weights for ', trim(op%label), '.'
      call messages_info(1, debug_only=.true.)
 
      do ipair = 1, op%npairs
        write(message(1), '(a,i3,f25.10,2(1x,i4))') '      ', ipair, op%wpair(ipair), op%ri_pos(ipair,1,1), op%ri_neg(ipair,1,1)
        call messages_info(1, debug_only=.true.)
      end do
    end if
 
    pop_sub(nl_operator_build_symmetric_weights)
  end subroutine nl_operator_build_symmetric_weights
 
  subroutine reallocate_array(ri, stencil_size, nri, old_size, new_size)
    integer,      allocatable, intent(inout) :: ri(:,:,:)
    integer,                   intent(in)    :: stencil_size, nri
    integer,                   intent(in)    :: old_size, new_size
 
    integer, allocatable :: tmp(:,:,:)
 
    safe_allocate_source_a(tmp, ri)
    safe_deallocate_a(ri)
    safe_allocate(ri(1:stencil_size, 1:nri, 1:new_size))
    ri(:,:,1:old_size) = tmp
    safe_deallocate_a(tmp)
  end subroutine reallocate_array
 
#include "undef.F90"
#include "real.F90"
#include "nl_operator_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "nl_operator_inc.F90"
 
end module nl_operator_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: