Fortran code snippets

21

Michael Wirth of the Craft of Coding blog explained that Fibonacci by linear recursion is much faster than with binary recursion, whose execution time grows exponentially with n, as illustrated by the code

module m
implicit none
integer, parameter :: ikind = selected_int_kind(13)
contains
pure recursive function fib_linear(a,b,n) result(fib)
integer (kind=ikind), intent(in) :: a,b,n
integer (kind=ikind)             :: fib
if (n <= 2) then
   fib = b
else
   fib = fib_linear(b,a+b,n-1)
end if
end function fib_linear
!
pure recursive function fib_binary(n) result(fib)
integer (kind=ikind), intent(in) :: n
integer (kind=ikind)             :: fib
if (n == 0) then
   fib = 0
else if (n == 1) then
   fib = 1
else
   fib = fib_binary(n-1) + fib_binary(n-2)
end if
end function fib_binary
end module m

program main
use m, only: ikind, fib_linear, fib_binary
implicit none
real :: t1,t2,t3
integer (kind=ikind), parameter :: n = 50_ikind
print*,"n =",n
call cpu_time(t1)
print*,"fib_linear(n) =",fib_linear(1_ikind,1_ikind,n)
call cpu_time(t2)
print*,"fib_binary(n) =",fib_binary(n)
call cpu_time(t3)
print*,"times for linear, binary recursion:",t2-t1,t3-t2
end program main

which gives output with gfortran -O3 of

 n =                   50
 fib_linear(n) =          12586269025
 fib_binary(n) =          12586269025
 times for linear, binary recursion:   0.00000000       31.7343750
22

Beatty’s theorem

Posted on 9 January 2022 by John D. Cook

Here’s a surprising theorem [1].

(Beatty’s theorem) Let a and b be any pair of irrational numbers greater than 1 with 1/ a + 1/ b = 1. Then the sequences { ⌊ na ⌋ } and { ⌊ nb ⌋ } contain every positive integer without repetition.

Fortran code (Cook used Python) 
module beatty_mod
private
public :: wp, floor_mult, check
integer , parameter :: wp = selected_real_kind(33, 4931) ! kind(1.0d0)
contains
pure function floor_mult(x,n) result(vec)
real(kind=wp), intent(in) :: x
integer      , intent(in) :: n
integer                   :: vec(n)
integer                   :: i
real(kind=wp)             :: xmult
xmult = 0.0_wp
do i=1,n
   xmult = xmult + x
   vec(i) = int(xmult)
end do
end function floor_mult
!
subroutine check(rvec,svec)
integer, intent(in) :: rvec(:),svec(:)
integer             :: i,j,n,r_old
logical, parameter  :: print_svec = .false.
n = size(rvec)
j = 0
do i=2,n
   r_old = rvec(i-1)
   if (rvec(i) /= r_old + 1) then
      j = j + 1
      if (print_svec) print*,svec(j),r_old+1
      if (svec(j) /= r_old + 1) then
         print*,i, j, r_old, rvec(i), svec(j)
         stop "did not work"
      end if
   end if 
end do
end subroutine check
end module beatty_mod
!
program test_beatty
use beatty_mod, only: wp, floor_mult,check
implicit none
real(kind=wp)    , parameter :: r = sqrt(2.0_wp), s = r/(r-1)
character (len=*), parameter :: fmt="(a,*(1x,i0))"
integer          , parameter :: n = 20 ! 10**7
logical          , parameter :: print_rs_vec = .true. ! .false.
integer                      :: rvec(n),svec(n)
rvec = floor_mult(r,n)
svec = floor_mult(s,n)
if (print_rs_vec) then
   print fmt,"{na}",rvec
   print fmt,"{nb}",svec
end if
print "('n, rvec(n) = ',i0,1x,i0)", rvec(n)
call check(rvec,svec)
end program test_beatty

Results:

{na} 1 2 4 5 7 8 9 11 12 14 15 16 18 19 21 22 24 25 26 28
{nb} 3 6 10 13 17 20 23 27 30 34 37 40 44 47 51 54 58 61 64 68
n, rvec(n) = 28

The numbers not in {na} are in {nb} and {na} and {nb} are disjoint.

Changing two lines in the main program to

integer          , parameter :: n = 10**8 ! n = 20
logical          , parameter :: print_rs_vec = .false.

the program, which uses quadruple precision, runs in a few seconds and verifies Beatty’s theorem up to that N.

1 Like
23

I wanted to fill a matrix from a vector by row, setting unfilled elements to a default value, and wrote a subroutine to do so.

The most common use of the reshape intrinsic may be to convert an array of one shape to one with the same number of elements but a different shape, but it has optional pad and order arguments that make the task mentioned above a one-liner:

program main
! demonstrate reshape
implicit none
integer          , parameter :: n1 = 3, n2 = 5, n = 8
integer                      :: v(n), m(n1,n2), i
character (len=*), parameter :: fmt_i = "(*(i4))"
v = [(2*i, i=1,n)]
print*,"v"
print fmt_i,v
m = reshape(v, shape(m), pad = [0])
call print_m("fill by column")
m = reshape(v, shape(m), pad = [0], order=[2,1])
call print_m("fill by row")
contains
subroutine print_m(title)
character (len=*), intent(in) :: title
print "(/,a)",title
do i=1,n1
   print fmt_i,m(i,:)
end do
end subroutine print_m
end program main
! output:
!  v
!    2   4   6   8  10  12  14  16
! 
! fill by column
!    2   8  14   0   0
!    4  10  16   0   0
!    6  12   0   0   0
! 
! fill by row
!    2   4   6   8  10
!   12  14  16   0   0
!    0   0   0   0   0
1 Like
24

First time sharing rather than asking for help, hope this is helpful for someone!

I want to have a good-looking terminal output for numerical result, so I wrote the boxprint subroutine:

    subroutine boxprint(str, symbol)

        character(len=*) :: symbol
        character(len=*), intent(in) :: str ! input str
        character(len=len_trim(str)+1) :: tmp ! tmp str for analysis
        character(len=1), parameter :: nl = char(10)
        integer, dimension(:), allocatable :: linelen ! newline position
        integer :: i, n, last

        allocate(linelen(len_trim(str) + 1), source = 0)
        tmp = trim(str) // nl

        last = 0
        n = 0
        do i = 1, len(tmp) + 1
            if (iachar(tmp(i:i)) == 10 .and. .not. iachar(tmp(i+1:i+1)) == 10) then
                linelen(i) = i - last - 1;
                last = i
                n = n + 1
            endif
        enddo

        write(*, '(A)') ''
        write(*, '(a)') repeat(symbol, maxval(linelen))
        write(*, '(a)') trim(str)
        write(*, '(a)') repeat(symbol, maxval(linelen))
        write(*, '(A)') ''

    end subroutine boxprint

To print any “paragraph” surrounded by symbol, you call this subroutine like the following, where sep is a character with long-enough memory allocated, and nl is the newline character:

        write(sep, '(a, I0, 2a, 2(a, F9.6), 2a, 3(a, F9.6), 2a, 2(a, F9.6))') &
                "Bisection Ends: iter = ", iter, &
                nl, spaces, &
                "norm(E) = ", norm(evalE), &
                ", longest = ", longest, &
                nl, spaces, &
                'w = ', wval, &
                ', q = ', qsell, &
                ', Q = ', Qbuy, &
                nl, spaces, &
                'w - psi*cagg = ', evalE(1), &
                ', usedInv - disInv = ', evalE(2)
        call boxprint(sep, '=')

With this calling, you will get output like this:

==========================================================
Bisection Ends: iter = 34
    norm(E) =  0.000010, longest =  0.000004
    w =  1.048213, q =  0.955162, Q =  0.993626
    w - psi*cagg = -0.000008, usedInv - disInv = -0.000006
==========================================================

Hope this is helpful!

2 Likes
25

Just a remark: You may allow for an optional argument to write to a particular unit (opened file) instead of only to screen (standard output)

26

@Arjen That’s a good idea!

27

Here is a program that computed 1000 digits of pi in 0.03s on my PC in 14 lines of FORTRAN:

      integer vect(3350), buffer(201)
      data vect/3350*2/, more/0/
      DO 2 n=1, 201
         karray = 0
         DO 3 L=3350, 1, -1
            num = 100000*vect(L) + karry*L
            karry = num/(2*L-1)
3        vect(L) = num - karry*(2*L-1)
         k = karry/100000
         buffer(n) = more+k
2     more = karry - k*100000
      write(*,100) buffer
100   format(1x,I1,'.'/(1x,10I5.5))
      end

or Fortran (using a variable m instead of l)

implicit none
integer vect(3350), buffer(201), m, n, karry, k, more, num
vect = 2
more = 0
do n=1, 201
   karry = 0
   do m=3350, 1, -1
      num = 100000*vect(m) + karry*m
      karry = num/(2*m-1)
      vect(m) = num - karry*(2*m-1)
   end do
   k = karry/100000
   buffer(n) = more + k
   more = karry - k*100000
end do
write (*,"(1x,I1,'.'/(1x,10I5.5))") buffer
end

with the theory explained in

A Spigot Algorithm for the Digits of π
by Stanley Rabinowitz and Stan Wagon
The American Mathematical Monthly, Vol. 102, No. 3 (Mar., 1995), pp. 195-203

giving output

 3.
 14159265358979323846264338327950288419716939937510
 58209749445923078164062862089986280348253421170679
 82148086513282306647093844609550582231725359408128
 48111745028410270193852110555964462294895493038196
 44288109756659334461284756482337867831652712019091
 45648566923460348610454326648213393607260249141273
 72458700660631558817488152092096282925409171536436
 78925903600113305305488204665213841469519415116094
 33057270365759591953092186117381932611793105118548
 07446237996274956735188575272489122793818301194912
 98336733624406566430860213949463952247371907021798
 60943702770539217176293176752384674818467669405132
 00056812714526356082778577134275778960917363717872
 14684409012249534301465495853710507922796892589235
 42019956112129021960864034418159813629774771309960
 51870721134999999837297804995105973173281609631859
 50244594553469083026425223082533446850352619311881
 71010003137838752886587533208381420617177669147303
 59825349042875546873115956286388235378759375195778
 18577805321712268066130019278766111959092164201989
3 Likes
28

password generator:

module rnd
  implicit none
  private
  public pwd
  contains
    function pwd(n) result(r)
      ! return a printable random string n chars long
      ! spaces are not considered i.e. only ascii 33 - 126
      ! be sure to call random_seed() in the main prog before calling this
      integer i, n
      character(n) r
      real h(n)
      call random_number(h)
      do i = 1, n
        r(i:i) = achar(33+floor(94*h(i)))
      end do
    end function pwd
end module rnd
29

A minor suggestion might be to use achar() here instead. Most compilers default to ascii anyway, so there may not be much practical difference, but who knows what might happen in the future with all of the fancy character encoding that are available.

1 Like
30

Here is an example of using R within a Fortran program to plot and describe data. When the plot window is closed, the program resumes. Admittedly CSV is not the right format for large data sets.

program xtrig
! write data to file, call R to plot it
implicit none
integer, parameter :: dp = kind(1.0d0), n = 10**4 + 1
real(kind=dp), parameter :: pi = 4*atan(1.0_dp), xk = 4.0_dp
real(kind=dp) :: x(n), y(n, 2)
integer :: i, outu
character (len=*), parameter :: out_file = "data.csv"
open (unit=outu, file=out_file, action="write", status="replace")
write (outu,"(*(a,:,','))") "x","sin(x)","cos(x)"
do i=1,n
   x(i) = xk*pi*(i-1)/(n-1.0_dp)
   y(i,:) = [sin(x(i)), cos(x(i))]
   write (outu,"(*(f0.6,:,','))") x(i), y(i,:)
end do
print*,"calling R, n =",n
print*
call execute_command_line("Rscript xplot_csv.r " // trim(out_file))
end program xtrig

R script xplot_csv.r

args <- commandArgs(trailingOnly = TRUE)
if (length(args) == 0) stop("You must specify a filename!")
filename <- args[1]
if(!interactive()) windows()
isWindowOpen <- function() {!is.null(dev.list()["windows"])}

plot_vs_first_column <- function(data) {
  # plot columns from the 2nd on against the 1st column
  num_columns <- ncol(data)
  # Determine y-axis limits
  y_limits <- c(min(sapply(data[-1], min)), max(sapply(data[-1], max)))
  # Plot the first data column against each of the others
  plot(data[[1]], data[[2]], type="l", col=2, xlab=names(data[1]), ylab="",
	  ylim=y_limits, main="plot")
  for(i in 3:num_columns) lines(data[[1]], data[[i]], col=i)
  legend("topright", legend=names(data)[-1], fill=2:num_columns, cex=0.8) # add a legend
}

data = read.csv(filename, header=TRUE)
print(summary(data))
plot_vs_first_column(data)
while(isWindowOpen()) Sys.sleep(1)
31

How many ways can you change N dollars into coins (cents, dimes, nickels, quarters, 50 cent and dollar coins)?

Here’s a nice compact solution using Fortran’s recursion feature.

program change_for_dollars
implicit none
integer ndollars,ncents

write(*,'("Enter number of dollars : ")',advance='NO')
read(*,*) ndollars
! convert to cents
ncents=ndollars*100

write(*,'("Number of ways ",I4," dollars can be changed to coins is")') ndollars
write(*,'(I8)') nways(ncents,6)

contains

recursive function nways(ncents,ncoins) result(nw)
implicit none
! arguments
integer, intent(in) :: ncents,ncoins
integer nw
! local variables
integer c,i
! cent, nickel, dime, quarter, 50 cent and dollar coin
integer, parameter :: coin(6)=[1,5,10,25,50,100]
if (ncoins.eq.0) then
  nw=0
  return
end if
nw=nways(ncents,ncoins-1)
c=coin(ncoins)
if (ncents.lt.c) return
do i=1,ncents/c
  nw=nw+nways(ncents-i*c,ncoins-1)
end do
if (mod(ncents,c).eq.0) nw=nw+1
end function

end program
2 Likes
32

An article from Quanta Magazine A New Generation of Mathematicians Pushes Prime Number Barriers says that the remainders of prime numbers when divided by another prime number are close to evenly distributed. The following program, which runs in a few seconds, verifies this.

program primes_and_remainders
    use iso_fortran_env, only: int64
    implicit none
    integer(kind=int64), parameter :: n = 10_int64**8, kmax = 10
    logical :: sieve(n)
    integer(kind=int64) :: primes(n/2), nprimes, i, j, k, remainder, freq(n)
    print*,"primes up to", n
    ! Initialize sieve for the Sieve of Eratosthenes
    sieve = .true.
    sieve(1) = .false.
    nprimes = 0
    ! Compute primes up to n using the Sieve of Eratosthenes
    do i = 2, n
        if (sieve(i)) then
            nprimes = nprimes + 1
            primes(nprimes) = i
            j = i * i
            do while (j <= n)
                sieve(j) = .false.
                j = j + i
            end do
        end if
    end do
    print*,"# of primes, last prime:", nprimes, primes(nprimes)
    ! For each prime number k up to n, compute frequency of remainders for divisor k of primes below n
    do i = 1, nprimes
        k = primes(i)
        if (k > kmax) exit
        if (k < 3) cycle
        freq = 0
        do j = 1, nprimes
            if (primes(j) < n) then
                remainder = mod(primes(j), k)
                freq(remainder+1) = freq(remainder+1) + 1
            end if
        end do
        ! Print the frequencies for the prime k
        print "(/,a,i0)", "Prime divisor: ", k
        do j = 1, k-1
            print *, "Remainder ", j, ":", freq(j+1)
        end do
    end do
end program primes_and_remainders

output:

 primes up to            100000000
 # of primes, last prime:              5761455             99999989

Prime divisor: 3
 Remainder                     1 :              2880517
 Remainder                     2 :              2880937

Prime divisor: 5
 Remainder                     1 :              1440298
 Remainder                     2 :              1440496
 Remainder                     3 :              1440474
 Remainder                     4 :              1440186

Prime divisor: 7
 Remainder                     1 :               960023
 Remainder                     2 :               960114
 Remainder                     3 :               960213
 Remainder                     4 :               960085
 Remainder                     5 :               960379
 Remainder                     6 :               960640
1 Like
33

Here is a module to pretty print a 2D arrays on the console. It also supports 1D arrays.

    module mod_show_matrix
    use mod_common
    implicit none
    
    interface show
        procedure show_vector_i, show_vector_r, show_vector_d
        procedure show_matrix_i, show_matrix_r, show_matrix_d
    end interface
    
    contains
    
        subroutine show_vector_i(v, w)
    ! Display the vector 'v' in a single column
    !   v : the array of real numbers
    !   w : the column width. default = 5
    !   s : sig. figures w-5 (calculated)
        integer, intent(in) :: v(:)
        integer, intent(in), optional :: w
        integer :: i,n,wt
        character(len=16) :: fmt
        if(present(w)) then
            wt = w
        else
            wt = 5
        end if
        n = size(v)
        write( fmt, "(a,g0,a)") "(*(1x,g",wt,".0))"
        write( * , fmt ) ( v(i), new_line("A"), i=1,n )    
    end subroutine
       
    subroutine show_vector_r(v, w)
    ! Display the vector 'v' in a single column
    !   v : the array of real numbers
    !   w : the column width. default = 12
    !   s : sig. figures w-5 (calculated)
        real(real32), intent(in) :: v(:)
        integer, intent(in), optional :: w
        integer :: i,n,dg,wt
        character(len=16) :: fmt
        if(present(w)) then
            wt = w
        else
            wt = 12
        end if
        dg = wt - 6
        n = size(v)
        write( fmt, "(a,g0,a,g0,a)") "(*(1x,g",wt,".",dg,"))"
        write( * , fmt ) ( v(i), new_line("A"), i=1,n )    
    end subroutine
    
    subroutine show_vector_d(v, w)
    ! Display the vector 'v' in a single column
    !   v : the array of real numbers
    !   w : the column width. default = 12
    !   s : sig. figures w-5 (calculated)
        real(real64), intent(in) :: v(:)
        real(real32), allocatable :: u(:)
        integer, intent(in), optional :: w
        u =real(v)
        where( abs(u)<1e-11 )
            u = 0.0
        end where
        call show_vector_r(u,w)
    end subroutine
    
    subroutine show_matrix_i(A, w)
    ! Display the matrix 'A' in columns
    !   A : the array of integers
    !   w : the column width. default = 5
        integer, intent(in) :: A(:,:)
        integer, intent(in), optional :: w
        integer :: i,j,n,m, wt
        character(len=16) :: fmt
        if(present(w)) then
            wt = w
        else
            wt = 5
        end if
        n = size(A,1)
        m = size(A,2)
        write( fmt, "(a,g0,a)") "(*(1x,g",wt,".0))"        
        write( * , fmt ) ( (A(i,j),j=1,m), new_line("A"), i=1,n )
    end subroutine
    
    subroutine show_matrix_r(A, w)
    ! Display the matrix 'A' in columns
    !   A : the array of real numbers
    !   w : the column width. default = 12
    !   s : sig. figures w-5 (calculated)
        real(real32), intent(in) :: A(:,:)
        integer, intent(in), optional :: w
        integer :: i,j,n,m,dg,wt
        character(len=16) :: fmt
        if(present(w)) then
            wt = w
        else
            wt = 12
        end if
        dg = wt - 6
        n = size(A,1)
        m = size(A,2)
        write( fmt, "(a,g0,a,g0,a)") "(*(1x,g",wt,".",dg,"))"
        write( * , fmt ) ( (A(i,j),j=1,m), new_line("A"), i=1,n )
    end subroutine
    
    subroutine show_matrix_d(A,w)
    ! Display the matrix 'A' in columns
    !   A : the array of dble numbers
    !   w : the column width. default = 12
    ! Converts 'A' into single precision and calls `show_matrix_r`
        real(real64), intent(in) :: A(:,:)
        real(real32), allocatable :: B(:,:)
        integer, intent(in), optional :: w
        B = real(A)
        where( abs(B)<1e-11 )
            B = 0.0
        end where
        call show_matrix_r(B,w)
    end subroutine
    
    end module

example code for testing, with integer, real and double vectors and matrices

    subroutine test_mod_show()
    use mod_show_matrix
    
    integer, parameter :: n = 12, m = 6
    
    integer :: iA(n,m), iV(n)
    real :: rA(n,m), rV(n)
    real(real64) :: dA(n,m), dV(n)
    
    call RANDOM_NUMBER(dV)    
    call RANDOM_NUMBER(dA)
    
    dV = -1000.0d0 + 2000.0d0 * dV
    dA = -1000.0d0 + 2000.0d0 * dA
    
    rV = REAL(dV)
    iV = NINT(dV)
    rA = REAL(dA)
    iA = NINT(dA)
    
    print *, "Vectors"
    call show(iV)
    call show(rV)
    call show(dV)
    
    print *, "Matrices"
    call show(iA)
    call show(rA)
    call show(dA)
    
    end subroutine

with output

 Vectors
  -812
  -300
  -725
   147
   439
   982
    92
   185
   420
   815
  -457
   186

 -811.931
 -300.402
 -724.711
  147.355
  439.164
  981.604
  91.8634
  184.538
  420.419
  815.408
 -456.781
  185.557

 -811.931
 -300.402
 -724.711
  147.355
  439.164
  981.604
  91.8634
  184.538
  420.419
  815.408
 -456.781
  185.557

 Matrices
   -63   667   699  -127  -706  -993
   771   -26  -703   395   815  -193
   417   347   648   618  -980   282
   418  -996  -201  -958   382   133
   863   546  -105   -85  -732  1000
   258   -19    71    54  -659   929
   733   763  -131  -150  -211   740
  -416   490  -259   458   -35  -244
   230  -285   773   999   481   797
   498  -204    67  -495  -629  -209
     7  -896   226   619   260  -364
  -451   886  -399   574  -569   130

 -62.7683      667.230      698.727     -127.294     -706.111     -992.679
  771.382     -26.1508     -703.281      394.592      815.305     -193.117
  416.630      346.530      647.946      618.252     -979.542      281.533
  418.474     -995.676     -200.940     -957.687      382.064      132.730
  862.885      546.431     -105.299     -85.1874     -732.206      999.754
  258.189     -18.6606      71.2264      54.2078     -658.914      928.694
  733.123      762.584     -130.860     -150.476     -211.242      739.778
 -415.526      490.129     -259.413      458.103     -35.1639     -244.046
  229.853     -284.816      772.724      998.592      481.294      796.775
  497.563     -203.624      66.5493     -495.422     -628.974     -209.398
  6.86585     -896.292      226.299      619.395      259.609     -364.034
 -450.501      885.879     -399.171      574.497     -569.222      129.725

 -62.7683      667.230      698.727     -127.294     -706.111     -992.679
  771.382     -26.1508     -703.281      394.592      815.305     -193.117
  416.630      346.530      647.946      618.252     -979.542      281.533
  418.474     -995.676     -200.940     -957.687      382.064      132.730
  862.885      546.431     -105.299     -85.1874     -732.206      999.754
  258.189     -18.6606      71.2264      54.2078     -658.914      928.694
  733.123      762.584     -130.860     -150.476     -211.242      739.778
 -415.526      490.129     -259.413      458.103     -35.1639     -244.046
  229.853     -284.816      772.724      998.592      481.294      796.775
  497.563     -203.624      66.5493     -495.422     -628.974     -209.398
  6.86585     -896.292      226.299      619.395      259.609     -364.034
 -450.501      885.879     -399.171      574.497     -569.222      129.725
2 Likes
34

Thanks. Testing it with

program xshow
use show_mod
use kind_mod, only: dp
implicit none
integer, parameter :: n1=5, n2=3
real(kind=dp)      :: x(n1, n2)
call random_seed()
call random_number(x)
call show(x)
end program xshow

I got sample output

0.454561    0.694803    0.750357               
0.903481E-010.520328    0.695920               
0.551235    0.509216E-010.886892               
0.326106    0.876039    0.466012               
0.172485    0.413229    0.484142

Ideally there is always at least one space between the printed numbers, which can be accomplished by replacing

    write( fmt, "(a,g0,a,g0,a)") "(*(g",wt,".",dg,"))"`

in subroutine show_matrix_r with

    write( fmt, "(a,g0,a,g0,a)") "(*(1x,g",wt,".",dg,"))"
1 Like
35

You are correct. I usually do force a space between columns, but in the code above I did not. I think the logic being that you control the spacing with the w argument if needed.

You can try to modify w=12to w=13 to the same effect.

In any case, I have updated my post to reflect your suggestion,

1 Like
36

with integer matrix: error comes!!

program xshow
   use show_mod
   implicit none
   integer, parameter:: sp = kind(0.d0)
   integer, parameter:: dp = kind(1.d0)
   integer, parameter :: n1 = 5, n2 = 3
   integer:: i
   real(kind=dp)      :: x(n1, n2)
   integer(kind=sp)   :: a(n1, n2)
   print*, ''
   print*, '3: '
   a = reshape( [(i, i=1,15)] , [5,3])
   print*, a
   call show(a)
!output: 
!|    call show(a)
!      |               1
!Error: There is no specific subroutine for the generic ‘show’ at (1)
37

The sp parameter is a real kind. You are defining an integer variable using a real kind.

1 Like
38

thanks, I missed it.

39

I have tested the code with native integers. Look at the updated code in my post and try again.

The issue is that you are using integer(kind=sp) where sp refers to a kind for real numbers. Try again with integer or integer(kind=4), or integer(int32)

I have not included a version for 8-byte integers as I do not use them often, but the show_matrix_i method can be used as a template to write a sub for this.

1 Like
40

Fortran code for direct inversion of 2×2, 3×3 and 4×4 matrices (native arrays). Actually the generic methods defined are

  • det(A) get the determinant of a matrix
  • inv(A) invert a matrtix
  • solve(A,b) solve the system A*x=b for b
module mod_array_inv
use mod_common
implicit none

    integer, parameter :: wp = real64
    
    interface det
        module procedure :: mat_det
    end interface
    interface inv
        module procedure :: mat_inv
    end interface
    interface solve
        module procedure :: mat_solve_vec, mat_solve_mat
    end interface

    contains
            
    function mat_det(A) result(d)
    real(wp), intent(in) :: A(:,:)
    real(wp) :: d
    integer :: n, m
    
        n = size(A, 1)
        m = size(A, 2)
        
        if(n /= m) then
            error stop "Expecting square matrix."
        end if
        
        select case(n)
        case(1)  
            d = A(1,1)
        case(2)
            d = mat2_det(A)
        case(3)
            d = mat3_det(A)
        case(4)
            d = mat4_det(A)
        case default
            error stop "Size is not supported."
        end select
    end function
    
    function mat_inv(A) result(B)
    real(wp), intent(in) :: A(:,:)    
    real(wp) :: B(size(A,1),size(A,2))
    integer :: n, m
    
        n = size(A, 1)
        m = size(A, 2)
        
        if(n /= m) then
            error stop "Expecting square matrix."
        end if
        
        select case(n)
        case (1)
            B = 1/A
        case (2)
            B = mat2_inv(A)
        case (3)
            B = mat3_inv(A)
        case (4)
            B = mat4_inv(A)
        case default
            error stop "Size is not supported."
        end select
    end function
    
    function mat_solve_vec(A, b) result(x)
    real(wp), intent(in) :: A(:,:), b(:)
    real(wp), allocatable :: x(:)
    integer :: n,m
    
        n = size(A, 1)
        m = size(A, 2)
        
        if(n /= m) then
            error stop "Expecting square matrix."
        end if

        select case(n)
        case (1)
            x = b/A(1,1)
        case (2)   
            x = mat2_solve_vec(A,b)
        case (3)
            x = mat3_solve_vec(A,b)
        case (4)
            x = mat4_solve_vec(A,b)
        case default
            error stop "Size is not supported."
        end select
    end function
    
    function mat_solve_mat(A, B) result(X)
    real(wp), intent(in) :: A(:,:), B(:,:)
    real(wp), allocatable :: X(:,:)
    real(wp), allocatable :: A_inv(:,:)
        A_inv = mat_inv(A)
        X = matmul(A_inv, B)
    end function
    
    pure function mat2_det(A) result(d)
    implicit real(wp) (T)
    real(wp) :: d
    real(wp), intent(in) :: A(2,2)
    real(wp) :: t2, t5
        t2 = A(1,1)*A(2,2)
        t5 = A(1,2)*A(2,1)
        d = t2-t5
    end function
    
    pure function mat2_inv(A) result(B)
    real(wp) :: B(2,2), d_inv
    real(wp), intent(in) :: A(2,2)        
        d_inv = 1/mat2_det(A)
        B(1,1) = A(2,2)*d_inv
        B(1,2) = -A(1,2)*d_inv
        B(2,1) = -A(2,1)*d_inv
        B(2,2) = A(1,1)*d_inv        
    end function
    
    pure function mat2_solve_vec(A,b) result(x)
    real(wp) :: x(2), d_inv
    real(wp), intent(in) :: A(2,2), b(2)
        d_inv = 1/mat2_det(A)
        x(1) = d_inv*(A(2,2)*b(1) - A(1,2)*b(2))
        x(2) = d_inv*(A(1,1)*b(2) - A(2,1)*b(1))    
    end function
    
    pure function mat3_det(A) result(d)
    implicit real(wp) (T)
    real(wp) :: d
    real(wp), intent(in) :: A(3,3)
    real(wp) :: t2, t3, t4, t7, t8, t9
        t2 = A(1,1)*A(2,2)*A(3,3)
        t3 = A(1,2)*A(2,3)*A(3,1)
        t4 = A(1,3)*A(2,1)*A(3,2)
        t7 = A(1,1)*A(2,3)*A(3,2)
        t8 = A(1,2)*A(2,1)*A(3,3)
        t9 = A(1,3)*A(2,2)*A(3,1)
        d = t2+t3+t4-t7-t8-t9
    end function
    
    pure function mat3_inv(A) result(B)
    real(wp) :: B(3,3), d_inv
    real(wp), intent(in) :: A(3,3)
    
        d_inv = 1/mat3_det(A)
        B(1,1) = d_inv*(A(2,2)*A(3,3)-A(2,3)*A(3,2))
        B(1,2) = -d_inv*(A(1,2)*A(3,3)-A(1,3)*A(3,2))
        B(1,3) = d_inv*(A(1,2)*A(2,3)-A(1,3)*A(2,2))
        B(2,1) = -d_inv*(A(2,1)*A(3,3)-A(2,3)*A(3,1))
        B(2,2) = d_inv*(A(1,1)*A(3,3)-A(1,3)*A(3,1))
        B(2,3) = -d_inv*(A(1,1)*A(2,3)-A(1,3)*A(2,1))
        B(3,1) = d_inv*(A(2,1)*A(3,2)-A(2,2)*A(3,1))
        B(3,2) = -d_inv*(A(1,1)*A(3,2)-A(1,2)*A(3,1))
        B(3,3) = d_inv*(A(1,1)*A(2,2)-A(1,2)*A(2,1))    
        
    end function
    
    pure function mat3_solve_vec(A,b) result(x)
    real(wp) :: x(3), d_inv
    real(wp), intent(in) :: A(3,3), b(3)
        d_inv = 1/mat3_det(A)
        
        x(1) = d_inv*(A(1,2)*(A(2,3)*b(3)-A(3,3)*b(2))+A(1,3)*(A(3,2)*b(2)-A(2,2)*b(3))+b(1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2)))
        x(2) = d_inv*(A(1,1)*(A(3,3)*b(2)-A(2,3)*b(3))+A(1,3)*(A(2,1)*b(3)-A(3,1)*b(2))-b(1)*(A(2,1)*A(3,3)-A(2,3)*A(3,1)))
        x(3) = d_inv*(A(1,1)*(A(2,2)*b(3)-A(3,2)*b(2))+A(1,2)*(A(3,1)*b(2)-A(2,1)*b(3))+b(1)*(A(2,1)*A(3,2)-A(2,2)*A(3,1)))
    end function
    
    pure function mat4_det(A) result(d)
    implicit real(wp) (T)
    real(wp) :: d
    real(wp), intent(in) :: A(4,4)
    real(wp) :: t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15
    real(wp) :: t16, t17, t18, t19, t20, t21, t22, t23, t24, t25, t26, t27
      t2 = A(1,1)*A(2,2)*A(3,3)*A(4,4)
      t3 = A(1,1)*A(2,3)*A(3,4)*A(4,2)
      t4 = A(1,1)*A(2,4)*A(3,2)*A(4,3)
      t5 = A(1,2)*A(2,1)*A(3,4)*A(4,3)
      t6 = A(1,2)*A(2,3)*A(3,1)*A(4,4)
      t7 = A(1,2)*A(2,4)*A(3,3)*A(4,1)
      t8 = A(1,3)*A(2,1)*A(3,2)*A(4,4)
      t9 = A(1,3)*A(2,2)*A(3,4)*A(4,1)
      t10 = A(1,3)*A(2,4)*A(3,1)*A(4,2)
      t11 = A(1,4)*A(2,1)*A(3,3)*A(4,2)
      t12 = A(1,4)*A(2,2)*A(3,1)*A(4,3)
      t13 = A(1,4)*A(2,3)*A(3,2)*A(4,1)
      t16 = A(1,1)*A(2,2)*A(3,4)*A(4,3)
      t17 = A(1,1)*A(2,3)*A(3,2)*A(4,4)
      t18 = A(1,1)*A(2,4)*A(3,3)*A(4,2)
      t19 = A(1,2)*A(2,1)*A(3,3)*A(4,4)
      t20 = A(1,2)*A(2,3)*A(3,4)*A(4,1)
      t21 = A(1,2)*A(2,4)*A(3,1)*A(4,3)
      t22 = A(1,3)*A(2,1)*A(3,4)*A(4,2)
      t23 = A(1,3)*A(2,2)*A(3,1)*A(4,4)
      t24 = A(1,3)*A(2,4)*A(3,2)*A(4,1)
      t25 = A(1,4)*A(2,1)*A(3,2)*A(4,3)
      t26 = A(1,4)*A(2,2)*A(3,3)*A(4,1)
      t27 = A(1,4)*A(2,3)*A(3,1)*A(4,2)
      d = t2+t3+t4+t5+t6+t7+t8+t9+t10+t11+t12+t13-t16-t17-t18-t19-t20-t21-t22-t23-t24-t25-t26-t27
    end function    
    
    pure function mat4_inv(A) result(B)
    real(wp) :: B(4,4), d_inv
    real(wp), intent(in) :: A(4,4)
    
      d_inv = 1/mat4_det(A)
      B(1,1) = d_inv*(A(2,2)*A(3,3)*A(4,4)-A(2,2)*A(3,4)*A(4,3)-A(2,3)*A(3,2)*A(4,4)+A(2,3)*A(3,4)*A(4,2)+A(2,4)*A(3,2)*A(4,3)-A(2,4)*A(3,3)*A(4,2))
      B(1,2) = -d_inv*(A(1,2)*A(3,3)*A(4,4)-A(1,2)*A(3,4)*A(4,3)-A(1,3)*A(3,2)*A(4,4)+A(1,3)*A(3,4)*A(4,2)+A(1,4)*A(3,2)*A(4,3)-A(1,4)*A(3,3)*A(4,2))
      B(1,3) = d_inv*(A(1,2)*A(2,3)*A(4,4)-A(1,2)*A(2,4)*A(4,3)-A(1,3)*A(2,2)*A(4,4)+A(1,3)*A(2,4)*A(4,2)+A(1,4)*A(2,2)*A(4,3)-A(1,4)*A(2,3)*A(4,2))
      B(1,4) = -d_inv*(A(1,2)*A(2,3)*A(3,4)-A(1,2)*A(2,4)*A(3,3)-A(1,3)*A(2,2)*A(3,4)+A(1,3)*A(2,4)*A(3,2)+A(1,4)*A(2,2)*A(3,3)-A(1,4)*A(2,3)*A(3,2))
      B(2,1) = -d_inv*(A(2,1)*A(3,3)*A(4,4)-A(2,1)*A(3,4)*A(4,3)-A(2,3)*A(3,1)*A(4,4)+A(2,3)*A(3,4)*A(4,1)+A(2,4)*A(3,1)*A(4,3)-A(2,4)*A(3,3)*A(4,1))
      B(2,2) = d_inv*(A(1,1)*A(3,3)*A(4,4)-A(1,1)*A(3,4)*A(4,3)-A(1,3)*A(3,1)*A(4,4)+A(1,3)*A(3,4)*A(4,1)+A(1,4)*A(3,1)*A(4,3)-A(1,4)*A(3,3)*A(4,1))
      B(2,3) = -d_inv*(A(1,1)*A(2,3)*A(4,4)-A(1,1)*A(2,4)*A(4,3)-A(1,3)*A(2,1)*A(4,4)+A(1,3)*A(2,4)*A(4,1)+A(1,4)*A(2,1)*A(4,3)-A(1,4)*A(2,3)*A(4,1))
      B(2,4) = d_inv*(A(1,1)*A(2,3)*A(3,4)-A(1,1)*A(2,4)*A(3,3)-A(1,3)*A(2,1)*A(3,4)+A(1,3)*A(2,4)*A(3,1)+A(1,4)*A(2,1)*A(3,3)-A(1,4)*A(2,3)*A(3,1))
      B(3,1) = d_inv*(A(2,1)*A(3,2)*A(4,4)-A(2,1)*A(3,4)*A(4,2)-A(2,2)*A(3,1)*A(4,4)+A(2,2)*A(3,4)*A(4,1)+A(2,4)*A(3,1)*A(4,2)-A(2,4)*A(3,2)*A(4,1))
      B(3,2) = -d_inv*(A(1,1)*A(3,2)*A(4,4)-A(1,1)*A(3,4)*A(4,2)-A(1,2)*A(3,1)*A(4,4)+A(1,2)*A(3,4)*A(4,1)+A(1,4)*A(3,1)*A(4,2)-A(1,4)*A(3,2)*A(4,1))
      B(3,3) = d_inv*(A(1,1)*A(2,2)*A(4,4)-A(1,1)*A(2,4)*A(4,2)-A(1,2)*A(2,1)*A(4,4)+A(1,2)*A(2,4)*A(4,1)+A(1,4)*A(2,1)*A(4,2)-A(1,4)*A(2,2)*A(4,1))
      B(3,4) = -d_inv*(A(1,1)*A(2,2)*A(3,4)-A(1,1)*A(2,4)*A(3,2)-A(1,2)*A(2,1)*A(3,4)+A(1,2)*A(2,4)*A(3,1)+A(1,4)*A(2,1)*A(3,2)-A(1,4)*A(2,2)*A(3,1))
      B(4,1) = -d_inv*(A(2,1)*A(3,2)*A(4,3)-A(2,1)*A(3,3)*A(4,2)-A(2,2)*A(3,1)*A(4,3)+A(2,2)*A(3,3)*A(4,1)+A(2,3)*A(3,1)*A(4,2)-A(2,3)*A(3,2)*A(4,1))
      B(4,2) = d_inv*(A(1,1)*A(3,2)*A(4,3)-A(1,1)*A(3,3)*A(4,2)-A(1,2)*A(3,1)*A(4,3)+A(1,2)*A(3,3)*A(4,1)+A(1,3)*A(3,1)*A(4,2)-A(1,3)*A(3,2)*A(4,1))
      B(4,3) = -d_inv*(A(1,1)*A(2,2)*A(4,3)-A(1,1)*A(2,3)*A(4,2)-A(1,2)*A(2,1)*A(4,3)+A(1,2)*A(2,3)*A(4,1)+A(1,3)*A(2,1)*A(4,2)-A(1,3)*A(2,2)*A(4,1))
      B(4,4) = d_inv*(A(1,1)*A(2,2)*A(3,3)-A(1,1)*A(2,3)*A(3,2)-A(1,2)*A(2,1)*A(3,3)+A(1,2)*A(2,3)*A(3,1)+A(1,3)*A(2,1)*A(3,2)-A(1,3)*A(2,2)*A(3,1))

    end function

    pure function mat4_solve_vec(A,b) result(x)
    real(wp) :: x(4)
    real(wp), intent(in) :: A(4,4), b(4)
        x = matmul(mat4_inv(A), b)
    end function
    
    
end module

and here is some sample output for n=1 all the way up to n=4

 Size =            1
 A=
 0.605321

 det(A) =   0.605321032671391
 inv(A)=
  1.65202

 b=
 x=
  1.55042

 max(error)=  0.000000000000000E+000
 Size =            2
 A=
 0.837911     0.753904
 0.262199E-01 0.359381

 det(A) =   0.281362058913013
 inv(A)=
  1.27729     -2.67948
 -.931892E-01  2.97805

 b=
 x=
  3.43630
 -2.86037

 max(error)=  2.220446049250313E-016
 Size =            3
 A=
 0.782167     0.695742     0.440628
 0.229894     0.439319     0.472592
 0.404291     0.615954     0.812874

 det(A) =   3.868436784268883E-002
 inv(A)=
  1.70655     -7.60369      3.49561
 0.108292      11.8307     -6.93686
 -.930825     -5.18290      4.74802

 b=
 x=
 -7.43101
  10.5217
 -4.62486

 max(error)=  3.330669073875470E-016
 Size =            4
 A=
 0.617587     0.715303     0.951393     0.570026
 0.201500     0.369288     0.805236E-01 0.722515
 0.594583     0.579624     0.292902     0.664950
 0.559947     0.366043E-01 0.914511     0.373112

 det(A) =  -0.108885730704336
 inv(A)=
 -1.61248     -2.23499      3.36355     0.797015
  1.70799     -.473662     0.294076E-01 -1.74459
  1.14364     0.492039     -1.75144     0.421352
 -.550738      2.19463     -.757887     0.622449

 b=
 x=
  1.43143
  2.01398
 -.357357
 -2.56754

 max(error)=  1.332267629550188E-015