elemental function question regarding use of array slices
in [Fortran]
|
Prev: How to control/prevent cast between from complex to realvariables
Next: Variable scope: how fortran deals with undeclared variables in subroutines?
From: robert.corbett on 17 Feb 2010 00:25 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 17 Feb 2010 09:18 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 17 Feb 2010 12:25 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 17 Feb 2010 12:46 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 17 Feb 2010 13:36
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 |