1 ![ {Systems with Five--Diagonal Matrices}*)
\r
2 SUBROUTINE FDIAG (N, DL2, DL1, DM, DU1, DU2, RS, X, MARK)
\r
4 !*****************************************************************
\r
6 ! Solving a system of linear equations *
\r
8 ! with a five-diagonal, strongly nonsingular matrix A via *
\r
9 ! Gauss algorithm without pivoting. *
\r
10 ! The matrix A is given as five N-vectors DL2, DL1, DM, DU1 *
\r
11 ! and DU2. The linear system has the form: *
\r
13 ! DM(1)*X(1)+DU1(1)*X(2)+DU2(1)*X(3) = RS(1) *
\r
14 ! DL1(2)*X(1)+DM(2)*X(2)+DU1(2)*X(3)+DU2(2)*X(4) = RS(2) *
\r
16 ! DL2(I)*X(I-2)+DL1(I)*X(I-1)+ *
\r
17 ! +DM(I)*X(I)+DU1(I)*X(I+1)+DU2(I)*X(I+2) = RS(I) *
\r
18 ! for I = 3, ..., N - 2, and *
\r
20 ! DL2(N-1)*X(N-3)+DL1(N-1)*X(N-2)+ *
\r
21 ! +DM(N-1)*X(N-1)+DU1(N-1)+X(N) = RS(N-1) *
\r
22 ! DL2(N)*X(N-2)+DL1(N)*X(N-1)+DM(N)*X(N) = RS(N) *
\r
26 ! INPUT PARAMETERS: *
\r
27 ! ================= *
\r
28 ! N : number of equations; N > 3 *
\r
29 ! DL2 : N-vector DL2(1:N); second lower co-diagonal *
\r
30 ! DL2(3), DL2(4), ... , DL2(N) *
\r
31 ! DL1 : N-vector DL1(1:N); lower co-diagonal *
\r
32 ! DL1(2), DL1(3), ... , DL1(N) *
\r
33 ! DM : N-vector DM(1:N); main diagonal *
\r
34 ! DM(1), DM(2), ... , DM(N) *
\r
35 ! DU1 : N-vector DU1(1:N); upper co-diagonal *
\r
36 ! DU1(1), DU1(2), ... , DU1(N-1) *
\r
37 ! DU2 : N-vector DU2(1:N); second upper co-diagonal *
\r
38 ! DU2(1), DU2(2), ... , DU2(N-2) *
\r
39 ! RS : N-vector RS(1:N); the right hand side of the *
\r
43 ! OUTPUT PARAMETERS: *
\r
44 ! ================== *
\r
45 ! DL2 :) overwritten with auxiliary vectors defining the *
\r
46 ! DL1 :) factorization of the cyclically tridiagonal *
\r
50 ! X : N-vector X(1:N); containing the solution of the *
\r
51 ! the system of equations *
\r
52 ! MARK : error parameter *
\r
53 ! MARK=-1 : condition N > 3 is not satisfied *
\r
54 ! MARK= 0 : numerically the matrix A is not strongly *
\r
56 ! MARK= 1 : everything is o.k. *
\r
58 ! NOTE: if MARK = 1, the determinant of A is given by: *
\r
59 ! DET A = DM(1) * DM(2) * ... * DM(N) *
\r
61 !----------------------------------------------------------------*
\r
63 ! subroutines required: FDIAGP, FDIAGS, MACHPD *
\r
65 !*****************************************************************
\r
67 ! author : Gisela Engeln-Muellges *
\r
68 ! date : 05.06.1988 *
\r
69 ! source : FORTRAN 77 *
\r
71 !*****************************************************************
\r
73 IMPLICIT DOUBLEPRECISION (A - H, O - Z)
\r
74 DOUBLEPRECISION DL1 (1:N), DL2 (1:N), DM (1:N)
\r
75 DOUBLEPRECISION DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N)
\r
79 ! Factor the matrix A
\r
81 CALL FDIAGP (N, DL2, DL1, DM, DU1, DU2, MARK)
\r
83 ! if MARK = 1, update and bachsubstitute
\r
85 IF (MARK.EQ.1) THEN
\r
86 CALL FDIAGS (N, DL2, DL1, DM, DU1, DU2, RS, X)
\r
89 END SUBROUTINE FDIAG
\r
92 SUBROUTINE FDIAGP (N, DL2, DL1, DM, DU1, DU2, MARK)
\r
94 !*****************************************************************
\r
96 ! Factor a five-diagonal, strongly nonsingular matrix A *
\r
97 ! that is defined by the five N-vectors DL2, DL1, DM, DU1 *
\r
98 ! and DU2, into its triangular factors L * R by applying *
\r
99 ! Gaussian elimination specialized for five-diagonal matrices*
\r
100 ! (without pivoting). *
\r
103 ! INPUT PARAMETERS: *
\r
104 ! ================= *
\r
105 ! N : number of equations; N > 3 *
\r
106 ! DL2 : N-vector DL2(1:N); second lower co-diagonal *
\r
107 ! DL2(3), DL2(4), ... , DL2(N) *
\r
108 ! DL1 : N-vector DL1(1:N); lower co-diagonal *
\r
109 ! DL1(2), DL1(3), ... , DL1(N) *
\r
110 ! DM : N-vector DM(1:N); main diagonal *
\r
111 ! DM(1), DM(2), ... , DM(N) *
\r
112 ! DU1 : N-vector DU1(1:N); upper co-diagonal *
\r
113 ! DU1(1), DU1(2), ... , DU1(N-1) *
\r
114 ! DU2 : N-vector DU2(1:N); second upper co-diagonal *
\r
115 ! DU2(1), DU2(2), ... , DU2(N-2) *
\r
118 ! OUTPUT PARAMETERS: *
\r
119 ! ================== *
\r
120 ! DL2 :) overwritten with auxiliary vectors that define *
\r
121 ! DL1 :) the factors of the five-diagonal matrix A; *
\r
122 ! DM :) the three co-diagonals of the lower triangular *
\r
123 ! DU1 :) matrix L are stored in the vectors DL2, DL1 and *
\r
124 ! DU2 :) DM. The two co-diagonals of the unit upper *
\r
125 ! triangular matrix R are stored in the vectors DU1 *
\r
126 ! and DU2, its diagonal elements each have the *
\r
128 ! MARK : error parameter *
\r
129 ! MARK=-1 : condition N > 3 is violated *
\r
130 ! MARK= 0 : numerically the matrix is not strongly *
\r
132 ! MARK= 1 : everything is o.k. *
\r
134 !----------------------------------------------------------------*
\r
136 ! subroutines required: MACHPD *
\r
138 !*****************************************************************
\r
140 ! author : Gisela Engeln-Muellges *
\r
141 ! date : 05.06.1988 *
\r
142 ! source : FORTRAN 77 *
\r
144 !*****************************************************************
\r
146 IMPLICIT DOUBLEPRECISION (A - H, O - Z)
\r
147 DOUBLEPRECISION DL2 (1:N), DL1 (1:N), DM (1:N), DU1 (1:N), &
\r
150 ! testing whether N > 3
\r
153 IF (N.LT.4) RETURN
\r
155 ! calculating the machine constant
\r
158 10 FMACHP = 0.5D0 * FMACHP
\r
159 IF (MACHPD (1.0D0 + FMACHP) .EQ.1) GOTO 10
\r
160 FMACHP = FMACHP * 2.0D0
\r
162 ! determining relative error bounds
\r
164 EPS = 4.0D0 * FMACHP
\r
166 ! initializing the undefined vector components
\r
172 DU2 (N - 1) = 0.0D0
\r
175 ! factoring the matrix A while checking for strong nonsingularity
\r
178 ROW = DABS (DM (1) ) + DABS (DU1 (1) ) + DABS (DU2 (1) )
\r
179 IF (ROW.EQ.0.0D0) THEN
\r
184 IF (DABS (DM (1) ) * D.LE.EPS) THEN
\r
188 DU1 (1) = DU1 (1) / DM (1)
\r
189 DU2 (1) = DU2 (1) / DM (1)
\r
190 ROW = DABS (DL1 (2) ) + DABS (DM (2) ) + DABS (DU1 (2) ) + DABS ( &
\r
192 IF (ROW.EQ.0.0D0) THEN
\r
197 DM (2) = DM (2) - DL1 (2) * DU1 (1)
\r
198 IF (DABS (DM (2) ) * D.LE.EPS) THEN
\r
202 DU1 (2) = (DU1 (2) - DL1 (2) * DU2 (1) ) / DM (2)
\r
203 DU2 (2) = DU2 (2) / DM (2)
\r
205 ! factoring A while checking for strong nonsingularity of A
\r
208 ROW = DABS (DL2 (I) ) + DABS (DL1 (I) ) + DABS (DM (I) ) &
\r
209 + DABS (DU1 (I) ) + DABS (DU2 (I) )
\r
210 IF (ROW.EQ.0.0D0) THEN
\r
215 DL1 (I) = DL1 (I) - DL2 (I) * DU1 (I - 2)
\r
216 DM (I) = DM (I) - DL2 (I) * DU2 (I - 2) - DL1 (I) * DU1 (I - 1)
\r
217 IF (DABS (DM (I) ) * D.LE.EPS) THEN
\r
222 DU1 (I) = (DU1 (I) - DL1 (I) * DU2 (I - 1) ) / DM (I)
\r
224 IF (I.LT. (N - 1) ) THEN
\r
225 DU2 (I) = DU2 (I) / DM (I)
\r
230 END SUBROUTINE FDIAGP
\r
233 SUBROUTINE FDIAGS (N, DL2, DL1, DM, DU1, DU2, RS, X)
\r
235 !*****************************************************************
\r
237 ! Solving a linear system of equations *
\r
239 ! for a five-diagonal, strongly nonsingular matrix A, once *
\r
240 ! the factor matrices L * R have been calculated by *
\r
241 ! SUBROUTINE FDIAGP. Here they are used as input arrays and *
\r
242 ! they are stored in the five N-vectors DL2, DL1, DM, DU1 *
\r
246 ! INPUT PARAMETERS: *
\r
247 ! ================= *
\r
248 ! N : number of equations; N > 3 *
\r
249 ! DL2 : N-vector DL2(1:N); ) lower triangular matrix L *
\r
250 ! DL1 : N-vector DL1(1:N); ) including the diagonal *
\r
251 ! DM : N-vector DM(1:N); ) elements *
\r
253 ! DU1 : N-vector DU1(1:N); ) unit upper triangular matrix *
\r
254 ! DU2 : N-vector DU2(1:N); ) R without its unit diagonal *
\r
256 ! RS : N-vector RS1(1:N); right side of the linear system *
\r
259 ! OUTPUT PARAMETERS: *
\r
260 ! ================== *
\r
261 ! X : N-vector X(1:N); the solution of the linear system *
\r
263 !----------------------------------------------------------------*
\r
265 ! subroutines required: none *
\r
267 !*****************************************************************
\r
269 ! author : Gisela Engeln-Muellges *
\r
270 ! date : 05.06.1988 *
\r
271 ! source : FORTRAN 77 *
\r
273 !*****************************************************************
\r
275 IMPLICIT DOUBLEPRECISION (A - H, O - Z)
\r
276 DOUBLEPRECISION DL2 (1:N), DL1 (1:N), DM (1:N)
\r
277 DOUBLEPRECISION DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N)
\r
281 RS (1) = RS (1) / DM (1)
\r
282 RS (2) = (RS (2) - DL1 (2) * RS (1) ) / DM (2)
\r
284 RS (I) = (RS (I) - DL2 (I) * RS (I - 2) - DL1 (I) * RS (I - 1) &
\r
288 ! backsubstitution
\r
291 X (N - 1) = RS (N - 1) - DU1 (N - 1) * X (N)
\r
292 DO 20 I = N - 2, 1, - 1
\r
293 X (I) = RS (I) - DU1 (I) * X (I + 1) - DU2 (I) * X (I + 2)
\r
296 END SUBROUTINE FDIAGS
\r