Hi Nvidia team,
this is another problem report from our test suite. A UDDT (particle_set_t) that has a UDDT (particle_t) as component gets a real component assigned but segfaults. Probably the construction of this data type has failed. The code is attached here:
main_ut.f90 (16.0 KB)
and below. The segfault (in the debugger):
Type "apropos word" to search for commands related to "word"...
Reading symbols from ./whizard_test...
(gdb) run
Starting program: /scratch/reuter/local/packages/whizard/debug/whizard-3.0.2/_build/tests/unit_tests/arsch/test_36/whizard_test
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
b
Program received signal SIGSEGV, Segmentation fault.
particles::particle_set_momentum () at main_ut.f90:315
315 prt%p2 = p2
(gdb) bt
#0 particles::particle_set_momentum () at main_ut.f90:315
#1 0x0000000000405319 in particles::particle_set_set_momentum_all () at main_ut.f90:455
#2 0x0000000000405a92 in isr_epa_handler_uti::isr_handler_1 () at main_ut.f90:508
#3 0x0000000000405cd7 in main_ut () at main_ut.f90:518
and the code:
module kinds
use, intrinsic :: iso_fortran_env
use, intrinsic :: iso_c_binding
implicit none
private
integer, parameter :: single = 4 ! 1.. 6 ! real32 ! c_float
integer, parameter :: double = 8 ! 7..15 ! real64 ! c_double
public :: single, double
public :: default
integer, parameter :: default = double
end module kinds
module constants
use kinds, only: default
implicit none
private
complex(default), parameter, public :: &
imago = (0._default, 1._default)
real(default), parameter, public :: &
zero = 0.0_default
real(default), parameter, public :: &
eps0 = epsilon (zero)
end module constants
module diagnostics
use kinds, only: default
implicit none
private
public :: pacify
interface pacify
module procedure pacify_real_default
end interface pacify
contains
elemental subroutine pacify_real_default (x, tolerance)
real(default), intent(inout) :: x
real(default), intent(in) :: tolerance
if (abs (x) < tolerance) x = 0._default
end subroutine pacify_real_default
end module diagnostics
module lorentz
use kinds, only: default, double
use constants, only: zero, eps0
use diagnostics
implicit none
private
public :: vector3_t
public :: vector3_null
public :: vector3_moving
public :: vector4_t
public :: vector4_null
public :: vector4_at_rest
public :: vector4_moving
public :: operator(*)
public :: operator(**)
type :: vector3_t
real(default), dimension(3) :: p
end type vector3_t
type :: vector4_t
real(default), dimension(0:3) :: p = &
[zero, zero, zero, zero]
end type vector4_t
type(vector3_t), parameter :: vector3_null = &
vector3_t ([ zero, zero, zero ])
type(vector4_t), parameter :: vector4_null = &
vector4_t ([ zero, zero, zero, zero ])
interface vector3_moving
module procedure vector3_moving_canonical
end interface
interface operator(*)
module procedure prod_integer_vector3, prod_vector3_integer
module procedure prod_real_vector3, prod_vector3_real
end interface
interface operator(*)
module procedure prod_vector3
end interface
interface vector4_moving
module procedure vector4_moving_canonical
end interface
interface operator(*)
module procedure prod_real_vector4, prod_vector4_real
module procedure prod_integer_vector4, prod_vector4_integer
end interface
interface operator(*)
module procedure prod_vector4
end interface
interface operator(**)
module procedure power_vector4
end interface
contains
elemental function vector3_canonical (k) result (p)
type(vector3_t) :: p
integer, intent(in) :: k
p = vector3_null
p%p(k) = 1
end function vector3_canonical
elemental function vector3_moving_canonical (p, k) result(q)
type(vector3_t) :: q
real(default), intent(in) :: p
integer, intent(in) :: k
q = vector3_null
q%p(k) = p
end function vector3_moving_canonical
elemental function prod_real_vector3 (s, p) result (q)
type(vector3_t) :: q
real(default), intent(in) :: s
type(vector3_t), intent(in) :: p
q%p = s * p%p
end function prod_real_vector3
elemental function prod_vector3_real (p, s) result (q)
type(vector3_t) :: q
real(default), intent(in) :: s
type(vector3_t), intent(in) :: p
q%p = s * p%p
end function prod_vector3_real
elemental function prod_integer_vector3 (s, p) result (q)
type(vector3_t) :: q
integer, intent(in) :: s
type(vector3_t), intent(in) :: p
q%p = s * p%p
end function prod_integer_vector3
elemental function prod_vector3_integer (p, s) result (q)
type(vector3_t) :: q
integer, intent(in) :: s
type(vector3_t), intent(in) :: p
q%p = s * p%p
end function prod_vector3_integer
elemental function prod_vector3 (p, q) result (s)
real(default) :: s
type(vector3_t), intent(in) :: p,q
s = dot_product (p%p, q%p)
end function prod_vector3
elemental function vector4_canonical (k) result (p)
type(vector4_t) :: p
integer, intent(in) :: k
p = vector4_null
p%p(k) = 1
end function vector4_canonical
elemental function vector4_at_rest (m) result (p)
type(vector4_t) :: p
real(default), intent(in) :: m
p = vector4_t ([ m, zero, zero, zero ])
end function vector4_at_rest
elemental function vector4_moving_canonical (E, p, k) result (q)
type(vector4_t) :: q
real(default), intent(in) :: E, p
integer, intent(in) :: k
q = vector4_at_rest(E)
q%p(k) = p
end function vector4_moving_canonical
elemental function prod_real_vector4 (s, p) result (q)
type(vector4_t) :: q
real(default), intent(in) :: s
type(vector4_t), intent(in) :: p
q%p = s * p%p
end function prod_real_vector4
elemental function prod_vector4_real (p, s) result (q)
type(vector4_t) :: q
real(default), intent(in) :: s
type(vector4_t), intent(in) :: p
q%p = s * p%p
end function prod_vector4_real
elemental function prod_integer_vector4 (s, p) result (q)
type(vector4_t) :: q
integer, intent(in) :: s
type(vector4_t), intent(in) :: p
q%p = s * p%p
end function prod_integer_vector4
elemental function prod_vector4_integer (p, s) result (q)
type(vector4_t) :: q
integer, intent(in) :: s
type(vector4_t), intent(in) :: p
q%p = s * p%p
end function prod_vector4_integer
elemental function prod_vector4 (p, q) result (s)
real(default) :: s
type(vector4_t), intent(in) :: p,q
s = p%p(0)*q%p(0) - dot_product(p%p(1:), q%p(1:))
end function prod_vector4
elemental function power_vector4 (p, e) result (s)
real(default) :: s
type(vector4_t), intent(in) :: p
integer, intent(in) :: e
s = p * p
if (e /= 2) then
if (mod(e, 2) == 0) then
s = s**(e / 2)
else if (s >= 0) then
s = sqrt(s)**e
else
s = -(sqrt(abs(s))**e)
end if
end if
end function power_vector4
end module lorentz
module particles
use kinds, only: default, double
use lorentz
implicit none
private
public :: particle_t
public :: particle_set_t
integer, parameter, public :: PRT_UNDEFINED = 0
integer, parameter, public :: PRT_BEAM = -9
integer, parameter, public :: PRT_INCOMING = 1
integer, parameter, public :: PRT_OUTGOING = 2
integer, parameter, public :: PRT_VIRTUAL = 4
integer, parameter, public :: PRT_BEAM_REMNANT = 9
type :: particle_t
!private
integer :: status = 0
type(vector4_t) :: p = vector4_null
real(default) :: p2 = 0
type(vector4_t), allocatable :: vertex
real(default), allocatable :: lifetime
integer, dimension(:), allocatable :: parent
integer, dimension(:), allocatable :: child
contains
generic :: init => init_particle
procedure :: init_particle => particle_init_particle
procedure :: reset_status => particle_reset_status
procedure :: set_momentum => particle_set_momentum
procedure :: set_children => particle_set_children
procedure :: set_parents => particle_set_parents
procedure :: get_status => particle_get_status
end type particle_t
type :: particle_set_t
! private !!!
integer :: n_beam = 0
integer :: n_in = 0
integer :: n_vir = 0
integer :: n_out = 0
integer :: n_tot = 0
type(particle_t), dimension(:), allocatable :: prt
contains
generic :: assignment(=) => init_particle_set
generic :: init => init_particle_set
procedure :: init_particle_set => particle_set_init_particle_set
procedure :: basic_init => particle_set_basic_init
procedure :: init_direct => particle_set_init_direct
generic :: set_momentum => set_momentum_single
generic :: set_momentum => set_momentum_indices
generic :: set_momentum => set_momentum_all
procedure :: set_momentum_single => particle_set_set_momentum_single
procedure :: set_momentum_indices => particle_set_set_momentum_indices
procedure :: set_momentum_all => particle_set_set_momentum_all
procedure :: reset_status => particle_set_reset_status
end type particle_set_t
contains
subroutine particle_init_particle (prt_out, prt_in)
class(particle_t), intent(out) :: prt_out
type(particle_t), intent(in) :: prt_in
prt_out%status = prt_in%status
prt_out%p = prt_in%p
prt_out%p2 = prt_in%p2
if (allocated (prt_in%vertex)) &
allocate (prt_out%vertex, source=prt_in%vertex)
if (allocated (prt_in%lifetime)) &
allocate (prt_out%lifetime, source=prt_in%lifetime)
end subroutine particle_init_particle
elemental subroutine particle_reset_status (prt, status)
class(particle_t), intent(inout) :: prt
integer, intent(in) :: status
prt%status = status
end subroutine particle_reset_status
elemental subroutine particle_set_momentum (prt, p, p2, on_shell)
class(particle_t), intent(inout) :: prt
type(vector4_t), intent(in) :: p
real(default), intent(in), optional :: p2
logical, intent(in), optional :: on_shell
prt%p = p
if (present (p2)) then
prt%p2 = p2
else
prt%p2 = p ** 2
end if
end subroutine particle_set_momentum
subroutine particle_set_children (prt, idx)
class(particle_t), intent(inout) :: prt
integer, dimension(:), intent(in) :: idx
if (allocated (prt%child)) deallocate (prt%child)
allocate (prt%child (count (idx /= 0)))
prt%child = pack (idx, idx /= 0)
end subroutine particle_set_children
subroutine particle_set_parents (prt, idx)
class(particle_t), intent(inout) :: prt
integer, dimension(:), intent(in) :: idx
if (allocated (prt%parent)) deallocate (prt%parent)
allocate (prt%parent (count (idx /= 0)))
prt%parent = pack (idx, idx /= 0)
end subroutine particle_set_parents
elemental function particle_get_status (prt) result (status)
integer :: status
class(particle_t), intent(in) :: prt
status = prt%status
end function particle_get_status
subroutine particle_set_init_particle_set (pset_out, pset_in)
class(particle_set_t), intent(out) :: pset_out
type(particle_set_t), intent(in) :: pset_in
integer :: i
pset_out%n_beam = pset_in%n_beam
pset_out%n_in = pset_in%n_in
pset_out%n_vir = pset_in%n_vir
pset_out%n_out = pset_in%n_out
pset_out%n_tot = pset_in%n_tot
if (allocated (pset_in%prt)) then
allocate (pset_out%prt (size (pset_in%prt)))
do i = 1, size (pset_in%prt)
pset_out%prt(i) = pset_in%prt(i)
end do
end if
end subroutine particle_set_init_particle_set
subroutine particle_set_basic_init (particle_set, n_beam, n_in, n_vir, n_out)
class(particle_set_t), intent(out) :: particle_set
integer, intent(in) :: n_beam, n_in, n_vir, n_out
particle_set%n_beam = n_beam
particle_set%n_in = n_in
particle_set%n_vir = n_vir
particle_set%n_out = n_out
particle_set%n_tot = n_beam + n_in + n_vir + n_out
allocate (particle_set%prt (particle_set%n_tot))
end subroutine particle_set_basic_init
subroutine particle_set_init_direct (particle_set, &
n_beam, n_in, n_rem, n_vir, n_out, pdg)
class(particle_set_t), intent(out) :: particle_set
integer, intent(in) :: n_beam
integer, intent(in) :: n_in
integer, intent(in) :: n_rem
integer, intent(in) :: n_vir
integer, intent(in) :: n_out
integer, dimension(:), intent(in) :: pdg
integer :: i, k, n
call particle_set%basic_init (n_beam, n_in, n_rem+n_vir, n_out)
n = 0
call particle_set%prt(n+1:n+n_beam)%reset_status (PRT_BEAM)
do i = n+1, n+n_beam
call particle_set%prt(i)%set_children &
([(k, k=i+n_beam, n+n_beam+n_in+n_rem, 2)])
end do
n = n + n_beam
call particle_set%prt(n+1:n+n_in)%reset_status (PRT_INCOMING)
do i = n+1, n+n_in
if (n_beam > 0) then
call particle_set%prt(i)%set_parents &
([i-n_beam])
end if
call particle_set%prt(i)%set_children &
([(k, k=n+n_in+n_rem+1, n+n_in+n_rem+n_vir+n_out)])
end do
n = n + n_in
call particle_set%prt(n+1:n+n_rem)%reset_status (PRT_BEAM_REMNANT)
do i = n+1, n+n_rem
if (n_beam > 0) then
call particle_set%prt(i)%set_parents &
([i-n_in-n_beam])
end if
end do
n = n + n_rem
call particle_set%prt(n+1:n+n_vir)%reset_status (PRT_VIRTUAL)
do i = n+1, n+n_vir
call particle_set%prt(i)%set_parents &
([(k, k=n-n_rem-n_in+1, n-n_rem)])
end do
n = n + n_vir
call particle_set%prt(n+1:n+n_out)%reset_status (PRT_OUTGOING)
do i = n+1, n+n_out
call particle_set%prt(i)%set_parents &
([(k, k=n-n_vir-n_rem-n_in+1, n-n_vir-n_rem)])
end do
end subroutine particle_set_init_direct
subroutine particle_set_set_momentum_single &
(particle_set, i, p, p2, on_shell)
class(particle_set_t), intent(inout) :: particle_set
integer, intent(in) :: i
type(vector4_t), intent(in) :: p
real(default), intent(in), optional :: p2
logical, intent(in), optional :: on_shell
call particle_set%prt(i)%set_momentum (p, p2, on_shell)
end subroutine particle_set_set_momentum_single
subroutine particle_set_set_momentum_indices &
(particle_set, indices, p, p2, on_shell)
class(particle_set_t), intent(inout) :: particle_set
integer, dimension(:), intent(in) :: indices
type(vector4_t), dimension(:), intent(in) :: p
real(default), dimension(:), intent(in), optional :: p2
logical, intent(in), optional :: on_shell
integer :: i
if (present (p2)) then
do i = 1, size (indices)
call particle_set%prt(indices(i))%set_momentum (p(i), p2(i), on_shell)
end do
else
do i = 1, size (indices)
call particle_set%prt(indices(i))%set_momentum &
(p(i), on_shell=on_shell)
end do
end if
end subroutine particle_set_set_momentum_indices
subroutine particle_set_set_momentum_all (particle_set, p, p2, on_shell)
class(particle_set_t), intent(inout) :: particle_set
type(vector4_t), dimension(:), intent(in) :: p
real(default), dimension(:), intent(in), optional :: p2
logical, intent(in), optional :: on_shell
call particle_set%prt%set_momentum (p, p2, on_shell)
end subroutine particle_set_set_momentum_all
subroutine particle_set_reset_status (particle_set, index, status)
class(particle_set_t), intent(inout) :: particle_set
integer, dimension(:), intent(in) :: index
integer, intent(in) :: status
integer :: i
if (allocated (particle_set%prt)) then
do i = 1, size (index)
call particle_set%prt(index(i))%reset_status (status)
end do
end if
particle_set%n_beam = &
count (particle_set%prt%get_status () == PRT_BEAM)
particle_set%n_in = &
count (particle_set%prt%get_status () == PRT_INCOMING)
particle_set%n_out = &
count (particle_set%prt%get_status () == PRT_OUTGOING)
particle_set%n_vir = particle_set%n_tot &
- particle_set%n_beam - particle_set%n_in - particle_set%n_out
end subroutine particle_set_reset_status
end module particles
module isr_epa_handler_uti
use kinds, only: default
use lorentz, only: vector4_t, vector4_moving, operator(*)
use particles, only: particle_set_t
implicit none
private
public :: isr_handler_1
contains
subroutine isr_handler_1 (u)
integer, intent(in) :: u
type(particle_set_t) :: pset
type(vector4_t), dimension(8) :: p
real(default) :: sqrts
real(default), dimension(2) :: x, xb
real(default) :: probability
call pset%init_direct (n_beam = 2, n_in = 2, n_rem = 2, n_vir = 0, n_out = 2, &
pdg = [11, -11, 11, -11, 22, 22, 13, -13])
sqrts = 100._default
x = [0.6_default, 0.9_default]
xb= 1 - x
p(1) = vector4_moving (sqrts/2, sqrts/2, 3)
p(2) = vector4_moving (sqrts/2,-sqrts/2, 3)
p(3:4) = x * p(1:2)
p(5:6) = xb * p(1:2)
p(7:8) = p(3:4)
call pset%set_momentum (p, on_shell = .false.)
end subroutine isr_handler_1
end module isr_epa_handler_uti
program main_ut
use isr_epa_handler_uti, only: isr_handler_1
implicit none
call isr_handler_1 (6)
end program main_ut