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!).
COLSOL
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.