Nested function calls in Fortran?

1

I’m writing a program that involves the Baker-Campbell-Hausdorff series. Suppose that I want to evaluate expressions with multiple nested commutators. (My actual application involves multiple nested Poisson brackets of 1D arrays, but it’s the same idea.) Suppose I’ve written a function cmt(a,b) that takes two matrices as input and returns their commutator,

cmt(a,b) = matmul(a,b) - matmul(b,a)

I need to evaluate very long expressions produced by Mathematica that look like,

e=cmt(a,b) + cmt(a,cmt(b,c))/3. - cmt(cmt(a,b),cmt(c,d))/10.

Is this acceptable in Fortran? gfortran doesn’t complain, but I’m getting strange results. I’m not sure if I’ve made a coding error or if this is non-standard Fortran. Is it legitimate for the actual argument in a function call (i.e., a user-supplied function, not an intrinsic) to contain the same function name?

1 Like
2

No problem. The function call is replaced by the returned result, and so on.

Can you post more details about the problem? And more code?

3

Maybe this can be useful,

program main
  use iso_fortran_env, only: wp => real64

  implicit none

  integer, parameter :: n=4
  real(kind=wp)      :: R(n), s(n), T(n), U(n)
  real(kind=wp)      :: W(n,n)

  R = [1.0, 2.0, 3.0, 4.0]
  S = [1.0, 2.0, 3.0, 4.0]
  T = [1.0, 2.0, 3.0, 4.0]
  U = [1.0, 2.0, 3.0, 4.0]

  write(*,*) '[R,S]'
  W = cmt(R,S)
  write(*,'(4f10.1)') W

  write(*,*) '[R,S] + [R,[R,S]]'
  W = cmt(R,S) + cmt(R,cmt(R,S))
  write(*,'(4f10.1)') W

  write(*,*) '[R,S] + [R,[R,S]]/3'
  W = cmt(R,S) + (cmt(R,cmt(R,S))/3.0)
  write(*,'(4f10.1)') W

  write(*,*) '[R,S] + [R,[R,S]]/3 - [[R,S],[T,U]]/10'
  W = cmt(R,S) + (cmt(R,cmt(R,S))/3.0) - (cmt(cmt(R,S),cmt(T,U))/10.0)
  write(*,'(4f10.1)') W


contains

  function cmt(A, B) result(C)
    real(wp), intent(in) :: A(n), B(n)
    real(wp) :: D(n,1), E(1,n)
    real(wp) :: C(n,n)

    D = reshape(A,(/n,1/))
    E = reshape(B,(/1,n/))

    C = matmul(D,E)

  end function cmt
end program main

results

1 Like
4

Thanks for these replies. I have been debugging on-and-off for the past couple of days. I mentioned that I had been seeing some strange results – the code did not give the correct result for a certain test problem. After doing some editing I began to get segmentation faults. And it was the frustrating kind of segfault where you can put diagnostics at the beginning of subroutines and the code seems to crash upon entry (without even print hello from subroutine blah) for no good reason. There is some legacy code mixed in with this so I spent a lot of time seeing if something in the legacy code might be trashing memory.

In the end it appears that all my problems go away when I take a very long statement – one with 66 continuation lines – and break it into several smaller statements. I am really surprised by this. Has anyone else here experienced this type of situation?

5

As far I as know, there was 39 continuation lines allowed before Fortran 2003. And since Fortran 2003 it is 255.

6

The problem with too many continuation lines causing a segfault was happening on my MacBook (MacOS Sierra) with an old version of gfortran:
gfortran --version
GNU Fortran (Homebrew GCC 7.2.0) 7.2.0

So I moved the code to another platform and tried it with 3 different compilers: Cray, Intel, and gfortran 8.3.0 . There was no problem in those 3 cases.

So the problem appears to be due to using an old version of gfortran on my Mac.

1 Like
7

And if you use gfortran -std=f95 it should also cause an error.

8

There are hundreds of compiler errors when I use -std=f95 because of the use of :: in declarations, etc. So I’m not able to try it with that compiler option.

9

I wonder what declarations with :: are disallowed by -std=f95. Certainly Fortran 95 and 90 have the :: syntax, but is there a special case that only F2003+ allows?

1 Like
10

Here are some examples of the compiler diagnostics:

   real*8 function mzran()
        1

Error: GNU Extension: Nonstandard type declaration REAL*8 at (1)
and

   integer :: i,k,ip,kp
                      1

Error: Unexpected data declaration statement at (1)

and

   integer :: n
              1

Error: Unexpected data declaration statement at (1)

and

   real*8, dimension(monoms) :: f3,f4,f5,f6,f7,f8,tmp,tmp2,tmp3,tmp4,tmp5,tmp6
        1

Error: GNU Extension: Nonstandard type declaration REAL*8 at (1)

and

   g5(210:461)=g(210:461)
  1

Error: Unclassifiable statement at (1)

The list of compiler complaints goes on and on. Since this is a bit off topic I decided to move on. The main thing, from my standpoint, is that now if I compile the code with Intel, Cray, or Gnu compilers (Gnu being a recent version and not the old version I was using on my Mac), it compiles fine and executes correctly.

11

Yes, *8 is non standard, although a classical extension. REAL(8) or REAL(kind=8) is standard, but read that post for still better solutions:

12

Thank you for pointing me to the “Best way to declare a double…” thread. I read it and the “It takes all kinds” article mentioned within it. As I move forward I would like to write code that adheres to the Fortran standard, and those discussions were very helpful.

1 Like
13

Personally, I have improved my practice by using:

use, intrinsic :: iso_fortran_env, only: wp=>real64
implicit none
real(wp) :: alpha
1 Like