1

I need your help. I'm rewriting a python program in fortran because of shorter calculation time.

There's one problem I really don't know what's wrong here. I want to use the arrays r,u,v in the main program, which I calculate in the subroutine.

Because the arrays change size, they have to be dynamic.

The following program has no errors in the building, but as soon I run it, it gets the error:

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Here's the main program: (u and v are not yet written, but should work analogous to r)

program single_sol2
implicit none

real :: dr =0.01
real :: chi0 = 1.0
integer ::order =1
real :: w=1
real :: phi0 = -2
real, allocatable :: r(:)

call iteration2(chi0, phi0, w, dr, order, r)
!print*,r

end program single_sol2



!And here the subroutine:


subroutine iteration2(chi0, phi0, w, dr, order, r)
implicit none

real :: dr
real :: chi0
integer :: order
real :: w
real :: phi0

real,intent(out):: r(:)
!real, dimension (:), allocatable:: r
real, dimension (:), allocatable :: u
real, dimension (:), allocatable :: v
real, parameter ::  pi   = 3.141593
integer ::  start, ende, rate, cmax
real :: time

integer :: i, solve_lim, j, k, l
real :: u_j, v_j, u_cond, rmax
real :: prec=0
logical :: behind_max
character (len=:), allocatable :: ToDo
character (len=:), allocatable :: last
integer :: max_iter =1000

call system_clock(start,rate,cmax)

do i=1, max_iter
    k=order-1
    behind_max = .false.
    rmax=500
    r=(/0.0,dr/)
    u=(/0.0,chi0*dr/)
    v=(/0.0,phi0*dr/)
    solve_lim = int(rmax/dr) +100
    !print*,i

    do j=1,solve_lim
       r=(/r,(j+1)*dr/)
       u_j= 2 * dr**2 * u(j+1) * (w + v(j+1) / r(j+1) - u(j+1)**2 / (8 * r(j+1)**2)) + 2 * u(j+1) - 
           u(j)
       u=(/u,u_j/)
       !print*,u
       v_j = 4 * pi * dr**2 * u(j+1)**2 / r(j+1) + 2 * v(j+1) - v(j)
       v=(/v,v_j/)


       if (k==0 .and. u(j+2) * u(j+1) < 0) then
                !print*,'firstif'
                ToDo = 'up'  ! crossing zero the "order"-th time
                !print*,'firstif'
                exit
       else
                !print*,'notfirstif'
       end if

       u_cond= u(j+1) * (u(j+2) - 2 * u(j+1) + u(j))
       !print*,u_cond
       !print*,abs(u(j+2))
       !print*,maxval(abs(u))

       if ((behind_max .and. abs(u(j+2)) > abs(u(j+1))) .or. (abs(u(j+2)) == maxval(abs(u)) .and. 
            u_cond > 0 )  ) then
                ToDo = 'down'  ! curving away from axis
                !print*, 'secondif'
                exit
       else
                !print*, 'notsecondif'
       end if


       if (.not. behind_max .and. k==0 .and. abs(u(j+2)) < abs(u(j+1)) ) then
                behind_max = .True.
                rmax = 10 * r(j+2) / order ! detecting last maximum
                !print*,rmax
                !print*, 'thirdif'
       else
                !print*,'notthirdif'
       end if

       if (k /= 0 .and. u(j+2) * u(j+1) < 0) then
                k = k -1  !detecting zero-crossing
       end if

       if (r(j+2) > rmax) then
                ToDo = 'exit'  !detecting , when finished
                exit
       end if
    end do

    !print*,u
    !print*,rmax

    if (ToDo == 'up') then
       if (last == 'down' .or. phi0 > -1.1 * 10**(-prec)) then
                phi0  = phi0 + 0.5 * 10**(-prec)
                prec = prec + 1  ! found current digit -> higher prec.
                last = 'higher prec'

            else
                phi0 = phi0 + 10**(-prec)
                last = 'up'  ! do what 'todo'
        end if

    elseif (ToDo == 'down') then
       if (last == 'up') then
                phi0 = phi0  -0.5 * 10**(-prec)
                prec = prec + 1  ! found current digit -> higher prec.
                last = 'higher prec'

        else
                phi0 = phi0 -10**(-prec)
                last = 'down'  !do what 'todo'
        end if
     end if

    !print*, prec

     if (prec == 7 .or. ToDo == 'exit') then
            exit
     end if
      ! finishing conditions
!print*,u
end do

open(10,file = 'D:\Masterarbeit\Arraytxtfiles\v_array_fortran.txt')
do l=1,size(v)
write(10,*) v(l)
end do

open(10,file = 'D:\Masterarbeit\Arraytxtfiles\r_array_fortran.txt')
do l=1,size(r)
write(10,*) r(l)
end do

open(10,file = 'D:\Masterarbeit\Arraytxtfiles\u_array_fortran.txt')
do l=1,size(u)
write(10,*) u(l)
end do
!s=size(u)
!print*,v
!print*,r(1000)
call system_clock(ende)
time=float(ende-start)/float(rate)   ! evtl. cmax beachten!
write(*,*) "Zeit in Sekunden: ",time
!print*,r
!print*,
return
end subroutine iteration2

Here is the working solution:

    program test
  implicit none
   real,dimension(:),allocatable :: myarray

   interface
      subroutine iteration (r)
      real, dimension(:),allocatable :: r
      end subroutine iteration
   end interface

   call iteration(myarray)
   print*,myarray

end program test

subroutine iteration(r)
implicit none
real,dimension(:),allocatable :: r

   real :: dr =0.01
   real :: chi0 = 1.0
   integer ::order =1
   real :: w=1
   real :: phi0 = -2

real, dimension (:), allocatable:: u
real, dimension (:), allocatable:: v
real, parameter ::  pi   = 3.141593
integer ::  start, ende, rate, cmax
real :: time

integer :: i, solve_lim, j, k, l
real :: u_j, v_j, u_cond, rmax
real :: prec=0
logical :: behind_max
character (len=:), allocatable :: ToDo
character (len=:), allocatable :: last
integer :: max_iter =1000

call system_clock(start,rate,cmax)

do i=1, max_iter
    k=order-1
    behind_max = .false.
    rmax=500
    r=(/0.0,dr/)
    u=(/0.0,chi0*dr/)
    v=(/0.0,phi0*dr/)
    solve_lim = int(rmax/dr) +100
    !print*,i

     .......
 
return
end subroutine iteration
fortran1
  • 19
  • 3
  • On very first glance I am tempted to suggest https://stackoverflow.com/questions/65393812/dynamic-array-and-allocation-in-subroutine as a duplicate. However I suspect there are more problems than just that and don't have time at the moment to dig deeper. – Ian Bush Dec 21 '20 at 14:24
  • Please compile with all debugging and checking options enabled and look at fixing issues that are reported. As the previous comment suggests, there is more than one problem with your code, so to help us focus our answer please try to explain exactly what you are wanting to do, and consider creating a much simpler program (see [mre]). For example, you have commented out another declaration for the argument `r` which potentially distracts from what you ultimately want to address. – francescalus Dec 21 '20 at 14:30
  • 1
    Please do not try the solution given in the answer you have. This will not fix your problem but you may be unfortunate to remove the segmentation fault and then have to spend more time working out why you are getting rubbish output. Instead, look at allocation inside the subroutine, as the linked question considers – francescalus Dec 21 '20 at 15:05
  • First, thank you for the quick comments. But they are not helpful at the moment. The other post I already saw, but i don't see how it should help me here. I can't do a simpler program, because of the time it takes. I only can say, that the subroutine as a program itself works fine(the values for r,u,v... are as expected). I want use the arrays r,u,v and the number phi0 for further calculation in the main program. This works fine for phi0, but it is a real number and not an array. The debugger says: suspendet breakpoint, but I don't know what it means here. – fortran1 Dec 21 '20 at 15:32
  • If you make `r` in the subroutine allocatable (the commented out declaration) and provide an explicit interface (preferably using a module) in the main program, what else goes wrong? (Details of both of those aspects are in answers to that question.) – francescalus Dec 21 '20 at 15:53
  • I solved it. first I had two source files. I really didn't thought it makes no difference. May someone can say what to do with 2 source files. And the other things: look at the code. – fortran1 Dec 21 '20 at 17:17
  • A better way, as francescalus suggests, is modules; you really shouldn't need to write interface blocks like this except for some corner cases, and modules removes this need. And multiple files should make no difference. – Ian Bush Dec 21 '20 at 18:18
  • Ah thx. Do you have an example how this could work? Because I tried a little but i didn't work. – fortran1 Dec 21 '20 at 23:56

0 Answers0