Fortran code snippets

1

As an experiment I’d like to start a thread where people can present snippets of Fortran code. I’ll start.

There are continuous symmetric statistical distributions that span (-Inf,Inf), such as the normal distribution, Student t, and Laplace. One may want to simulate data with skew, for example to simulate daily stock market returns, which have some negative skew. In a paper On Bayesian Modeling of Fat Tails and Skewness, Fernandez and Steel suggested scaling positive and negative variates by different amounts to create skew. Here is a function that takes input data with mean and variance assumed to be 0 and 1 and creates skewed data with the same expected mean and variance but with nonzero skew.

module asymm_mod
implicit none
integer, parameter :: dp = kind(1.0d0)
private
public :: asymm
contains
pure function asymm(x,c,abs_x) result(y)
! create asymmetry for symmetric zero-mean, unit-variance data by differentially scaling positive and negative values of x(:)
real(kind=dp), intent(in) :: x(:)
real(kind=dp), intent(in) :: c     ! asymmetry parameter: > 1 to create positively skewed data, 0 < c < 1 for negatively skewed data
real(kind=dp), intent(in) :: abs_x ! known absolute mean
real(kind=dp)             :: y(size(x))
real(kind=dp)             :: ymean,y2,yvar
ymean = (c-(1/c)) * abs_x/2  ! mean of scaled data
y2    = (c**2 + 1/(c**2))/2  ! 2nd moment of scaled data
yvar  = y2 - ymean**2        ! variance of scaled data
y     = (merge(c,1/c,x>0) * x - ymean)/sqrt(yvar) ! shift and scale y to have zero mean and unit variance
end function asymm
end module asymm_mod

Here is a driver (the modules USEd are big and not provided, but the c function to evade Fortran’s string limitations is shown here).

program xasymm
! 10/24/2021 08:49 AM driver for asymm, which generates skewed data from symmetric data
use kind_mod         , only: dp
use random_normal_mod, only: random_normal
use statistics_mod   , only: print_stats
use util_mod         , only: c,pi
use asymm_mod        , only: asymm
integer      , parameter   :: n = 10**6, niter = 5
real(kind=dp), allocatable :: x(:,:)
integer                    :: iter
call random_seed()
allocate (x(n,2))
write (*,*) "#obs =",n
do iter=1,niter
   x(:,1) = random_normal(n) 
   x(:,2) = asymm(x(:,1),3.0_dp,sqrt(2/pi))
   call print_stats(c("median","mean","sd","skew","kurt","perc_05","perc_95","min","max"),x)
end do
end program xasymm

Sample output:

 #obs =     1000000

                 var     median       mean         sd       skew       kurt    perc_05    perc_95        min        max
             original     0.0004    -0.0006     1.0010    -0.0028    -0.0082    -1.6479     1.6433    -5.3463     4.6752
               skewed    -0.5743     0.0002     1.0000     1.4770     1.8428    -0.8718     2.0894    -1.5381     7.0050

                 var     median       mean         sd       skew       kurt    perc_05    perc_95        min        max
             original    -0.0001     0.0000     0.9999     0.0007     0.0009    -1.6425     1.6449    -4.5028     4.7293
               skewed    -0.5750     0.0000     1.0002     1.4812     1.8653    -0.8708     2.0920    -1.3861     7.0928

                 var     median       mean         sd       skew       kurt    perc_05    perc_95        min        max
             original     0.0006     0.0007     1.0001    -0.0000    -0.0094    -1.6463     1.6469    -4.5609     4.6480
               skewed    -0.5740     0.0010     1.0002     1.4763     1.8457    -0.8715     2.0952    -1.3966     6.9610

It can be seen that skewed variates with zero mean and unit variance are generated.

4 Likes
2

To try out the idea in What Is the Complex Step Approximation? by Nick Higham I wrote a program to compute the derivative of exp(x) at x = 0 using the complex step and the usual two-point finite difference approximation.

program xcomplex_diff
! illustrate numerical differentiation with a complex step approximation
implicit none
integer, parameter :: wp = kind(1.0)
real(kind=wp) :: h
integer :: i
write (*,"(3a14)") "h","err_cmplx","err_real"
do i=1,7
   h = 0.1_wp**i
   write (*,"(3f14.10)") h,abs(1.0_wp-[aimag(exp(cmplx(0.0_wp,h))/h),(exp(h)-exp(-h))/(2*h)])
end do
end program xcomplex_diff

output:

             h     err_cmplx      err_real
  0.1000000015  0.0016657710  0.0016676188
  0.0100000007  0.0000166893  0.0000168085
  0.0010000000  0.0000001192  0.0000169277
  0.0001000000  0.0000000000  0.0001658201
  0.0000100000  0.0000000000  0.0013579130
  0.0000010000  0.0000000000  0.0165235996
  0.0000001000  0.0000000000  0.1920926571

The complex step gives a smaller error. Do people use this method in practice to compute gradients of functions? Do you need to define a complex function of complex variables to apply the complex step, even though the function of interest is a real function of real variables?

1 Like
3

In principle, yes, but it may be possible to avoid using complex variables in the program implementation for some instances of f(x).

Take the example problem used in your link, namely, evaluating the derivative of f(x) = exp(x) at x = 1. Let us pretend that we do not know the analytical derivative of f(x), and simply use the approximate formula, f’(x) = Im(x+ih)/h. For our function, this gives f’(x) = Im exp(x+ih)/h = [exp(x) sin h]/h, so the derivative at x = 1 is e.sin h/h, and I know that sin h/h → 1 as h → 0, so the derivative is e = exp(1).

There, I have evaluated the derivative without

(i) knowing the analytical derivative of the function,

(ii) evaluating the sin function, which came from the imaginary part of exp(x + ih).

For more details of the method, including a harder example, see https://hal.archives-ouvertes.fr/hal-01483287/document .

1 Like
4

I had never heard of it, but the mathematical explanation is clear enough. What I do find puzzling in the exposé by Higham is this remark:

It is important to note that the method used to evaluate f must not itself use complex arithmetic (as methods based on the Schur decomposition do); if it does, then the interaction of those complex terms with the much smaller \mathrm{i}hE term can lead to damaging subtractive cancellation.

Any idea what is meant there? It appears after the demonstration that you can use the method for the derivative of matrix functions

1 Like
5

I think I have read about this complex trick before, it’s very clever. @Arjen I think the idea behind that remark is that when you are evaluating \exp(x+iz), then if z is very small compared to x, as it is in our case above, then evaluating the value of \exp(x+iz) must not add x and z together internally, or something like that. In other words, that perhaps there are ways of evaluating complex functions that are not accurate if the imaginary part is very small compared to the real part. Or to rephrase it, perhaps some complex function implementations have different level of accuracy in different areas of the complex plane. It seems he is saying to avoid such implementations, that the method relies on being able to accurately evaluate all complex functions everywhere in the complex plane.

6

More info about the complex numerical differentiation (with refs):

1 Like
7

I thought I’d dabble a bit with a few simple functions myself. It works well with several very simple ones I tried. But it fails to give meaningful results in this particular case:

    real(kind=kind(1.0d0))    :: h
    complex(kind=kind(1.0d0)) :: x, y
    integer :: i

    do i = 1,60
        h = 0.1d0 ** (3*i)
        x = cmplx(1.0d0,h)

        y = cos(x)/(1.0d0 - sin(x))

        write(*,*) imag( y ) / h, h
    enddo

The output is:

   6.3079749720536427        1.0000000000000002E-003
   6.3079935004987480        1.0000000000000006E-006
   6.3079933380414968        1.0000000000000009E-009
   6.3079934912382347        1.0000000000000014E-012
   6.3079935393259392        1.0000000000000017E-015
   6.3079938054362819        1.0000000000000018E-018
   6.3079933162628521        1.0000000000000024E-021
   6.3079936397112633        1.0000000000000025E-024
   6.3079936883062775        1.0000000000000027E-027
   6.3079935364468538        1.0000000000000031E-030
   6.3079936662095397        1.0000000000000035E-033
   6.3079937431475148        1.0000000000000040E-036
   6.3079948745883598        1.0000000000000042E-039
   6.3113184820803223        1.0000000000000045E-042
   8.8393816275634745        1.0000000000000050E-045
  -0.0000000000000000        1.0000000000000051E-048
  -0.0000000000000000        1.0000000000000053E-051
  -0.0000000000000000        1.0000000000000058E-054
  -0.0000000000000000        1.0000000000000060E-057
  ...

So it definitely is not a panacea

8

I tried it with gfortran 10.2.0 and Intel Fortran oneAPI 2021.4.0

9

The problem is with this statement. The expression to the right of ‘=’ is of type default real. To make it double precision, replace the line by

x = cmplx(1.0d0,h,kind(1d0))
1 Like
10

Before the Fortran 2008 standard, elemental procedures were a subset of pure ones. Fortran 2008 allows elemental functions to be further declared pure or impure. I have not used impure elemental procedures in production codes, but I can see some applications:

(1) For a derived type I often define a display subroutine interface with underlying module procedures display_scalar and display_vec. Instead I can just define a single impure elemental subroutine display, since call display([dt1,dt2]) has the effect of
call display(dt1)
call display(dt2)

(2) Sometimes you want to write the same data to multiple units, but write ([unit1, unit2],*) foo is not legal – the unit number must be a scalar. However, the write statements can be placed in an impure elemental subroutine with a unit number argument, and that subroutine can then be called with an array of unit numbers. You can get the effect of the Unix tee command, which writes to standard output and one or more files. As usual, @FortranFan had already thought of this .

(3) If you want to write each derived type to a separate file, or read each derived type from a separate file, that can be done with call foo([dt1,dt2],[unit1,unit2])

Here is an illustration.

module m
implicit none
contains
impure elemental subroutine disp(i,output_unit)
integer, intent(in) :: i,output_unit
write (output_unit,*) "i=",i
write (output_unit,*) "i**2=",i**2
end subroutine disp
end module m

program main
use m, only: disp
implicit none
call disp([4,9],6) ! call disp with argument i successively 4 and 9
call disp(3,[6,7]) ! call disp with argument output_unit successively 6 (standard output) and 7
call disp([10,11],[10,11]) ! call disp with arguments (10,10) and then (11,11)
end program main
11

Oh, yes, you are quite correct - I completely missed that bit! And that definitely explains the results.

12

My earlier program also should have specified the kind argument in the cmplx function. Note to self: do not use the real and cmplx functions without a kind argument. A version of the code using double precision shows that with the complex step you can use very small step sizes and not have problems with roundoff error, as opposed to the conventional approach with dydx = (f(x+h)-f(x-h))/(2*h)

program xcomplex_diff
! illustrate numerical differentiation with a complex step approximation
implicit none
integer, parameter :: wp = kind(1.0d0)
real(kind=wp) :: h
integer :: i
write (*,"(3a18)") "h","err_cmplx","err_real"
do i=1,14
   h = 0.1_wp**i
   write (*,"(3f18.14)") h,abs(1.0_wp-[aimag(exp(cmplx(0.0_wp,h,kind=wp))/h),(exp(h)-exp(-h))/(2*h)])
end do
end program xcomplex_diff

output:

                 h         err_cmplx          err_real
  0.10000000000000  0.00166583353172  0.00166750019844
  0.01000000000000  0.00001666658333  0.00001666674999
  0.00100000000000  0.00000016666666  0.00000016666668
  0.00010000000000  0.00000000166667  0.00000000166689
  0.00001000000000  0.00000000001667  0.00000000001210
  0.00000100000000  0.00000000000017  0.00000000002676
  0.00000010000000  0.00000000000000  0.00000000052636
  0.00000001000000  0.00000000000000  0.00000000607747
  0.00000000100000  0.00000000000000  0.00000002722922
  0.00000000010000  0.00000000000000  0.00000008274037
  0.00000000001000  0.00000000000000  0.00000008274037
  0.00000000000100  0.00000000000000  0.00003338943111
  0.00000000000010  0.00000000000000  0.00024416632505
  0.00000000000001  0.00000000000000  0.00079927783736
1 Like
13

An analysis of the total error (discretization error + error in computing the function itself) for the central difference approximation shows that the total error is minimized when h is of the order of (machine_epsilon)^(1/3). Your results illustrate that conclusion.

1 Like
14

This program simulates the birthday problem and confirms that in a group of only 23 people there is about a 50% chance that two people share a birthday.

program xbirthday_paradox
implicit none
integer, parameter :: ndays = 365, niter = 10**7
integer            :: i,iter,bdays(ndays),nsame(niter),ntimes,ntimes_cumul
real               :: xran
call random_seed()
nsame = ndays + 1
do_iter: do iter=1,niter
   do i=1,ndays
      call random_number(xran)
      bdays(i) = 1 + ndays*xran
      if (i==1) cycle
      if (any(bdays(i) == bdays(:i-1))) then
         nsame(iter) = i
         cycle do_iter
      end if
   end do
end do do_iter
ntimes_cumul = 0
write (*,"('#sim = ',i0)") niter
write (*,"(4a9)") "#kids","#times","prob","cumul"
do i=2,maxval(nsame)
   ntimes = count(nsame == i)
   ntimes_cumul = ntimes_cumul + ntimes
   if (ntimes > 0) write (*,"(2i9,2f9.6)") i,ntimes,[ntimes,ntimes_cumul]/real(niter)
end do
end program xbirthday_paradox

It takes about 6s on my PC using gfortran -O3 to give sample output below. It would be more efficient to keep bdays sorted and do a bisection search for matching birthdays. If someone extends the program to simulate how many kids are needed to have a 50% chance that n = 2 to 5 kids have the same birthday, the program can be submitted to Rosetta Code

#sim = 10000000
    #kids   #times     prob    cumul
        2    27304 0.002730 0.002730
        3    54287 0.005429 0.008159
        4    81544 0.008154 0.016314
        5   107983 0.010798 0.027112
        6   133332 0.013333 0.040445
        7   157121 0.015712 0.056157
        8   180527 0.018053 0.074210
        9   203142 0.020314 0.094524
       10   223524 0.022352 0.116876
       11   241319 0.024132 0.141008
       12   258106 0.025811 0.166819
       13   273717 0.027372 0.194191
       14   286781 0.028678 0.222869
       15   298633 0.029863 0.252732
       16   307366 0.030737 0.283469
       17   314169 0.031417 0.314885
       18   319913 0.031991 0.346877
       19   322705 0.032270 0.379147
       20   322507 0.032251 0.411398
       21   322835 0.032283 0.443682
       22   320950 0.032095 0.475776
       23   315935 0.031594 0.507370
       24   310578 0.031058 0.538428
       25   303176 0.030318 0.568745
       26   294716 0.029472 0.598217
       27   285857 0.028586 0.626803
       28   275502 0.027550 0.654353
       29   265974 0.026597 0.680950
       30   253333 0.025333 0.706284
       31   241172 0.024117 0.730401
       32   228939 0.022894 0.753295
       33   216330 0.021633 0.774928
       34   203090 0.020309 0.795237
       35   190968 0.019097 0.814333
       36   178320 0.017832 0.832165
       37   166397 0.016640 0.848805
       38   153077 0.015308 0.864113
       39   141690 0.014169 0.878282
       40   130051 0.013005 0.891287
       41   119239 0.011924 0.903211
       42   108679 0.010868 0.914079
       43    98849 0.009885 0.923964
       44    89730 0.008973 0.932937
       45    80893 0.008089 0.941026
       46    73108 0.007311 0.948337
       47    65053 0.006505 0.954842
       48    58332 0.005833 0.960675
       49    51772 0.005177 0.965852
       50    45802 0.004580 0.970433
       51    40625 0.004062 0.974495
       52    35479 0.003548 0.978043
       53    31458 0.003146 0.981189
       54    27397 0.002740 0.983929
       55    23779 0.002378 0.986306
       56    20843 0.002084 0.988391
       57    17842 0.001784 0.990175
       58    15364 0.001536 0.991711
       59    13097 0.001310 0.993021
       60    11343 0.001134 0.994155
       61     9487 0.000949 0.995104
       62     8162 0.000816 0.995920
       63     6904 0.000690 0.996611
       64     5798 0.000580 0.997190
       65     4953 0.000495 0.997686
       66     4127 0.000413 0.998098
       67     3453 0.000345 0.998444
       68     2770 0.000277 0.998721
       69     2357 0.000236 0.998957
       70     1933 0.000193 0.999150
       71     1603 0.000160 0.999310
       72     1321 0.000132 0.999442
       73     1093 0.000109 0.999551
       74      902 0.000090 0.999642
       75      700 0.000070 0.999712
       76      649 0.000065 0.999777
       77      471 0.000047 0.999824
       78      380 0.000038 0.999862
       79      323 0.000032 0.999894
       80      239 0.000024 0.999918
       81      198 0.000020 0.999938
       82      144 0.000014 0.999952
       83      105 0.000010 0.999963
       84       70 0.000007 0.999970
       85       71 0.000007 0.999977
       86       56 0.000006 0.999982
       87       40 0.000004 0.999986
       88       28 0.000003 0.999989
       89       18 0.000002 0.999991
       90       21 0.000002 0.999993
       91       11 0.000001 0.999994
       92       16 0.000002 0.999996
       93        8 0.000001 0.999996
       94       10 0.000001 0.999997
       95        9 0.000001 0.999998
       96        2 0.000000 0.999999
       97        5 0.000000 0.999999
       98        4 0.000000 1.000000
       99        4 0.000000 1.000000
      102        1 0.000000 1.000000
15

Monte Carlo simulations can save the day when analytical solutions are not known or are difficult to use, but when the latter are available they are usually far more efficient. Here is a shorter program (run time < 0.1 s) to evaluate the smallest value of n for which pbar(n) falls below 0.5.

! https://en.wikipedia.org/wiki/Birthday_problem
program bday
implicit none
double precision :: pbar
integer :: i
!
pbar = 1
do i = 1, 364
   pbar = pbar*(1d0-i/365d0)     ! Formula for pbar(n) in Wikipedia article
   if(pbar < 0.5d0)exit
end do
print 10,i+1
10 format('pbar(',i3,') is < 0.5')
end program
16

Nice!
Hypothetical birthday problem on other planets, :grin:

by_solar_system

Year in other planets (in Earth days)

Mercury     88 days
Venus      225 days
Earth      365 days
Mars       687 days
Jupiter   4333 days
Saturn   10759 days
Uranus   30687 days
Neptune  60190 days
3 Likes
17

Venus has a period of rotation of 243 Earth Days, a period of revolution of 225 Earth Days, and a surface temperature over 400 C, so an Earthling placed on Venus would not live to see the end of its first Venusian day!

1 Like
18

I used the complex step method recently to compute Jacobians for the Inviscid Euler equations.
You can see this report for more information:

https://www.researchgate.net/publication/338399112_Comparing_Matrix-based_and_Matrix-free_Discrete_Adjoint_Approaches_to_the_Euler_Equations

One of the reasons for investigating the complex step method was that it was easy to convert the original Fortran code using real to using complex by simple string substitution.
(Sadly I can’t share the Fortran code currently, but I hope to get it open source in the future.)

There are several benefits to the complex step method. Most notably, and as demonstrated above, it is second order accurate and it can achieve effectively zero truncation error for sufficiently small step sizes since it involves no differencing of similar-sized terms.

Conceptually, I find the complex step method very similar to the dual number approach using operator overloading. However, as I understand it, the dual number approach is more computationally efficient than the complex step method since it only computes terms necessary for the derivative term whereas the complex arithmetic incurs redundant computation.

2 Likes
19

I wrote a program to simulate the probability that one’s savings will last N years given a spending rate and assumptions about portfolio returns. Such simulations have been discussed in the Wall Street Journal and elsewhere.

Main Program 
program xruin
use kind_mod    , only: dp
use ziggurat_mod, only: rnor, jsrset
use stats_mod   , only: mean, median, mode
implicit none
integer, parameter :: nsteps = 40, nsim = 10000, nnfut = 4, h = nsteps/nnfut, nspend = 3, nlev = 3
real(kind=dp)      :: avg_return, sd_return, spend, wealth, wealth_init, wealth_final(nsim), &
                      spent(nsim), portfolio_return, leverage
integer            :: i, ifut, ilev, isim, ispend, tsurv(nsim), nfut(nnfut)
real(kind=dp), parameter :: xsim = real(nsim,kind=dp), leverage_min = 0.0_dp, leverage_h = 0.5_dp, &
                      spend_min = 0.02_dp, spend_h = 0.01_dp
character (len=20) :: label_fut(nnfut)
avg_return  = 0.06_dp ! average annual stock market return
sd_return   = 0.15_dp ! standard deviation of stock market returns
wealth_init = 1.00_dp
write (*,"(*(a10))") "#sim","avg_ret","sd_ret"
write (*,"(i10,2f10.2)") nsim,avg_return,sd_return
nfut = [(ifut*h,ifut=1,nnfut)] ! time horizons for which survival probabilities shown
do ifut=1,nnfut
   write (label_fut(ifut),"('p',i0)") nfut(ifut) ! labels for survival probabilities
end do
write (*,"(/,a6,10x,2a7,2a13)") "spend","spent","spent","years_surv","years_surv"
write (*,"(a6,a10,2a7,4a13,*(a7))") "rate","leverage","median","mean","median", &
       "mean","wealth_avg","wealth_surv",(trim(label_fut(ifut)),ifut=1,nnfut)
do ispend=1,nspend
   do ilev=1,nlev
      spend = spend_min + (ispend-1)*spend_h ! spending rate
      call jsrset(123456)    ! reset the seed to get the same normal variates for each level of spending and leverage
      leverage = leverage_min + (ilev-1)*leverage_h ! allocation to stock market
      spent = 0.0_dp
      do isim=1,nsim
         wealth = wealth_init
         do i=1,nsteps ! loop over time periods
            portfolio_return = leverage * (avg_return + sd_return*rnor())
            wealth = (wealth - spend) * (1 + portfolio_return) ! withdraw "spend" at beginning of period, invest the remainder
            spent(isim) = spent(isim) + spend ! cumulative amount spent
            if (wealth < spend) exit ! ran out of money
         end do
         wealth_final(isim) = wealth
         tsurv(isim) = i ! time that wealth survived
      end do
      write (*,"(f6.3,f10.4,2f7.3,4(1x,f12.4),*(1x,f6.4))") spend,leverage,median(spent),mean(spent), &
        median(tsurv),sum(tsurv)/xsim,mean(wealth_final),mean(pack(wealth_final,tsurv > nsteps)), &
        [(count(tsurv > nfut(ifut)),ifut=1,nnfut)]/xsim
   end do
   write (*,*)
end do
end program xruin
Results and Discussion 
      #sim   avg_ret    sd_ret
     10000      0.06      0.15

 spend            spent  spent   years_surv   years_surv
  rate  leverage median   mean       median         mean   wealth_avg  wealth_surv    p10    p20    p30    p40
 0.020    0.0000  0.800  0.800      41.0000      41.0000       0.2000       0.2000 1.0000 1.0000 1.0000 1.0000
 0.020    0.5000  0.800  0.800      41.0000      40.9867       1.7227       1.7285 1.0000 1.0000 0.9999 0.9966
 0.020    1.0000  0.800  0.796      41.0000      40.7936       7.1893       7.3493 1.0000 0.9987 0.9917 0.9782

 0.030    0.0000  0.990  0.990      33.0000      33.0000       0.0100       0.0000 1.0000 1.0000 1.0000 0.0000
 0.030    0.5000  1.200  1.179      41.0000      40.1650       0.9740       1.1111 1.0000 1.0000 0.9763 0.8749
 0.030    1.0000  1.200  1.168      41.0000      39.8441       5.4993       6.1338 1.0000 0.9905 0.9469 0.8963

 0.040    0.0000  0.960  0.960      24.0000      24.0000       0.0400       0.0000 1.0000 1.0000 0.0000 0.0000
 0.040    0.5000  1.600  1.419      40.0000      35.9466       0.3888       0.7836 1.0000 0.9873 0.7713 0.4828
 0.040    1.0000  1.600  1.473      41.0000      37.5720       3.9787       5.3301 0.9997 0.9539 0.8361 0.7455

Row 6 of the table above means that if the investor spends 3% of the initial portfolio value, adjusted for inflation, and invests the 100% of savings in the stock market, which has average after-inflation returns of 6% with standard deviation of 15%, that the probability of the savings lasting 30 years (column p30) is 94.7%. If the annual spending rate is 4%, the 30-year survival probability falls to 83.6%. The investor should decide whether spending 33% more per year is worth a higher risk of running out of money. Using different return assumptions, Morningstar recommends 3.3% as a safe withdrawal rate.

3 Likes
20

This is exactly right. Using complex numbers for derivatives is a cute toy and a decent hack in languages that don’t let you use custom types easily, but dual number based approaches are pretty much strictly superior.

2 Likes