 Fortran code snippets

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.

2 Likes

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

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

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

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.

1 Like

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

I tried it with gfortran 10.2.0 and Intel Fortran oneAPI 2021.4.0

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

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

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

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

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

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.

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

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

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

Nice!
Hypothetical birthday problem on other planets,  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

1 Like

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

I used the complex step method recently to compute Jacobians for the Inviscid 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

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

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