From: robert.corbett on
On Feb 16, 2:14 pm, Paul van Delst <paul.vande...(a)noaa.gov> wrote:
> glen herrmannsfeldt wrote:
> > Paul van Delst <paul.vande...(a)noaa.gov> wrote:
> > (snip)
>
> >>  ELEMENTAL FUNCTION Cloud_Add( cld1, cld2 ) RESULT( cldsum )
> > (snip)
> >> that is overloaded with OPERATOR(+),
> > (snip)
>
> >>    IF ( atm1%n_Clouds > 0 ) THEN
> >>      nc = atm1%n_Clouds
> >>      atmsum%Cloud(1:nc) = atmsum%Cloud(1:nc) + atm2%Cloud(1:nc)
> >>    END IF
> > (snip)
> >> When invoked, the line
>
> >>  atmsum%Cloud(1:nc) = atmsum%Cloud(1:nc) + atm2%Cloud(1:nc)
>
> >> causes a segfault. When I change it to the more usual
>
> >>  atmsum%Cloud = atmsum%Cloud + atm2%Cloud
>
> >> everything is fine. No segfault and results are correct.
>
> > What is the value of nc at that point?  If it is nowhere near
> > the value it should have, then segfault is likely.  If it is,
> > then I would wonder about compiler bugs.
>
> > I notice that nc comes from atm1, but atm2 is used in the sum.
> > I would print nc as a debugging aid until the problem was fixed.
>
> Well, I can't print it out since the function is elemental,

You can if you use a technique I call "lying." Declare the function
to be elemental
in the code that calls it, but do not declare it elemental in the code
that defines
it. You might need to break the function out into a file of its own
to compile it.
Of course, you should do this only for debugging, never for production
code.

Bob Corbett
From: Paul van Delst on
Paul van Delst wrote:
> Hello,
>
> I have a module function like so
[description deleted]

Just to follow up:

Down below is a (relatively) simple working example that indicates the issue.

If I compile and run through the debugger I get:

$ xlf2003 -qversion
IBM XL Fortran for AIX, V12.1
Version: 12.01.0000.0001

$ xlf2003 -qdbg test_elementalslice.f90
** b_define === End of Compilation 1 ===
** a_define === End of Compilation 2 ===
** test_elementalslice === End of Compilation 3 ===
1501-510 Compilation successful for file test_elementalslice.f90.

$ dbx a.out
Type 'help' for help.
reading symbolic information ...
(dbx) run
n_a = 1
a%n = 4
a%x = 3.141590118 3.141590118 3.141590118 3.141590118
n_b = 1
b%n = 4
b%x = 0.6931470037 0.6931470037 0.6931470037 0.6931470037
n_b = 2
b%n = 4
b%x = 0.6931470037 0.6931470037 0.6931470037 0.6931470037
n_b = 3
b%n = 4
b%x = 0.6931470037 0.6931470037 0.6931470037 0.6931470037
n_a = 2
a%n = 4
a%x = 3.141590118 3.141590118 3.141590118 3.141590118
n_b = 1
b%n = 4
b%x = 0.6931470037 0.6931470037 0.6931470037 0.6931470037
n_b = 2
b%n = 4
b%x = 0.6931470037 0.6931470037 0.6931470037 0.6931470037
n_b = 3
b%n = 4
b%x = 0.6931470037 0.6931470037 0.6931470037 0.6931470037

Segmentation fault in . at 0xf458
0x000000000000f458 f8e30008 std r7,0x8(r3)
(dbx) up
__a_define_NMOD_a_add(??, ??, ??), line 126 in "test_elementalslice.f90"
(dbx) quit


The line in question is
asum%b(1:n_b) = asum%b(1:n_b) + a2%b(1:n_b)

The next commented out line
! asum%b = asum%b + a2%b

gets around the problem, but this doesn't allow for the structure arrays to be allocated
to a larger dimension, and then used for a smaller "n_b" -- for example, if I want to
allocate the "b" array component of "a" to some maximum size and then only operate on a
subset.

I've submitted the test case to our local IBM folks for their assessment.

cheers,

paulv

<-----begin----->
module b_define
implicit none

private
public :: operator(+)
public :: b_type
public :: b_destroy
public :: b_create
public :: b_inspect

interface operator(+)
module procedure b_add
end interface operator(+)

type :: b_type
integer :: n
real, allocatable :: x(:)
end type

contains

! Public

elemental subroutine b_destroy(b)
type(b_type), intent(out) :: b
end subroutine b_destroy

elemental subroutine b_create(b,n)
type(b_type), intent(out) :: b
integer, intent(in) :: n
allocate(b%x(n))
b%n = n
b%x = 0.0
end subroutine b_create

subroutine b_inspect(b)
type(b_type), intent(in) :: b
print *, ' b%n = ',b%n
print *, ' b%x = ',b%x
end subroutine b_inspect

! Private

elemental function b_add(b1,b2) result(bsum)
type(b_type), intent(in) :: b1, b2
type(b_type) :: bsum
integer :: n
if ( b1%n /= b2%n ) return
bsum = b1
n = b1%n
bsum%x(1:n) = bsum%x(1:n) + b2%x(1:n)
end function b_add

end module b_define



module a_define
use b_define
implicit none

private
public :: operator(+)
public :: a_type
public :: a_destroy
public :: a_create
public :: a_inspect


interface operator(+)
module procedure a_add
end interface operator(+)

type :: a_type
integer :: n, n_b
real, allocatable :: x(:)
type(b_type), allocatable :: b(:)
end type

contains

! Public

elemental subroutine a_destroy(a)
type(a_type), intent(out) :: a
end subroutine a_destroy

elemental subroutine a_create(a,n,n_b)
type(a_type), intent(out) :: a
integer, intent(in) :: n, n_b
allocate(a%x(n))
if ( n_b > 0 ) then
allocate(a%b(n_b))
call b_create(a%b,n)
end if
a%n = n
a%n_b = n_b
a%x = 0.0
end subroutine a_create

subroutine a_inspect(a)
type(a_type), intent(in) :: a
integer :: i
print *, 'a%n = ',a%n
print *, 'a%x = ',a%x
if ( a%n_b > 0 ) then
do i = 1, a%n_b
print *, ' n_b = ',i
call b_inspect(a%b(i))
end do
end if
end subroutine a_inspect

! Private

elemental function a_add(a1,a2) result(asum)
type(a_type), intent(in) :: a1, a2
type(a_type) :: asum
integer :: n, n_b
if ( a1%n /= a2%n ) return
asum = a1
n = a1%n
asum%x(1:n) = asum%x(1:n) + a2%x(1:n)
if ( a1%n_b > 0 ) then
n_b = a1%n_b
asum%b(1:n_b) = asum%b(1:n_b) + a2%b(1:n_b)
! asum%b = asum%b + a2%b
end if
end function a_add

end module a_define



program test_elementalslice
use a_define
integer, parameter :: n = 4
integer, parameter :: n_a = 2, n_b = 3
integer :: i, j
type(a_type) :: a(n_a), a2(n_a)
call a_create(a,n,n_b)
do i = 1, n_a
a(i)%x = 3.14159
do j = 1, a(i)%n_b
a(i)%b(j)%x = 0.693147
end do
print *, 'n_a = ',i
call a_inspect(a(i))
end do
a2 = a
a = a + a2
print *, '*** After summation ***'
do i = 1, n_a
print *, 'n_a = ',i
call a_inspect(a(i))
end do
call a_destroy(a)
call a_destroy(a2)
end program test_elementalslice
<-----end----->
From: glen herrmannsfeldt on
Paul van Delst <paul.vandelst(a)noaa.gov> wrote:
(snip)

> Down below is a (relatively) simple working example that
> indicates the issue.

> If I compile and run through the debugger I get:

> $ xlf2003 -qversion
> IBM XL Fortran for AIX, V12.1
> Version: 12.01.0000.0001

(snip)

> Segmentation fault in . at 0xf458
> 0x000000000000f458 f8e30008 std r7,0x8(r3)
> (dbx) up
> __a_define_NMOD_a_add(??, ??, ??), line 126 in "test_elementalslice.f90"
> (dbx) quit

I would also print out the register contents at this point.

It is not so easy to tell sometimes, but if r3 has a really
strange value then sometimes it is obvious. A few instructions
before and after, which I presume would be a loop, would also help.

I had thought we were discussing IA32, as that is what most
people use. It seems, though, to be PowerPC. That likely means
that the compiler hasn't been tested as much.

> The line in question is
> asum%b(1:n_b) = asum%b(1:n_b) + a2%b(1:n_b)

> The next commented out line
> ! asum%b = asum%b + a2%b

You could also compile just the one function with -S and
post the generated assembly code. Most likely, though, IBM
will find it.

(snip of the whole program listing)

-- glen
From: Richard Maine on
glen herrmannsfeldt <gah(a)ugcs.caltech.edu> wrote:

> I had thought we were discussing IA32, as that is what most
> people use.

For XLF?

> It seems, though, to be PowerPC. That likely means
> that the compiler hasn't been tested as much.

That seems an odd conclusion. XLF has a history going back... gee... I
don't recall how far, but well longer than some of the IA32 compilers
have existed. I don't think that the popularity of the platform is a
very good measure of compiler robustness. There have sure been plenty of
junk compilers on popular platforms. I think I'll avoid the distraction
of mentioning compiler names.

--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
From: Paul van Delst on
glen herrmannsfeldt wrote:
> I had thought we were discussing IA32, as that is what most
> people use. It seems, though, to be PowerPC. That likely means
> that the compiler hasn't been tested as much.

No, xlf. I didn't point it out in my original post, but I did in a followup.
Linux+gfortran works fine. AIX+xlf does not.

>> The line in question is
>> asum%b(1:n_b) = asum%b(1:n_b) + a2%b(1:n_b)
>
>> The next commented out line
>> ! asum%b = asum%b + a2%b
>
> You could also compile just the one function with -S and
> post the generated assembly code. Most likely, though, IBM
> will find it.

It is being submitted to the IBM compiler people. I was told it will likely be a low
priority since there is a simple workaround, i.e. do something like
do i = 1, n_b
asum%b(i) = asum%b(i) + a2%b(i)
end do
rather than
asum%b(1:n_b) = asum%b(1:n_b) + a2%b(1:n_b)

cheers,

paulv