Introduction and theoretical background - Earth Dynamics: Deformations and Oscillations of the Rotating Earth - page 84

SUBROUTINE SPMAT(N1,N2,N3,C,R,B,M1,M2,M3)

C
C SPMAT computes the N1 by N1 matrix C connecting
C second derivatives to function values at N1 nodal points R(I).
C B is the N3=N1-2 by N2=2*N1-2 augmented matrix used internally
C in the Gaussian elimination procedure.
C M1, M2, M3 are the respective maximum dimensions for N1, N2, N3
C specified in the main programme.
C

IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION C(M1,M1),R(M1),B(M3,M2)
N=N3+1
C Construct augmented matrix A.
DO 11 I=1,N3
DO 10 J=1,N2
B(I,J)=0.D0
10 CONTINUE
B(I,I)=2.D0*(R(I+2)-R(I))
B(I,I+N3)=1.D0
IF(I.EQ.1.OR.I.EQ.N3)GO TO 11
B(I,I-1)=R(I+1)-R(I)
B(I,I+1)=R(I+2)-R(I+1)
11 CONTINUE
C Apply conditions that third derivatives are constant
C in first two and last two intervals.
B(1,2)=R(3)-R(2)-((R(2)-R(1))**2.D0)/(R(3)-R(2))
B(N3,N3-1)=R(N3+1)-R(N3)-((R(N1)-R(N1-1))**2.D0)/(R(N1-1)-R(N1-2))
B(1,1)=B(1,1)+((R(3)-R(1))*(R(2)-R(1)))/(R(3)-R(2))
B(N3,N3)=B(N3,N3)
1 +((R(N1)-R(N1-1))*(R(N1)-R(N1-2)))/(R(N1-1)-R(N1-2))
NM1=N3-1
C Eliminate sub-diagonal.
DO 13 I=1,NM1
DIAG=B(I,I)
ELL=B(I+1,I)
DO 12 J=1,N2
B(I,J)=B(I,J)/DIAG
B(I+1,J)=B(I+1,J)-B(I,J)*ELL
12 CONTINUE
13 CONTINUE
C Eliminate super-diagonal.
DO 15 I=1,NM1
NO=N-I
DIAG=B(NO,NO)
ELL=B(NO-1,NO)
DO 14 J=1,N2
B(NO,J)=B(NO,J)/DIAG
B(NO-1,J)=B(NO-1,J)-B(NO,J)*ELL
14 CONTINUE
15 CONTINUE
C Move inverse to position formerly occupied by A.
DO 17 I=1,N3
DO 16 J=1,N3
B(I,J)=B(I,J+N3)
16 CONTINUE
17 CONTINUE
C Construct matrix B in N-2 by N augmentation space.
DO 19 I=1,N3
DO 18 J=1,N1
B(I,J+N3)=0.0D0

18 CONTINUE