Another day, another article showing ancient Fortran code and how it can be modernized by replacing it with another language, with no mention of modern Fortran.
Note: looks like the Fortran 77 update to this routine (from 1996) can be found here (see page 1003).
My list of Fortran-related books has some related codes, and there are 16 books there with “finite element” in the title.
Bathe, Klaus-Jurgen (2014). Finite Element Procedures, 2nd ed.. Prentice Hall. Text here Related code ADINA, OpenSTAP, and STAP90
The platform that hosts this blog post seems interesting: https://www.forem.com/. Apparently, it’s the platform behind https://dev.to/. I wonder if it would be a good fit for Fortran-lang to provide a platform to up-and-coming Fortran writers to build an audience.
Call me a dinosaur if you like, but I don’t see anything in the Julia code rewrite that is any more readable than a complete refactoring of the original code into Modern Fortran would generate. Frankly, I don’t see how the Julia code is any closer to the “original pseudo-code or formulas” than the original F77 code.
Yes! We need more Fortran bloggers! 
I agree, and that highlights the importance of having posts illustrating how legacy Fortran code can be modernized.
The forem accepts comments, by the way, and a contribution showing how a modern Fortran version could look like is probably welcome.
Yes, people that know modern Fortran have to be showing a better way. Otherwise, all the general computing public will see are articles like this one.
In all fairness I think people like starting afresh better than refactoring old stuff,
that gives your work more visibility. “Fortran” feels too old, maybe the commettee should rename it “Fortranly” ore something like that (just kidding)…
I think it’s one of the reasons these newer languages/libraries started off number crunching as mostly Fortran-77 wrappers, and they’re gradually replacing them with new code instead of bugfixing/refactoring them.
Second, I think when these languages take off, they do from the universities, where everyone desperately needs to maximize productivity / minimize time-to-results. In this respect, what helps is:
This is my first take on it:
module Factorization
implicit none
integer, parameter :: dp = kind(1.d0)
contains
subroutine dense_ldlt(M,F)
integer :: M, j, i, r
real(dp) :: s, t
real(dp), dimension(M,M) :: F
do j = 2, M
do i = 1, j-1
F(i, j) = F(i, j) - dot_product(F(1:i-1, i), F(1:i-1, j))
end do
s = F(j, j)
do r = 1, j-1
t = F(r, j) / F(r, r)
s = s - t * F(r, j)
F(r, j) = t
end do
F(j, j) = s
end do
end subroutine dense_ldlt
end module Factorization
program main
use Factorization
implicit none
integer, parameter :: M = 10
real(dp), dimension(M,M) :: F
call random_number(F)
call dense_ldlt(M,F)
print *, F
end program main
My doubts are:
Do I really need to pass M as a parameter to dense_ldlt in order to declare F with the correct dimension?
Does fortran does not accept the more concise -= notation for mutations (as in s = s - t*F(r,j)?
And, harder: the original post also shows how to create a new data type to represent sparse matrices, and by defining how indexing is done on those matrices, the code can be used with slight modifications in these cases. A new type in fortran could be used for that, more or less similarly?
IMO, most of the point of the blog post was showing how similar the dense and sparse versions are. What does the sparse code with a custom sparse matrix look like?
I’m also curious. I think is not bad. It is an interesting exercise in modern Fortran. If nobody comes with a solution I will try it when I can.
This could be an opportunity to reach out for those with experience with both Fortran and finite element analysis. E.g., Cyprien Rusu’s beginner videos about Fortran (on his site, mirror on youtube) uploaded by 2018 always hinted to his web site FEA for all where he continues to blog to the present day.
The nice thing about Fortran is backward compatibility. Granted, the routine doesn’t conform to “modern” programming conventions, but it still does the job.
Here’s the original routine from Bathe’s textbook (hopefully without any transcription errors!).
It’s best to compile COLSOL separately without warning flags. With gfortran use the -std=legacy flag.
SUBROUTINE COLSOL (A,V,MAXA,NN,NWK,NNM,KKK,IOUT)
C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C. .
C. P R O G R A M .
C. TO SOLVE FINITE ELEMENT STATIC EQUILIBRIUM EQUATIONS IN .
C. CORE, USING COMPACTED STORAGE AND COLUMN REDUCTION SCHEME .
C. .
C. - - INPUT VARIABLES -- .
C. A(NWK) = STIFFNESS MATRIX STORED IN COMPACTED FORM .
C. V(NN) = RIGHT-HAND-SIDE LOAD VECTOR .
C. MAXA(NNM) = VECTOR CONTAINING ADDRESSES OF DIAGONAL .
C. ELEMENTS OF STIFFNESS MATRIX IN A .
C. NN = NUMBER OF EQUATIONS .
C. NWK = NUMBER OF ELEMENTS BELOW SKYLINE OF MATRIX .
C. NNM = NN + 1 .
C. KKK = INPUT FLAG .
C. EQ. 1 TRIANGULARIZATION OF STIFFNESS MATRIX .
C. EQ. 2 REDUCTION AND BACK-SUBSTITUTION OF LOAD VECTOR .
C. IOUT = UNIT NUMBER USED FOR OUTPUT .
C. .
C. - - OUTPUT - - .
C. A(NWK) = D AND L - FACTORS OF STIFFNESS MATRIX .
C. V(NN) = DISPLACEMENT VECTOR .
C. .
C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C. THIS PROGRAM US USED IN SINGLE PRECISION ARITHMETIC ON CRAY
C. EQUIPMENT AND DOUBLE PRECISION ARITHMETIC ON IBM MACHINES,
C. ENGINEERING WORKSTATIONS AND PCS. DEACTIVATE ABOVE LINE FOR
C. SINGLE PRECISION ARITHMETIC.
C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
DIMENSION A(NWK),V(NN),MAXA(NNM)
C
C PERFORM L*D*L(T) FACTORIZATION OF STIFFNESS MATRIX
C
IF (KKK-2) 40,150,150
40 DO 140 N=1,NN
KN=MAXA(N)
KL=KN + 1
KU=MAXA(N+1) - 1
KH=KU - KL
IF (KH) 110,90,50
50 K=N - KH
IC=0
KLT=KU
DO 80 J=1,KH
IC=IC + 1
KLT=KLT - 1
KI=MAXA(K)
ND=MAXA(K+1) - KI - 1
IF (ND) 80,80,60
60 KK=MIN0(IC,ND)
C=0.
DO 70 L=1,KK
70 C=C + A(KI+L)*A(KLT+L)
A(KLT)=A(KLT) - C
80 K=K+1
90 K=N
B=0.
DO 100 KK=KL,KU
K=K - 1
KI=MAXA(K)
C=A(KK)/A(KI)
B=B + C*A(KK)
100 A(KK)=C
A(KN)=A(KN) - B
110 IF (A(KN)) 120,120,140
120 WRITE (IOUT,2000) N,A(KN)
GO TO 800
140 CONTINUE
GO TO 900
C
C REDUCE RIGHT-HAND-SIDE LOAD VECTOR
C
150 DO 180 N=1,NN
KL=MAXA(N) + 1
KU=MAXA(N+1) - 1
IF (KU-KL) 180,160,160
160 K=N
C=0.
DO 170 KK=KL,KU
K=K - 1
170 C=C + A(KK)*V(K)
V(N)=V(N) - C
180 CONTINUE
C
C BACK-SUBSTITUTE
C
DO 200 N=1,NN
K=MAXA(N)
200 V(N)=V(N)/A(K)
IF (NN.EQ.1) GO TO 900
N=NN
DO 230 L=2,NN
KL=MAXA(N) + 1
KU=MAXA(N+1) - 1
IF (KU-KL) 230,210,210
210 K=N
DO 220 KK=KL,KU
K=K - 1
220 V(K)=V(K) - A(KK)*V(N)
230 N=N - 1
GO TO 900
C
800 STOP
900 RETURN
C
2000 FORMAT(//' STOP - STIFFNESS MATRIX NOT POSITIVE DEFINITE',//,
1 ' NONPOSITIVE PIVOT FOR EQUATION ',I8,//,
2 ' PIVOT = ',E20.12 )
END
I’ve also written a few usage examples:
program main
use skyline
implicit none
type(skyline_type(:)), allocatable :: A
call example_8p4(A)
call example_8p5(A)
call reaction_difffusion(A)
contains
! EXAMPLE 8.4, pg. 710, Finite Element Procedures, Bathe
!
! Evaluate the triangular factors D and L^T of the stiffness matrix
!
! | 5 -4 1 |
! | 6 -4 1 |
! K = | 6 -4 |
! | Sym 5 |
!
subroutine example_8p4(A)
type(skyline_type(:)), allocatable, intent(out) :: A
real(dp) :: D(4)
integer :: i
allocate(skyline_type(4) :: A)
A%nnz = 9
A%dptr = [1,2,4,7,10]
A%data = [5._dp,6._dp,-4._dp,6._dp,-4._dp,1._dp,5._dp,-4._dp,1._dp ]
! Known solution for D
D = [5._dp,14._dp/5._dp,15._dp/7._dp,5._dp/6._dp]
call factorize(A)
write(*,'(9X,A20,2X,A20)') "found", "expected"
do i = 1, A%n
write(*,'("D(",I1,",",I1,") = ",E20.12,2X,E20.12)') &
i, i, A%data(A%dptr(i)), D(i)
end do
end subroutine
! EXAMPLE 8.5, pg. 710, Finite Element Procedures, Bathe
!
! Evaluate the triangular factors D and L^T of the stiffness matrix
!
! | 2 -2 -1 |
! | 3 -2 0 |
! K = | 5 -3 0 |
! | 10 4 |
! | Sym 10 |
!
subroutine example_8p5(A)
type(skyline_type(:)), allocatable, intent(out) :: A
real(dp) :: D(5)
integer :: i
allocate(skyline_type(5) :: A)
A%nnz = 12
A%dptr = [1,2,4,6,8,13]
A%data = [2,3,-2,5,-2,10,-3,10,4,0,0,-1]
! Known solution for D
D = [2._dp,1._dp,1._dp,1._dp,0.5_dp]
call factorize(A)
write(*,'(9X,A20,2X,A20)') "found", "expected"
do i = 1, A%n
write(*,'("D(",I1,",",I1,") = ",E20.12,2X,E20.12)') &
i, i, A%data(A%dptr(i)), D(i)
end do
end subroutine
! REACTION-DIFFUSION EQUATION
!
! We solve the PDE
!
! \Delta C - (Th)^2 C = 0
!
! where \Delta is the Laplace operator, and Th is Thiele's modulus.
! The boundary conditions are given by
!
! \partial_x C = 0, at x = 0 (Neumann)
! C = 1, at x = 1 (Dirichlet)
!
! We solve the equation using a finite-difference method on a grid of
! five points. The analytical solution is given by
!
! C(x) = cosh(Th*x)/cosh(Th)
!
subroutine reaction_difffusion(A)
type(skyline_type(:)), allocatable, intent(out) :: A
integer, parameter :: n = 5
real(dp), parameter :: sigma = 1.0_dp ! sigma := (Thiele)^2
real(dp), allocatable :: b(:)
real(dp) :: h, x
integer :: i
allocate(skyline_type(n-1) :: A)
allocate(b(n))
h = 1.0_dp / real(n - 1,dp)
A%nnz = 7
A%dptr = [1,2,4,6,8]
allocate(A%data(A%nnz))
A%data = [ 1._dp + 0.5_dp*h**2*sigma, &
2._dp + h**2*sigma, -1.0_dp, &
2._dp + h**2*sigma, -1.0_dp, &
2._dp + h**2*sigma, -1.0_dp ]
call factorize(A)
b(1:n-1) = [0, 0, 0, 1]
call solve(A,b)
b(n) = 1
write(*,'(8X,A20,2X,A20)') "found", "true"
do i = 1, n
x = (i-1)*h
write(*,'("c(",I2,") = ",E20.12,2X,E20.12)') i, b(i), &
cosh(sqrt(sigma)*x)/(cosh(sqrt(sigma)))
end do
end subroutine
end program
As you can see, I’m not calling COLSOL directly, but working instead on a derived type representing a Sparse matrix in Skyline storage. The derived type and procedures are defined in the following module:
module skyline
use, intrinsic :: iso_fortran_env, only: error_unit
implicit none
private
public :: skyline_type, dp
public :: factorize, solve
integer, parameter :: dp = kind(1.0d0)
type :: skyline_type(n)
integer, len :: n
integer :: nnz = -1
integer :: dptr(n + 1)
real(dp), allocatable :: data(:)
end type
interface
subroutine colsol(a,v,maxa,nn,nwk,nnm,kkk,iout)
import dp
implicit none
real(dp) :: a(*), v(*)
integer :: maxa(nnm)
integer :: nn, nwk, nnm, kkk, iout
end subroutine
end interface
contains
!> Factorize symmetric positive-definite matrix in-place as LDL^T
subroutine factorize(A)
type(skyline_type(*)), intent(inout) :: A
integer, parameter :: kkk = 1
real(dp) :: vdummy(1)
call colsol(A%data,vdummy,A%dptr,a%n,A%nnz,a%n+1,kkk,error_unit)
end subroutine
!> Solve Ax = b for x using pre-computed LDL^T factorization.
!> The solution x overwrites the RHS vector b.
subroutine solve(A,b)
type(skyline_type(*)), intent(in) :: A
real(dp), intent(inout) :: b(:) ! size should match A%n
integer, parameter :: kkk = 2
call colsol(A%data,b,A%dptr,a%n,A%nnz,a%n+1,kkk,error_unit)
end subroutine
end module
One thing that took me a good hour to realize was that the storage format used by Bathe is not the storage format shown on Wikipedia. The columns are stored upwards, starting from the diagonal entries as shown in this snapshot from the textbook:
The programs above have been tested with the Intel Fortran Compiler Classic. Unfortunately, gfortran and ifx lag behind with support of parameterized derived types. It’s not difficult however to modify the example to use a non-parameterized derived type. Some helper routines for matrix construction would be desirable in both cases.
M as a parameter to dense_ldlt in order to declare F with the correct dimension?No, you can use an assumed-shape array as the argument.
-= notation for mutations (as in s = s - t*F(r,j)?Fortran does not have such updating operators.
I don’t understand the Julia code.
The idea here is that one can define a data type and define how indexing and mutation of the elements of that data type are defined.
For a sparse matrix, this means defining a type that contains only the non-zero values and (for simplicity) the row and column indexes of these non-zero elements, as 3 vectors.
Then, one would define what it means to do val = M(i,j) for an instance that type, and what it means to do M(i,j) = val for an instance of that type.
With those definitions, an algorithm that was designed for “regular” matrices can work for this newly defined type.
I think a library with efficient representations of sparse matrices and special matrix types like diagonal matrices, permutation matrices etc. would be a great addition to the Fortran ecosystem. Something similar to what Eigen and others does for C++.
It would probably not be as elegant or aesthetic as any Julia implementation, but I think Fortran has most language features required to make it work. For example one would need to overload arithmetic operations for the matrix types.
It would be nice to be able to overload an indexing operator as well, but this is not possible in Fortran at the moment. One could work around this limitation by using a function that returns a pointer. Example:
module matrix_mod
use iso_fortran_env, only: dp => real64
implicit none
private
public fake_sparse_matrix_t
! A sparse matrix implementation in Fortran...
type fake_sparse_matrix_t
private
! ... that's really just a dense matrix because I'm lazy
real(dp), allocatable :: data(:,:)
contains
procedure el
end type
interface fake_sparse_matrix_t
module procedure init
end interface
contains
type(fake_sparse_matrix_t) pure function init(n, m) result(this)
integer, intent(in) :: n
integer, intent(in) :: m
allocate(this%data(n, m))
end function
function el(this, i, j) result(ptr)
class(fake_sparse_matrix_t), target, intent(in) :: this
integer, intent(in) :: i
integer, intent(in) :: j
real(dp), pointer :: ptr
ptr => this%data(i, j)
end function
end module
program main
use iso_fortran_env, only: dp => real64
use matrix_mod, only: fake_sparse_matrix_t
implicit none
type(fake_sparse_matrix_t) :: m
m = fake_sparse_matrix_t(4, 4)
! Insert data
m%el(1,1) = 1.0_dp
m%el(2,1) = 3.0_dp
m%el(2,2) = 4.0_dp
! Read data
write(*,*) 'm(1,1) = ', m%el(1,1)
write(*,*) 'm(2,1) = ', m%el(2,1)
write(*,*) 'm(2,2) = ', m%el(2,2)
end program
This prints out:
m(1,1) = 1.0000000000000000
m(2,1) = 3.0000000000000000
m(2,2) = 4.0000000000000000
The downside is that by using a pointer one can’t use this in any pure procedures which is a shame.
I don’t know how the efficiency would be for a custom implementation like this would be though. One benefit with Eigen and other C++ libraries is that they are able to evaluate matrix expressions like E = A + B + C * D without creating any temporaries in between. I’m not sure if the same would be possible to achieve in Fortran.
Random insertion of elements in a sparse matrix like in my example above is definitely not the best for performance, but one could also make alternative and more performant interfaces for specific tasks like that (which is what for example Eigen does).