lookup.F90

Go to the documentation of this file.

!! Copyright (C) 2009 X. Andrade
!!
!! 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 lookup_oct_m
  use debug_oct_m
  use global_oct_m
  use messages_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public ::                &
    lookup_t,              &
    lookup_init,           &
    lookup_end,            &
    lookup_get_list
 
  type lookup_t
    private
    integer            :: nobjs = 0
    integer            :: dim = 0
    real(real64), allocatable :: pos(:, :)
  end type lookup_t
 
contains
 
  subroutine lookup_init(this, dim, nobjs, pos)
    type(lookup_t), intent(inout) :: this
    integer,        intent(in)    :: dim
    integer,        intent(in)    :: nobjs
    real(real64),   intent(in)    :: pos(:, :)
 
    push_sub(lookup_init)
 
    this%nobjs = nobjs
    this%dim = dim
    safe_allocate(this%pos(1:this%dim, 1:this%nobjs))
 
    this%pos(1:this%dim, 1:this%nobjs) = pos(1:this%dim, 1:this%nobjs)
 
    pop_sub(lookup_init)
  end subroutine lookup_init
 
  ! -----------------------------------------
  subroutine lookup_end(this)
    type(lookup_t), intent(inout) :: this
 
    push_sub(lookup_end)
    safe_deallocate_a(this%pos)
 
    pop_sub(lookup_end)
  end subroutine lookup_end
 
  ! ------------------------------------------
  subroutine lookup_get_list(this, npoint, points, radius, nlist, list)
    type(lookup_t),                 intent(in)   :: this
    integer,                        intent(in)   :: npoint
    real(real64),                   intent(in)   :: points(:, :)
    real(real64),                   intent(in)   :: radius
    integer,                        intent(out)  :: nlist(:)
    integer, optional, allocatable, intent(out)  :: list(:, :)
 
    real(real64) :: r2
    integer :: ii, ipoint
 
    ! No PUSH SUB, called too often.
 
    if (present(list)) then
      safe_allocate(list(1:this%nobjs, 1:npoint))
    end if
 
    nlist(1:npoint) = 0
 
    do ii = 1, this%nobjs
      do ipoint = 1, npoint
        r2 = sum((this%pos(1:this%dim, ii) - points(ipoint, 1:this%dim))**2)
        if (r2 < radius**2) then
          nlist(ipoint) = nlist(ipoint) + 1
!This is a PGI pragma to force the optimization level of this file to -O0.
!-O2 or below is needed for 10.5. -O1 or below is needed for 10.8.
!The line after the pragma causes a segmentation fault otherwise.
!pgi$r opt=0
          if (present(list)) list(nlist(ipoint), ipoint) = ii
        end if
      end do
    end do
 
  end subroutine lookup_get_list
 
end module lookup_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: