Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

[Fortran] A piece of code that confuses me

  1. Sep 22, 2014 #1

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    Hello all,

    There is a subroutine in the code I'm working on that looks like this:

    Code (Text):

    subroutine onestep(iter,iegy,icos,jcos,iphi,rinit,rfin, &
      &  v2,v2bar,hel,helbar,drdis, &
      &  psib,psie,psibarb,psibare,err)
      use params
      use helectron
      use angles
      use neutrinodist
      use parallel
      implicit real(kind=reel8) (a-h,o-z)
      integer,intent(in) :: iegy,icos,jcos
      real(kind=reel16),intent(in) :: rinit,rfin
      real(kind=reel16),dimension(ncosth),intent(in) :: drdis
      real(kind=reel8),dimension(nflavor,nflavor,negy,nphi),intent(IN) :: hel,helbar
      complex(kind=reel8),dimension(nflavor,nflavor,ncosth),intent(in)::v2,v2bar
      complex(kind=reel8),dimension(nflavor,nflavor) :: v,exphel,vbar,exphbarel
      complex(kind=reel8),dimension(nflavor,nflavor):: psib,psie,psibarb,psibare
      complex(kind=reel8),dimension(nflavor):: psio,psibaro
      complex(kind=reel8) :: ph1,ph2,csldr
      real(kind=reel8) :: xl,xldr,xli
      real(kind=reel8) :: cldr,sldr,xbs,xas
      real(kind=reel8),intent(out),dimension(nflavor,2) :: err
      real(kind=reel8) :: xa,xb,xc
      complex(kind=reel8) :: hnow(nflavor,nflavor)

    etc.
     
    So there are some inputs psib, psie, etc., which are arrays as declared in the subroutine. However, when the code goes to call this subroutine, the call command looks like this:

    Code (Text):

      call onestep(iter,iegy,icos,jcos,iphi,rinit,rfin,v2,v2bar,hel,helbar,drdis, &
      &  psi(1,1,iegy,jcos,iphi,ibeg),  &
      &  psi(1,1,iegy,jcos,iphi,iend),  &
      &  psibar(1,1,iegy,jcos,iphi,ibeg), &
      &  psibar(1,1,iegy,jcos,iphi,iend), &
      &  err(1,1,iegy,jcos,iphi))
     
    All the iegy, jcos, etc. are integer counters for nested do-loops.

    This seems that they passed one single complex number psi(1,1,iegy,jcos,iphi,ibeg) into a slot in the subroutine designated for a matrix psib which is an nflavor by nflavor (2-d) matrix.

    Can anyone guess what's going on? I would have assumed that the correct way to pass the matrix into the subroutine would be to use something like psi( :,:,iegy,jcos,iphi,ibeg) so that it's passing the first 2 slots into the 2 dimensional array psib.
     
  2. jcsd
  3. Sep 23, 2014 #2

    Mark44

    Staff: Mentor

    You haven't given us enough information to answer your question. We see the declarations for psib, psie, psibarb, psibare and a few others, but you're asking about psi( ... ) and psibar( ... ), neither of which appears in the code you show. I would guess that these are functions, but without seeing how they are defined, that's only a guess.
     
  4. Sep 23, 2014 #3

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    Sorry, psi and psibar are 6 dimensional arrays. They are declared here:

    Code (Text):

    complex(kind=reel8),dimension(nflavor,nflavor,negy,ncospproc,nphi,2):: psi
    complex(kind=reel8),dimension(nflavor,nflavor,negy,ncospproc,nphi,2):: &
      &  psibar
     
    The first 2 dimensions of psi and psibar match psib and psie (which are psi beginning and psi end) and the last dimension in psi and psibar is storing both the beginning and ending psi's, the negy, ncospproc, and nphi are to store all the energy and angle bins. I understand onestep works by evolving one specific energy/angle bin, but I don't understand why they pass only psi(1,1,etc) to it instead of psi( :,:,etc) to it.
     
  5. Sep 23, 2014 #4

    Mark44

    Staff: Mentor

    The first two dimensions apparently give the size of a square matrix, the smallest of which would be 1 X 1, which is essentially a scalar. It seems far superior to pass values in the call to psi( ... ) than to omit the first two parameters, which would likely cause the compiler to generate an error.

    When you use a function, you need to supply as many parameters of the appropriate types as are listed in the function's header.
     
  6. Sep 23, 2014 #5

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    I'm not sure I understand. In the subroutine, I am passing psi -> psib and psie, but psib and psie are nflavor by nflavor arrays, shouldn't I pass an array of dimension nflavor by nflavor to the subroutine instead of 1 number like psi(1,1,1,1,1,1)?
     
  7. Sep 23, 2014 #6

    mathman

    User Avatar
    Science Advisor
    Gold Member

    I haven't worked in Fortran for a while so I may be missing something. However, shouldn't nflavor, ncosth, etc. be part of the argument list for the subroutine?
     
  8. Sep 23, 2014 #7

    Mark44

    Staff: Mentor

    OK, I think I finally understand what you're asking.

    In the definition of the onestep subroutine, the 13th argument is psib, which is declared in the sub as a two-dimensional array of type complex.

    However, in the call to onestep, the actual parameter is psi(...), which is defined as a six-dimensional array of type complex.

    It's been years since I did any Fortran programming, but I suspect is happening is that what is actually passed is the address in memory of the beginning of the array, whether 2-D or 6-D, and that the elements in the array are used by the called routine to figure out the actual shape of the array. Those details would probably be evident in the subroutine implementation, where you show "etc."
     
  9. Sep 23, 2014 #8

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    Where I show "etc.", the subroutine treats psi and psib as 2-D arrays and evolves them forward in time, there's not much else happening. I can copy paste that part of the code:

    Code (Text):

    subroutine onestep(iter,iegy,icos,jcos,iphi,rinit,rfin, &
      &  v2,v2bar,hel,helbar,drdis, &
      &  psib,psie,psibarb,psibare,err)
      use params
      use helectron
      use angles
      use neutrinodist
      use parallel
      implicit real(kind=reel8) (a-h,o-z)
      integer,intent(in) :: iegy,icos,jcos
      real(kind=reel16),intent(in) :: rinit,rfin
      real(kind=reel16),dimension(ncosth),intent(in) :: drdis
      real(kind=reel8),dimension(nflavor,nflavor,negy,nphi),intent(IN) :: hel,helbar
      complex(kind=reel8),dimension(nflavor,nflavor,ncosth),intent(in)::v2,v2bar
      complex(kind=reel8),dimension(nflavor,nflavor) :: v,exphel,vbar,exphbarel
      complex(kind=reel8),dimension(nflavor,nflavor):: psib,psie,psibarb,psibare
      complex(kind=reel8),dimension(nflavor):: psio,psibaro
      complex(kind=reel8) :: ph1,ph2,csldr
      real(kind=reel8) :: xl,xldr,xli
      real(kind=reel8) :: cldr,sldr,xbs,xas
      real(kind=reel8),intent(out),dimension(nflavor,2) :: err
      real(kind=reel8) :: xa,xb,xc
      complex(kind=reel8) :: hnow(nflavor,nflavor)

    !  the full Hamiltonian is Hermitian and of the form:
    !
    !  hel(1,1)+v2(1,1,icos)  hel(1,2)+v2(1,2,icos)
    !  hel(2,1)+v2(2,1,icos)  hel(2,2)+v2(2,2,icos)
    !
    !  a  b+ic
    !  b-ic  -a
    !
    !  the diagonal terms are real, off-diagonal are complex
    !
    !  diagonalize neutrino matrix

      hnow(:,:)=hel(:,:,iegy,iphi)+v2(:,:,icos)
       
      call exphcal( hnow,drdis(icos),exphel)

      psio(:)=psie(:,1)
      do ifli=1,nflavor

      psie(:,ifli)=czero

      do i=1,nflavor
      do j=1,nflavor
      psie(i,ifli)=psie(i,ifli)+exphel(i,j)*psib(j,ifli)
      enddo
      enddo

      enddo
      do i=1,nflavor
      err(i,1)=1.d2*( (dble(psie(i,1))-dble(psio(i)))**2  &
      &  + (imag(psie(i,1))-imag(psio(i)))**2)
      enddo
       

    !  diagonalize anti-neutrino matrix

      hnow(:,:)=helbar(:,:,iegy,iphi)+v2bar(:,:,icos)
      call exphcal( hnow,drdis(icos),exphel)

      psibaro(:)=psibare(:,1)
      do ifli=1,nflavor

      psibare(:,ifli)=czero
      do i=1,nflavor
      do j=1,nflavor
      psibare(i,ifli)=psibare(i,ifli)+exphel(i,j)*psibarb(j,ifli)
      enddo
      enddo

      enddo
      do i=1,nflavor
      err(i,2)=1.d2*((dble(psibare(i,1))-dble(psibaro(i)))**2 &
      &  +(imag(psibare(i,1))-imag(psibaro(i)))**2)
      enddo




      if(debug.gt.5) then
      if(iter.eq.3.and.iegy.eq.1.and.icos.eq.1) then
      write(6,*) '--------------------------'
      write(6,*) ' End of Onestep'
      write(6,*) ' rinit,rfin = ',rinit,rfin
      write(6,*) ' psio = ',psio
      write(6,*) ' psiele = ',psie
      write(6,*) ' psielb = ',psib
      write(6,*) ' v = ',v
      write(6,*) ' v2 = ',v2(:,:,1)
      write(6,*) ' xl ',xl
      write(6,*) ' xldr ',xldr
      write(6,*) ' coa ',coa
      write(6,*) '--------------------------'
      endif
      endif

      return
      end
     
    exphcal is a subroutine that finds the eigenvalues and eigenvectors of the Hamiltonian (and evolves with ##e^{iH\Delta t}## by one step).

    I don't see anywhere in this subroutine that is talking about finding the memory spaces and such? I don't even know how that kind of code would look like. @_@ Maybe Fortran is automatic in doing that? In which case, how does it know how to put a 6-D array into a 2-D array correctly? I.e. how does it know that it is the first 2 dimensions of the 6-D array which must fit into the 2-D array, keeping the other 4 dimensions of the 6-D array constant?
     
  10. Sep 23, 2014 #9

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    I believe they are in the argument list for the subroutine...but I'm not sure where you're going with this?
     
  11. Sep 23, 2014 #10

    Mark44

    Staff: Mentor

    AFAIK, yes. In the call to onestep, the 13th argument isn't actually an array - it's the address in memory of the start of the array. I'm reasonably sure that that's how Fortran does things.
    Again, and I think this is true, it's not actually working with 6-D arrays or 2-D arrays. The thing that actually gets passed is, I believe, a single number, the address in memory where the array starts. The code for onestep is able to pull out the first two numbers (nflavor and nflavor) to figure out the shape of the array and go from there.

    BTW, I wouldn't have coded things this way. For one thing, it's Gawdawful style to write a subroutine with 16 (or whatever) parameters. It would be better to package them into a single structure and pass that.

    For another thing, all the arrays involved here are really 2-D arrays (I think) and the other stuff is just cruft that was needed because the original author didn't seem to know about structures (STRUCTURE keyword) or records (RECORD keyword).

     
  12. Sep 23, 2014 #11

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    Interesting! This would explain it. In that case, I must ask, do you know if the passing of the memory address thing is the only way to pass an array in a Fortran code? If I wrote something like

    Code (Text):

    program test
    real,dimension(3,3)::r
    call setzero(r)
    end program test

    contains

    subroutine setzero(a)
    real,dimension(3,3)::a
    a=0
    return
    end subroutine setzero

     
    would that work, or do I have to wrote something like "call setzero(r(1,1))"?
     
  13. Sep 23, 2014 #12

    Mark44

    Staff: Mentor

    The first is the way to go. In the line just above you are passing r(1, 1) which is a single element in the array.
     
  14. Sep 23, 2014 #13

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    In old Fortran, all parameters were "pass by address". Then some computers/compilers started using a little "pass by value" for speed in special cases. I don't know what the current situation is. For large matrices, it is clearly faster to just pass one address than to copy the values of the entire array.
     
  15. Sep 23, 2014 #14

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    But isn't passing r(1,1) the same as what the other code was doing by passing psi(1,1,iegy,etc)? I thought by writing r(1,1) in the call command, it would pass the address to the first element?
     
  16. Sep 23, 2014 #15

    jim mcnamara

    User Avatar

    Staff: Mentor

    reel8 - Does your compiler allow that datatype? - used to be real*8, so I suppose it is a kind of typedef. Or a typo. Have not done FORTRAN since F77. So I may be next door to useless.
     
  17. Sep 23, 2014 #16

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    It's actually a defined parameter that's in a separate module for this code. It's not a Fortran standard.
     
  18. Sep 23, 2014 #17

    Mark44

    Staff: Mentor

    No, r(1, 1) is the element in row 1 column 1, which is not at all the same as its address. r all by itself represents the whole array, and is most likely the address of the first element in the array.
     
  19. Sep 24, 2014 #18

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    So I don't understand why r(1,1) is the element in row 1 column 1, but psi(1,1,iegy,jcos,iphi,ibeg) is the address of the first element?
     
  20. Sep 24, 2014 #19

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    In FORTRAN, that depends on where r(1,1) appears. If it is a parameter in a function call, it represents the address of that element of the array. Any variable in a function call parameter list represents the address of that variable. Elsewhere, it is the value of the variable.
     
  21. Sep 24, 2014 #20

    jtbell

    User Avatar

    Staff: Mentor

    This is true even for single non-array variables. A type of bug that I ran into regularly when I programmed in Fortran 77 was to accidentally try to pass a single variable to a subprogram that was expecting an array. When the subprogram modified the contents of the array, it modified not only the corresponding variable in the calling (sub)program, it also modified other variables that happened to be allocated following it in memory! This usually produced run-time error messages referring to sections of the program different from the location of the subprogram call, making the bug difficult to track down.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: [Fortran] A piece of code that confuses me
  1. Help with fortran code (Replies: 17)

  2. Help with fortran code (Replies: 4)

Loading...