A Comparison of Programming Languages in Economics

I agree with this point, and I believe this issue can be solved with better documentation. This paper and its author are highly regarded in computational economics. Even he misuses double precision real literals, no wonder so many others make this mistake.

1 Like

At least it’s not C++…

2 Likes

Wish it was a little more specific about the compiler options used and so on. Since the code
was so short it was easy enough to get some numbers. On a very modest Linux laptop old enough to be the same vintage as the report, using a “before” and “after” version of the code making no other changes running each 100 times I got

100 runs
       Average (secs)    | Shortest   | Longest
Before 1.102803080400001 | 1.08397698 | 1.24383700 
After  1.1099996502      | 1.09448707 | 1.21478796 

with gfortran, which is quite a bit different than the 1.6 seconds in the article, but it was a very
different platform (a MSWindows box).

so any performance difference is small. I would have to change the code or instrument it to determine it any better, but there is some chance it could be up to a 1% but I doubt it.

It does produce slightly different answers, as can be expected from the above discussion

I just used fpm(1) and copied the file over app/main.f90 and fired it up using the default debug release and then --profile release. The results shown are for --profile release. The interesting part is the two profiles ran significantly different speeds (debug mode was 1.9 seconds, where the release mode was 1.1 seconds). So not sure what compiler options were used, and the report saying ifort speed changed by something like 7 or 9 percent could be because multiple runs were not averaged and so on.

Note the test is so small, that simple changes like running to the screen versus writing to an unbuffered file give around a five percent change on my platform, so I do not think my results change the general conclusions of the paper significantly, it just points out there might be some questions about the small changes compared to the previous results and such.

The only code change was to define a “dp” type and change all the constants to “dp” and the real(8) to real(dp).

PS: There is an odd transpose of an array and then a repeated transpose back that looked odd, but I did not look into it but it caught my eye as I did not notice the array as being changed anywhere and the way it was repeated would slow down an interpreted language.

2 Likes

It’s even more convenient to write
integer, parameter :: dp = kind(1d0)

If the link to Fortran code provided upthread is truly reflective of their performance data in the first paper, the results mentioned in the paper for Fortran paper are unreliable. It’s rather poor instrumentation of the code, it’s absolutely at a beginner level. Immediately the lack of the following is glaring:

  • in situ the code can be instrumented to run a certain number of trials, say NTRIALS where it is a reasonable number such as a dozen. Discarding then the fastest and the slowest times, the CPU usage then be marked as the median of the rest. Then repeat the same with certain number of execution attempts, say NRUNS, and get the overall median,
  • use SYSTEM_CLOCK to get the tick counts using an integer type with a range of at least 10 (colloquially a 64-bit integer or better) rather than CPU_TIME to measure timing,
  • separate in code with suitable instrumentation the measurement of strictly numerical computing instructions from those involving output (and input) involving PRINT (and READ) statements. And report them distinctly. If the times toward IO are on the same order of magnitude as numerical operations, there is little to no reason for those keen on Fortran to pay much attention to such studies, the IO time will not be that competitive; besides, that will be of little practical interest any way.

It would be worth experimenting with a BLAS call to ?gemm for the operation:

     expectedValueFunction = matmul(mValueFunction,transpose(mTransition))

I also agree with @FortranFan’s observation to separate the I/O part. Currently, the printing is performed quite frequently, every 10 iterations.

The textbook Introduction to Computational Economics using Fortran also uses real*8 to declare double precision real numbers.
https://www.ce-fortran.com/source-codes/

Again, the changes do not immediately change the general conclusions, but
it seems like a good place to show how

  • the warning flags of the compiler should be used to their fullest
  • a few compiler options can affect speed significantly,
  • the often overlooked gprof(1) and gcov(1) tools can tell you a
    lot about how your code is used and where it spends time with
    very little effort.

So for users with the GCC programs gfortran, gprof, and gcov you should
give something like this a try if you are interested in information for
improving performance or testing coverage …

#!/bin/bash
(
exec 2>&1
# do everything in a scratch directory
SCR=scr_$(uuidgen)
mkdir -p $SCR
cp econ1.f90 $SCR/
cd $SCR
gfortran -g -pg -fprofile-arcs -ftest-coverage -O0 econ1.f90 -o econ1
./econ1
# better when you have a lot of procedures
gprof --brief econ1 gmon.out  
# better for short programs
gprof --brief --line econ1 gmon.out 
# coverage report
gcov --stdout econ1.f90
) >reports.out
exit

and simply compiling with -Warn -Wextra identified a few easily corrected
issues. Always a good practice to use all those optional warnings!

And I was off a bit on my guess about performance and using 64-bit
constants instead of 32-bit. Using a more rigorious timing test and 1 000,
executions it does slow it down on my platform 1.36 percent. I am still
a bit surprised, as almost all the usage of 32-bit values was combined
with 64-bit values.

By the way, for fun I built it using 128-bit values and got about the
same answer but it was 45 times slower (with gfortran). I expected a
slowdown but not of that order either. Live and learn.

So just using those example commands you can tell what lines of your
code were executed via gcov(1) (in this trivial case not so interesting
except for the counts but if building a test suite that is good info.).

And you can tell that three lines account for over 1/2 of the CPU time
consumed.

In this case the procedure report is essentially null, but in a big
program it quickly can identify which procedures are consuming the
CPU time.

gprof(1) and gcov(1) have their limitations but are often very useful. You
should read up on them and similar tools if this is interesting.

And as a footnote using unrolling and a few other switches on my platform
took ifort(1) to the top of the heap at 1`04``.

reports program output
 Steady State values
 Output:   0.56273143387317726      Capital:   0.17819828739139082      Consumption:   0.38453314648178644     
 Iteration:           1 Sup Diff:   5.2741593406882385E-002
 Iteration:          10 Sup Diff:   3.1346949265582680E-002
 Iteration:          20 Sup Diff:   1.8703459893196328E-002
 Iteration:          30 Sup Diff:   1.1165512033874614E-002
 Iteration:          40 Sup Diff:   6.6685417080749598E-003
 Iteration:          50 Sup Diff:   3.9842927486826163E-003
 Iteration:          60 Sup Diff:   2.3813118039122116E-003
 Iteration:          70 Sup Diff:   1.4236586450859789E-003
 Iteration:          80 Sup Diff:   8.5133977471318900E-004
 Iteration:          90 Sup Diff:   5.0920517522479170E-004
 Iteration:         100 Sup Diff:   3.0462324421209885E-004
 Iteration:         110 Sup Diff:   1.8226485632144573E-004
 Iteration:         120 Sup Diff:   1.0906950872535681E-004
 Iteration:         130 Sup Diff:   6.5276437362649098E-005
 Iteration:         140 Sup Diff:   3.9071082119535028E-005
 Iteration:         150 Sup Diff:   2.3388077119768091E-005
 Iteration:         160 Sup Diff:   1.4008644636964718E-005
 Iteration:         170 Sup Diff:   8.3913172028715621E-006
 Iteration:         180 Sup Diff:   5.0264745379280384E-006
 Iteration:         190 Sup Diff:   3.0108998635425266E-006
 Iteration:         200 Sup Diff:   1.8035522479920019E-006
 Iteration:         210 Sup Diff:   1.0803409159487742E-006
 Iteration:         220 Sup Diff:   6.4713169423136208E-007
 Iteration:         230 Sup Diff:   3.8763619381043668E-007
 Iteration:         240 Sup Diff:   2.3219657918627234E-007
 Iteration:         250 Sup Diff:   1.3908720930544405E-007
 Iteration:         257 Sup Diff:   9.7160356538061876E-008
  
 My check:  0.14654914369569541     
  
 Elapsed time is    1.73320293    

gprof by procedure

Flat profile:

Each sample counts as 0.01 seconds.
  %   cumulative   self              self     total           
 time   seconds   seconds    calls  ms/call  ms/call  name    
100.10      0.75     0.75        1   750.73   750.73  MAIN__

                        Call graph


granularity: each sample hit covers 2 byte(s) for 1.33% of 0.75 seconds

index % time    self  children    called     name
                0.75    0.00       1/1           main [2]
[1]    100.0    0.75    0.00       1         MAIN__ [1]
-----------------------------------------------
                                                 <spontaneous>
[2]    100.0    0.00    0.75                 main [2]
                0.75    0.00       1/1           MAIN__ [1]
-----------------------------------------------
 
Index by function name

   [1] MAIN__ (econ1.f90)
Flat profile:

Each sample counts as 0.01 seconds.
  %   cumulative   self              self     total           
 time   seconds   seconds    calls  Ts/call  Ts/call  name    
 28.69      0.22     0.22                             MAIN__ (econ1.f90:75 @ 2d76)
 21.35      0.38     0.16                             MAIN__ (econ1.f90:76 @ 2de1)
 16.02      0.50     0.12                             MAIN__ (econ1.f90:74 @ 2d22)
 10.68      0.58     0.08                             MAIN__ (econ1.f90:88 @ 2eb8)
  4.67      0.61     0.04                             MAIN__ (econ1.f90:73 @ 2cf5)
  4.00      0.64     0.03                             MAIN__ (econ1.f90:85 @ 2e71)
  2.67      0.66     0.02                             MAIN__ (econ1.f90:71 @ 2cba)
  2.67      0.68     0.02                             MAIN__ (econ1.f90:73 @ 2e21)
  2.67      0.70     0.02                             MAIN__ (econ1.f90:71 @ 2ea4)
  2.00      0.72     0.02                             MAIN__ (econ1.f90:77 @ 2def)
  2.00      0.73     0.02                             MAIN__ (econ1.f90:84 @ 2e3e)
  1.33      0.74     0.01                             MAIN__ (econ1.f90:78 @ 2df9)
  0.67      0.75     0.01                             MAIN__ (econ1.f90:81 @ 2ded)
  0.67      0.75     0.01                             MAIN__ (econ1.f90:76 @ 2e2a)
  0.00      0.75     0.00        1     0.00     0.00  MAIN__ (econ1.f90:6 @ 2429)

                        Call graph


granularity: each sample hit covers 2 byte(s) for 1.33% of 0.75 seconds

index % time    self  children    called     name
                0.00    0.00       1/1           main (econ1.f90:104 @ 3680) [70]
[15]     0.0    0.00    0.00       1         MAIN__ (econ1.f90:6 @ 2429) [15]
-----------------------------------------------

Index by function name

  [15] MAIN__ (econ1.f90:6 @ 2429) [2] MAIN__ (econ1.f90:76 @ 2de1) [14] MAIN__ (econ1.f90:76 @ 2e2a)
   [7] MAIN__ (econ1.f90:71 @ 2cba) [13] MAIN__ (econ1.f90:81 @ 2ded) [11] MAIN__ (econ1.f90:84 @ 2e3e)
   [5] MAIN__ (econ1.f90:73 @ 2cf5) [10] MAIN__ (econ1.f90:77 @ 2def) [6] MAIN__ (econ1.f90:85 @ 2e71)
   [3] MAIN__ (econ1.f90:74 @ 2d22) [12] MAIN__ (econ1.f90:78 @ 2df9) [9] MAIN__ (econ1.f90:71 @ 2ea4)
   [1] MAIN__ (econ1.f90:75 @ 2d76) [8] MAIN__ (econ1.f90:73 @ 2e21) [4] MAIN__ (econ1.f90:88 @ 2eb8)

coverage report. Note the lines are numbered and match up to the gprof lines

econ1.gcno:version 'B03*', prefer 'A93*'
econ1.gcda:version 'B03*', prefer version 'A93*'
        -:    0:Source:econ1.f90
        -:    0:Graph:econ1.gcno
        -:    0:Data:econ1.gcda
        -:    0:Runs:1
        -:    1:!============================================================================
        -:    2:! Name        : RBC_F90.f90
        -:    3:! Description : Basic RBC model with full depreciation
        -:    4:! Date        : July 21, 2013
        -:    5:!============================================================================
        1:    6:program RBC_F90
        -:    7:!----------------------------------------------------------------
        -:    8:! 0. variables to be defined
        -:    9:!----------------------------------------------------------------
        -:   10:implicit none
        -:   11:integer,  parameter :: dp=kind(0.0d0)
        -:   12:integer,  parameter :: nGridCapital = 17820
        -:   13:integer,  parameter :: nGridProductivity = 5
        -:   14:real(dp), parameter :: tolerance = 0.0000001_dp
        -:   15:real     :: time_begin, time_end
        -:   16:integer  :: nCapital, nCapitalNextPeriod, gridCapitalNextPeriod, nProductivity
        -:   17:integer  :: iteration
        -:   18:real(dp) :: aalpha, bbeta, capitalSteadyState, outputSteadyState, consumptionSteadyState
        -:   19:real(dp) :: valueHighSoFar, valueProvisional, consumption, capitalChoice
        -:   20:real(dp) :: maxDifference
        -:   21:real(dp), dimension(nGridProductivity)                   :: vProductivity
        -:   22:real(dp), dimension(nGridProductivity,nGridProductivity) :: mTransition
        -:   23:real(dp), dimension(nGridCapital)                        :: vGridCapital
        -:   24:real(dp), dimension(nGridCapital,nGridProductivity)      :: mOutput, mValueFunction, mValueFunctionNew, mPolicyFunction
        -:   25:real(dp), dimension(nGridCapital,nGridProductivity)      :: expectedValueFunction
        1:   26:  call CPU_TIME(time_begin)
        -:   27:  !----------------------------------------------------------------
        -:   28:  ! 1. Calibration
        -:   29:  !----------------------------------------------------------------
        1:   30:  aalpha = 0.33333333333_dp; ! Elasticity of output w.r.t
        1:   31:  bbeta  = 0.95_dp ! Discount factor
        -:   32:  ! Productivity value
        1:   33:  vProductivity = (/0.9792_dp, 0.9896_dp, 1.0000_dp, 1.0106_dp, 1.0212_dp/)
        -:   34:  ! Transition matrix
        -:   35:  mTransition = reshape( (/ &
        -:   36:   0.9727_dp,  0.0273_dp,  0._dp,      0._dp,      0._dp,      &
        -:   37:   0.0041_dp,  0.9806_dp,  0.0153_dp,  0.0_dp,     0.0_dp,     &
        -:   38:   0.0_dp,     0.0082_dp,  0.9837_dp,  0.0082_dp,  0.0_dp,     &
        -:   39:   0.0_dp,     0.0_dp,     0.0153_dp,  0.9806_dp,  0.0041_dp,  &
        1:   40:   0.0_dp,     0.0_dp,     0.0_dp,     0.0273_dp,  0.9727_dp   /),  (/5,5/))
        -:   41:  !----------------------------------------------------------------
        -:   42:  ! 2. Steady state
        -:   43:  !----------------------------------------------------------------
        1:   44:  capitalSteadyState     = (aalpha*bbeta)**(1.0_dp/(1.0_dp-aalpha))
        1:   45:  outputSteadyState      = capitalSteadyState**(aalpha)
        1:   46:  consumptionSteadyState = outputSteadyState-capitalSteadyState
        1:   47:  print *, 'Steady State values'
        1:   48:  print *, 'Output: ', outputSteadyState, 'Capital: ', capitalSteadyState, 'Consumption: ', consumptionSteadyState
        -:   49:  ! Grid for capital
    17821:   50:  do nCapital = 1, nGridCapital
    17821:   51:     vGridCapital(nCapital) = 0.5_dp*capitalSteadyState+0.00001_dp*(nCapital-1)
        -:   52:  enddo
        -:   53:  !----------------------------------------------------------------
        -:   54:  ! 3. Pre-build Output for each point in the grid
        -:   55:  !----------------------------------------------------------------
        6:   56:  do nProductivity = 1, nGridProductivity
    89106:   57:     do nCapital = 1, nGridCapital
    89105:   58:        mOutput(nCapital, nProductivity) = vProductivity(nProductivity)*(vGridCapital(nCapital)**aalpha)
        -:   59:     enddo
        -:   60:  enddo
        -:   61:  !----------------------------------------------------------------
        -:   62:  ! 4. Main Iteration
        -:   63:  !----------------------------------------------------------------
        1:   64:  maxDifference = 10.0_dp
        1:   65:  iteration     = 0
      258:   66:  do while (maxDifference>tolerance)
      257:   67:     expectedValueFunction = matmul(mValueFunction,mTransition);
     1542:   68:     do nProductivity = 1,nGridProductivity
        -:   69:        ! We start from previous choice (monotonicity of policy function)
     1285:   70:        gridCapitalNextPeriod = 1
 22900242:   71:        do nCapital = 1,nGridCapital
 22898700:   72:           valueHighSoFar = -100000.0_dp
 60486291:   73:           do nCapitalNextPeriod = gridCapitalNextPeriod,nGridCapital
 60486291:   74:              consumption = mOutput(nCapital,nProductivity)-vGridCapital(nCapitalNextPeriod) !88
 60486291:   75:              valueProvisional = (1.0_dp-bbeta)*log(consumption)+bbeta*expectedValueFunction(nCapitalNextPeriod,nProductivity) !89
60486291*:   76:              if (valueProvisional>valueHighSoFar) then ! we break when we have achieved the max !90
 37587591:   77:                 valueHighSoFar        = valueProvisional
 37587591:   78:                 capitalChoice         = vGridCapital(nCapitalNextPeriod)
 37587591:   79:                 gridCapitalNextPeriod = nCapitalNextPeriod
        -:   80:              else
 22898700:   81:                 exit
        -:   82:              endif
        -:   83:           enddo
 22898700:   84:           mValueFunctionNew(nCapital,nProductivity) = valueHighSoFar
 22899985:   85:           mPolicyFunction(nCapital,nProductivity)   = capitalChoice
        -:   86:        enddo
        -:   87:     enddo
22900242*:   88:     maxDifference = maxval((abs(mValueFunctionNew-mValueFunction)))
      257:   89:     mValueFunction = mValueFunctionNew
      257:   90:     iteration = iteration+1
      257:   91:     if (mod(iteration,10)==0 .OR. iteration==1) then
       26:   92:        print *, 'Iteration:', iteration, 'Sup Diff:', MaxDifference
        -:   93:     endif
        -:   94:  enddo
        -:   95:  !----------------------------------------------------------------
        -:   96:  ! 5. PRINT RESULTS
        -:   97:  !----------------------------------------------------------------
        1:   98:  print *, 'Iteration:', iteration, 'Sup Diff:', MaxDifference
        1:   99:  print *, ' '
        1:  100:  print *, 'My check:', mPolicyFunction(1000,3)
        1:  101:  print *, ' '
        1:  102:  call CPU_TIME (time_end)
        1:  103:  print *, 'Elapsed time is ', time_end - time_begin
        1:  104:end program RBC_F90

I do find it a bit depressing that the, apparently respected, author of a textbook for no obvious reason except perhaps ignorance doesn’t use Standard Fortran.

If the author had called it “Introduction to Computational Fortran using a language similar to Fortran but non-standard” that might be been acceptable, but otherwise it is not. Don’t textbooks go through any form of review process any more?

There is something like 12 languages in the paper. Probably none of the 2 or 3 reviewers were familiar with Fortran…

Yes, trying to cover multiple languages with a real-world significant case that covers circumstances where specialists may not even be available is handled relatively well by the paper and it’s conclusions are reasonable for the intended audience, but I do wish examples like this were reviewed. Perhaps it is up to the Fortran Discourse community to emphasize that it welcomes reviewing such examples before they are committed to posterity. Not all forums would. But, as noted above this was not in a Fortran textbook and very specifically was not even concentrating on Fortran. The bright side is Fortran was one of the languages considered required to be included.
I am sure if a major review of the other examples were done they could also be improved. I don’t think the audience was looking at this as a resource for finding good Fortran material; I find this much more disturbing when I find similar issues in major Fortran projects and references though.

An unintended consequence of the non-optimal use is that a transpose of a transpose was being done (essentially, a no-op) and one of them was in a loop and all the compilers I tried managed to identify and eliminate that so that correcting it had no significant change in the results. Quite a few interpreted languages with less optimization would likely not have done that; so an unstated conclusion for the paper might be that Fortran is still fast even for someone not an expert in the language; which some of the other results indicate is not true for other languages. It would be an interesting follow up if this were done with domain experts optimizing each example. I suspect that the really slow ones would improve dramatically but would require in-depth knowledge of the language; but that Fortran was near-optimal written in a intuitive form. Just my hunch, would be interesting to see that put to the test. That would be a plus for Fortran for a good sized population of people using small off-the-cuff or single-user programs.

This reminds me of something. To calculate the covariance or correlation matrix of the columns of matrix A, you need to compute

matmul(transpose(A),A)

which results in a symmetric matrix. I wonder if there is a code that speeds up the calculation by recognizing this symmetry. I have tried to do so by writing loops that calculate only the upper diagonal elements of the result, but that function was slower than
the matrix operations that ignore the symmetry.

One answer to this has been posted here:

The Julia code generated by the Linnea tool is

using LinearAlgebra.BLAS
using LinearAlgebra

"""
    algorithm0(ml0::Array{Float64,2})

Compute
A = (X^T X).

Requires at least Julia v1.0.

# Arguments
- `ml0::Array{Float64,2}`: Matrix X of size 100 x 100 with property Non-singular.
"""                    
function algorithm0(ml0::Array{Float64,2})
    # cost: 1e+06 FLOPs
    # X: ml0, full
    ml1 = Array{Float64}(undef, 100, 100)
    # tmp1 = (X^T X)
    syrk!('L', 'T', 1.0, ml0, 0.0, ml1)

    # tmp1: ml1, symmetric_lower_triangular
    for i = 1:100-1;
        view(ml1, i, i+1:100)[:] = view(ml1, i+1:100, i);
    end;

    # tmp1: ml1, full
    # A = tmp1
    return (ml1)
end

It shouldn’t be too hard to adapt this to Fortran. The main differences are the call to ?syrk has a few more arguments.

subroutine tranposeAtimesA(k,n,A,C)
   integer, intent(in) :: k, n
   real(dp), intent(in) :: A(k,n)
   real(dp), intent(inout) :: C(n,n)

   real(dp), parameter :: alpha = 1.0d0, beta = 0.0d0
   integer :: i

  call dsyrk('L','T',n,k,alpha,A,k,beta,C,n)
 
  ! the lower triangular part of C should now contain matmul(transpose(A),A)
  ! we need to copy it into the upper triangular part
   do i = 1, n
      C(i,i+1:n) = C(i+1:n,i)
   end do
end subroutine

:warning: Disclaimer: I haven’t attempted to compile this or verify if it works correctly. I also don’t know if the view calls in Julia are somehow accessing the arrays in a more favorable order. The Fortran version attempts to copy the lower-triangular columns (excluding the diagonal) into the upper-triangular rows.

2 Likes

Even after fixing the accuracy issues introduced in the original code by the authors, I observed other issues in the measured times. Running the C++ and the Fortran codes on my machine, I found the Intel compiler coming on top at 2X the speed of C++. These are my results (allowing highest optimization options for all compilers):

! My Check = 0.1465491436962635
 ! Elapsed time is = 1.24415     C++
 ! Elapsed time is   1.26562500  gfortran 
 ! Elapsed time is   0.6093750   ifort        
 ! Elapsed time is   0.6677110   nvfortran
 ! Elapsed time is   0.861564    Julia 1.8.0
1 Like

I find problems more with 64-bit integer constants, where 1024 is typically used when 1024_dp can be required.
eg

integer*8 ::i ; i = 1024*1024*1024*4  can lead to a 32-bit calculation with some compilers, while
"i = 12345678901" will be flaged as a problem with most compilers.
1 Like

Julia community has this practice of reviewing new manuscripts of books. I am actually curious about how many textbooks using Fortran follow the modern Fortran standard.
The FORTRAN FOR SCIENTISTS & ENGINEERS book by Stephen Chapman only uses upper cases for Fortran codes.

The 4th edition of Chapman’s book uses upper case for keywords but not variable names or intrinsic functions. Here is one of the first programs in the book.

PROGRAM roots
! Purpose:
! This program solves for the roots of a quadratic equation of the form
! A * X**2 + B * X + C = 0. It calculates the answers regardless of the
! type of roots that the equation possesses (Fortran 95/2003 style).
!
IMPLICIT NONE
! Declare the variables used in this program
REAL :: a ! Coefficient of X**2 term of equation
REAL :: b ! Coefficient of X term of equation
REAL :: c ! Constant term of equation
REAL :: discriminant ! Discriminant of the equation
REAL :: imag_part ! Imaginary part of equation (for complex roots)
REAL :: real_part ! Real part of equation (for complex roots)
REAL :: x1 ! First solution of equation (for real roots)
REAL :: x2 ! Second solution of equation (for real roots)
! Prompt the user for the coefficients of the equation
WRITE (*,*) 'This program solves for the roots of a quadratic '
WRITE (*,*) 'equation of the form A * X**2 + B * X + C = 0. '
WRITE (*,*) 'Enter the coefficients A, B, and C:'
READ (*,*) a, b, c
! Echo back coefficients
WRITE (*,*) 'The coefficients A, B, and C are: ', a, b, c
! Calculate discriminant
discriminant = b**2 - 4. * a * c
! Solve for the roots, depending upon the value of the discriminant
IF ( discriminant > 0. ) THEN ! there are two real roots, so...
X1 = ( -b + sqrt(discriminant) ) / ( 2. * a )
X2 = ( -b - sqrt(discriminant) ) / ( 2. * a )
WRITE (*,*) 'This equation has two real roots:'
WRITE (*,*) 'X1 = ', x1
WRITE (*,*) 'X2 = ', x2
ELSE IF ( discriminant == 0. ) THEN ! there is one repeated root, so...
x1 = ( -b ) / ( 2. * a )
WRITE (*,*) 'This equation has two identical real roots:'
WRITE (*,*) 'X1 = X2 = ', x1
ELSE ! there are complex roots, so ...
real_part = ( -b ) / ( 2. * a )
imag_part = sqrt ( abs ( discriminant ) ) / ( 2. * a )
WRITE (*,*) 'This equation has complex roots:'
WRITE (*,*) 'X1 = ', real_part, ' +i ', imag_part
WRITE (*,*) 'X2 = ', real_part, ' -i ', imag_part
END IF
END PROGRAM roots

Recent Fortran textbooks usually do present codes that follow recent standards. It’s the books that use Fortran but are not primarily about Fortran that often do not, for example here (the book is not published yet – maybe the codes will be revised).

Do you suggest that being wrong? The F2018 Standard itself uses upper case in all code snippets shown in NOTEs (both for keywords and (most of) variables) and in all references to the keywords throughout the text. Personally I find the mixed mode (uppercase keywords, lowercase objects) quite useful, especially if the code is intended to be shown as an example (say, to students)

I am not suggesting that is wrong. It is just not easy to read.
A textbook can influence thousands of potential users, which requires the coding style to be as good as it can be. Based on my personal experience, those who code Fortran keywords in upper cases also tend to code the other parts in upper cases. I totally agree that mixed cases are useful.

I tend so use uppercase for global variables in a module or for parameter values only, and not always even for those, but the trial offer of fortPLUS/spag had a large amount of options for case, allowing for parameters of procedures to have a leading uppercase character, for fixed values and/or keywords to be uppercase, and so on. It was interesting to run code through it and look at the various results. Without a tool like that I find it rather tedious to follow a lot of rules not enforced by any compilers (that I know of) so I rarely use case as being significant.

Many people preferred uppercase even when the compiler had an extension to not require it in the past, because uppercase was far more legible on low-res terminals, so the comment about it being unpleasant to read evoked “how times have changed” .

Several people have told me they are used to uppercase meaning “shout” and so find it very disturbing to see all of a code in uppercase.

WIth a lot of utilities that do colored syntax highlighting I find case even less useful than I used to find it (for “highlighting” certain words);to the point I made a little filter program to change everything not in a comment to lowercase that I use quite frequently (only works on free-format Fortran, and does not do anything to cpp(1) lines, which can cause problems with macros if you use them a lot).

Are there other tools other than plusFORT/spag that apply case rules consistently to code?