What do we mean by common?

1

I have two include files that define variables by using the COMMON declaration.

An example of this is…

  COMMON / ZONES / SFILE,OFILE,PFILE,LFILE,FILEZ,radM,radR,Nrad,       &
     &              unitz,start_zones,increment_zones,iunit,lunit

I believe that this groups all of the variables on the right under a common label of ‘zones’ (in the example). Yet, the code I have that includes these two files references the variables on the right directly. So, what’s the point?

As far as I can see, this is a way of having a centrally-defined group of variables that can then be used throughout a program (as long as you include the files into the subroutines).

Yet, the data type isn’t defined, nor are they initialised with any values.

The other file contains lines such as…

PARAMETER (MJ=4000,MH=4,ME=9,MTAU=20000)

Here I can see an initial value, but no data type and what is meant by PARAMETER?

2

For new code you should avoid common at all costs. It is an obsolescent feature (officially declared as such in the standard). However, assuming you’re trying to understand some old code, common is a way of having “global” variables. All variables declared in a common block are accessible to any program unit that declares the same common block. However, there are lots of pitfalls here.

The compiler does not (in many cases can not) check that the different declarations of the common block are the same. This has tricky implications. As you’ve noticed, the types of the variables are not declared within the common block. If the types between declarations are different, you’ll get weird behavior (i.e. interpreting the bits of a real value as though it were an integer). It does not check that arrays have the same shapes, and thus may inadvertently overlap/shift the values of other variables. It does not check that the variables are in the same order, and thus potentially accessing each other’s data instead. The common practice was to declare a common block and its variables in a single file, and include that file wherever it was used, to prevent chances of mismatches.

The parameter statement says that these variables are constants with these values.

If you’re variables do not have type declarations, that probably means they are using implicit typing, and thus variables that start with i-n are integer, and all others are real. Thus comes the old joke, God is real, because it’s undeclared, and undeclared variables starting with G are real.

Working with these old codes can at times be confusing and frustrating, so good luck and don’t hesitate to come back here asking questions. Also, @rouson and I work on “legacy” codes like these professionally, so don’t hesitate to reach out if you need it.

2 Likes
3

Thank you SO much for this.
Yes, I am taking a large program (well, actually not that bad at around 8000 lines of code) and first getting my head around it, then tidying it up before redesigning it to make it more efficient (if possible).
I am thinking of shifting this to free-form (F90?) code, since I really would like to move away from the line length issue AND stop having to type things like .EQ.
One question is, though; is the efficiency of the free-form compiler as good (or better) than the older approach?

Your information on the use of the two files to help centralise the COMMON and PARAMETER fits with what I thought it was. My question here is; if I am to remove the COMMON (and, yes, I hate the fact that this is using implicit definitions - way too dangerous), is there an alternative I should be implementing? Or is this a case of defining these in the normal way (INTEGER :: i = 0), executing a save (if I need the data to be static) and parsing everything through subroutine parameters?

Sorry to be extending this further (you did offer!) but having a single point where I define a value (let’s say the gravitational constant) makes good sense. To have to redefine it in every relevant subroutine introduces new risks of mistakes but there is no global structure within Fortran.

Should I define a file that simply defines all my ‘global’ variables in the standard way, and just include this in the relevant subroutines?

4

Working with these old codes can at times be confusing and frustrating

Yes, but isn’t that where the fun lies? :laughing:

1 Like
5

Your looking for modules. Define your constants in a module(s), and use that module, possibly (IMHO preferably) with the only clause. I.e.

module my_constants
  implicit none
  real, parameter :: pi = 3.14
  integer, parameter :: the_answer = 42
end module

program main
  use my_constants, only: the_answer
  implicit none

  print *, "The answer to life, the universe and everything is ", the_answer
end program

To the extent possible you should avoid “globally” accessible, mutable data (i.e. module variables that aren’t parameters) and instead pass things as arguments. Of course, starting from an existing code that doesn’t do that, you should take baby steps to move in that direction.

6

Yep. The “how does this work?” and “what were they thinking?” puzzles are sometimes quite intriguing.

7

You can use == in fixed form, and the efficiency of a program does not depend on whether the source code is free or fixed form. The compiler will convert the code to a representation that is independent of such details.

2 Likes
8

My concern is that any change to the syntax of a language would require a rewrite of the compiler (at least in part) and this could easily lead to inefficient translation.

One of the strengths of Fortran, it appears, is the efficiency of the compiler. I take it from your reply that this is not a valid concern. :slight_smile:

Thank you.

9

And this is the EXACT reason why I like to use forums. I would not have found this solution easily without your help!

Yes, just a single change can cause the code to go awry (as it is currently doing).

10

Short answer: yes, although technically as Beliavsky has explained, it’s not the source form (fixed- vs free-) but rather the type of language constructs used historically vs today.

There’s a nice benchmark example from Michael Hirsch, comparing a Fortran II vs modern Fortran inverse tangent function. For the example subprogram, the difference can be a whole order of magnitude.

The example is taken from the IBM Fortran II manual (pg. 29). The subprogram evaluates arctan(x) accurate to four decimal places, for any positive input x. The expansions used are:

\mathrm{arctan}(x) = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + \dots

for the region 0 \leq x \leq 1, and

\mathrm{arctan}(x) = \frac{\pi}{2} - \frac{1}{x} + \frac{\left(\frac{1}{x}\right)^3}{3} - \frac{\left(\frac{1}{x}\right)^5}{5} + \frac{\left(\frac{1}{x}\right)^7}{7} - \dots

where x > 1.

:warning: Note: the example below only serves to demonstrate different styles of programming. To evaluate the inverse tangent in practice, use the intrinsic function atan(x).

Here is a slightly modified example of the Fortran II code:

      function arctan(x)
      if (x) 2,3,3
2     arctan = -99   ! instead of stop,
      return         ! return a value outside plausible range
3     arctan = 0.0
      if (x - 1.0) 10, 10, 5
5     term = -1.0/x
      arctan = 1.57079    ! pi / 2
      goto 11
10    term = x
11    prevxp = 1.0
      y = term ** 2.0
12    arctan = arctan + term
      presxp = prevxp + 2.0
13    term = -prevxp / presxp * y * term
      prevxp = presxp
14    if(term - 0.00005) 15, 12, 12
15    if(-term - 0.00005) 16, 12, 12
16    return
      end

Here’s what the new code looks like:

elemental real function arctan(x)
  use, intrinsic :: ieee_arithmetic, only: &
    ieee_value, ieee_quiet_nan

  real, intent(in) :: x
  real :: prevxp, presxp, term, y
  real :: nan
  real, parameter :: pi = 4.*atan(1.)
  integer :: i

  nan = ieee_value(0., ieee_quiet_nan)

  if (x < 0.) then
    arctan = nan
    return
  endif

  arctan = 0.

  if (x - 1. > 0.) then
    term = -1. / x
    arctan = pi/2
  else
    term = x
  endif

  prevxp = 1.
  y = term ** 2.

  do while(term-0.00005 >= 0. .or. -term - 0.00005 >= 0.)
    arctan = arctan + term
    presxp = prevxp + 2.

    term = -prevxp / presxp * y * term
    prevxp = presxp
  enddo

end function arctan
1 Like
11

@ivanpribec any idea why the modern version would be faster than the older version? The code seems equivalent at first look.

12

Strictly speaking, a PARAMETER attribute on a name makes it a “named constant”, not a variable. It does not have an existence at runtime itself, but it can be used to form expressions.

If you are feeling mean, give a compiler a huge PARAMETER array whose only reference is in a SIZE() intrinsic and watch it struggle needlessly.

2 Likes
13

Both versions cheat, since they directly use the value of π for the case x > 1, rather than doing the calculation using the series. This reminds me of a former professor who kept a plaque in his office on which was engraved Old age and treachery beats youthful exuberance.

Both versions fail to work at all for negative arguments, and do so without an error message (unless we think of “-99” or “NaN” as an informative message). The old version has the last digit wrong (not correctly rounded) in the value 1.57079.

I recognize that you are comparing different styles of programming, but let the OP note that using the Leibniz series is one of the most inefficient and slow methods of calculating π, since one requires about five billion terms to obtain a result with ten correct decimal digits when x = 1. It is an instructive exercise to estimate how long it would take a human to perform the same calculation using pencil and paper. Or on the adorable IBM 1620 computer.

A subtle point to note is that sometimes what we do before writing the program, such as considering several algorithms before selecting one, can be more important than how we implement and run the selected algorithm.

4 Likes
14

The general approach to solve this is to use a taylor expansion of arctan(x-1) for inputs near 1. Julia uses a slightly more complicated approach, with the following brances

[0,7/16]      atan(x) = polynomial approximation
[7/16,11/16]  atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
[11/16,19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
[19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
[39/16,Inf]   atan(x) = atan(Inf) + atan( -1/t )

These branches provide <1ulp with a 11 degree polynomial for Float64 and a 5 degree polynomial for Float32.

2 Likes
15

Thanks for the warnings, I will edit my previous post with a note.

Edit: here’s a plot of the poor convergence of the Leibniz series near x = 1:

Click to show plot

image

For the record, the IBM Fortran II manual was aware this is not the best way to evaluate the arctangent, quoting:

Example 2: Series Evaluation

The following subprogram is an example of a series evaluation with IF-type branching. The function defined by this subprogram cannot be defined in a single statement. The subprogram computes the value of \mathrm{arctan}\,x, correct to 4 decimal places, for any given x greater than or equal to zero. Actually, the arctangent function is available on the FORTRAN II library tape and, in practice, would not normally be rewritten as a function subprogram. [emphasis added] It will, however, serve to illustrate a type of problem.

2 Likes
16

I must confess, I misinterpreted the benchmark timings. The intrinsic function atan(x) is roughly one order of magnitude faster. As Michael states on GitHub:

Using the same algorithm, modern Fortran 95+ techniques are approximately the same speed as awkward Fortran 77 style coding.

My answer to the OP remains, the “free-form” compiler is as good as the old approach.

1 Like
17

Conclusion: the new approach is as good as the old approach. Also the old approach is as good as the new approach. :slight_smile:

18

In that case, here is a compact version of the old version, in a “write-only” style that the compiler is happy to digest, but not human readers:

      FUNCTIONARCTAN(X);IF(X)2,3,3
    2 ARCTAN=-99;RETURN
    3 ARCTAN=0.0;IF(X-1.0)10,10,5
    5 TERM=-1.0/X;ARCTAN=1.57079;GOTO11
   10 TERM=X
   11 PREVXP=1.0;Y=TERM**2.0
   12 ARCTAN=ARCTAN+TERM;PRESXP=PREVXP+2.0
   13 TERM=-PREVXP/PRESXP*Y*TERM;PREVXP=PRESXP
   14 IF(TERM-0.00005)15,12,12
   15 IF(-TERM-0.00005)16,12,12
   16 RETURN;END

By the way, let us note that at one time humans wrote/scribed text with no punctuation, not even spaces between words, a style named Scriptio Continua.

1 Like
19

If runtime is the only measure. In many cases, readability is more important.

2 Likes
20

My intention is to slowly move to arguments; though the subroutines already have LOTS of arguments! Hahaha!

The code I have just opens up into the initial routine (no program main); should I change this?

Equally, before I do, can I just place the module at the top and reference it (using the use and only statements) in the existing code (given that there is no initial subroutine defined)?