[Fortran] A piece of code that confuses me

[Matterwave]

Hello all,

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

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:

  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.

 

[Mark44]

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.

 

[Matterwave]

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.

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

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.

 

[Mark44]

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

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.

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.

 

[Matterwave]

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.

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)?

 

[mathman]

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?

 

[Mark44]

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)?

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."

 

[Matterwave]

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."

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:

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?

 

[Matterwave]

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?

I believe they are in the argument list for the subroutine...but I'm not sure where you're going with this?

 

[Mark44]

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?

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.

In which case, how does it know how to put a 6-D array into a 2-D array correctly?

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).

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?

 

[Matterwave]

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).

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

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))"?

 

[Mark44]

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

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))"?

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.

 

[FactChecker]

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.

 

[Matterwave]

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.

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?

 

[jim mcnamara]

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.

 

[Matterwave]

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.

It's actually a defined parameter that's in a separate module for this code. It's not a Fortran standard.

 

[Mark44]

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?

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.

 

[Matterwave]

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.

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?

 

[FactChecker]

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.

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.

 

[jtbell]

Any variable in a function call parameter list represents the address of that variable.

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.

 

[mathman]

I believe they are in the argument list for the subroutine...but I'm not sure where you're going with this?

"subroutine onestep(iter,iegy,icos,jcos,iphi,rinit,rfin, &
& v2,v2bar,hel,helbar,drdis, &
& psib,psie,psibarb,psibare,err)"

I don't see them (nflavor,ncosth) here. In the subroutine they are just parameters - specific values need to be assigned by the call.

 

[Matterwave]

"subroutine onestep(iter,iegy,icos,jcos,iphi,rinit,rfin, &
& v2,v2bar,hel,helbar,drdis, &
& psib,psie,psibarb,psibare,err)"

I don't see them (nflavor,ncosth) here. In the subroutine they are just parameters - specific values need to be assigned by the call.

Oh, I think that's because nflavor and ncosth are in fact parameters defined in a different module, their value is already assigned at compile time.