# [Fortran] A piece of code that confuses me

1. Sep 22, 2014

### Matterwave

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. Sep 23, 2014

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

3. Sep 23, 2014

### Matterwave

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.

4. Sep 23, 2014

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

5. Sep 23, 2014

### Matterwave

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

6. Sep 23, 2014

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

7. Sep 23, 2014

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

8. Sep 23, 2014

### Matterwave

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?

9. Sep 23, 2014

### Matterwave

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

10. Sep 23, 2014

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

11. Sep 23, 2014

### Matterwave

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

12. Sep 23, 2014

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

13. Sep 23, 2014

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

14. Sep 23, 2014

### Matterwave

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?

15. Sep 23, 2014

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

16. Sep 23, 2014

### Matterwave

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

17. Sep 23, 2014

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

18. Sep 24, 2014

### Matterwave

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?

19. Sep 24, 2014

### FactChecker

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.

20. Sep 24, 2014

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