5 C[FE{F 4.12.1}{Systems with Five-Diagonal Matrices}
\r
6 C[ {Systems with Five-Diagonal Matrices}*)
\r
8 SUBROUTINE FDIAG (N,DL2,DL1,DM,DU1,DU2,RS,X,MARK)
\r
11 C*****************************************************************
\r
13 C Solving a system of linear equations *
\r
15 C with a five-diagonal, strongly nonsingular matrix A via *
\r
16 C Gauss algorithm without pivoting. *
\r
18 C The matrix A is given as five N-vectors DL2, DL1, DM, DU1 *
\r
19 C and DU2. The linear system has the form: *
\r
21 C DM(1)*X(1)+DU1(1)*X(2)+DU2(1)*X(3) = RS(1) *
\r
22 C DL1(2)*X(1)+DM(2)*X(2)+DU1(2)*X(3)+DU2(2)*X(4) = RS(2) *
\r
24 C DL2(I)*X(I-2)+DL1(I)*X(I-1)+ *
\r
25 C +DM(I)*X(I)+DU1(I)*X(I+1)+DU2(I)*X(I+2) = RS(I) *
\r
26 C for I = 3, ..., N - 2, and *
\r
28 C DL2(N-1)*X(N-3)+DL1(N-1)*X(N-2)+ *
\r
29 C +DM(N-1)*X(N-1)+DU1(N-1)+X(N) = RS(N-1) *
\r
30 C DL2(N)*X(N-2)+DL1(N)*X(N-1)+DM(N)*X(N) = RS(N) *
\r
34 C INPUT PARAMETERS: *
\r
35 C ================= *
\r
36 C N : number of equations; N > 3 *
\r
37 C DL2 : N-vector DL2(1:N); second lower co-diagonal *
\r
38 C DL2(3), DL2(4), ... , DL2(N) *
\r
39 C DL1 : N-vector DL1(1:N); lower co-diagonal *
\r
40 C DL1(2), DL1(3), ... , DL1(N) *
\r
41 C DM : N-vector DM(1:N); main diagonal *
\r
42 C DM(1), DM(2), ... , DM(N) *
\r
43 C DU1 : N-vector DU1(1:N); upper co-diagonal *
\r
44 C DU1(1), DU1(2), ... , DU1(N-1) *
\r
45 C DU2 : N-vector DU2(1:N); second upper co-diagonal *
\r
46 C DU2(1), DU2(2), ... , DU2(N-2) *
\r
47 C RS : N-vector RS(1:N); the right hand side of the *
\r
51 C OUTPUT PARAMETERS: *
\r
52 C ================== *
\r
53 C DL2 :) overwritten with auxiliary vectors defining the *
\r
54 C DL1 :) factorization of the cyclically tridiagonal *
\r
58 C X : N-vector X(1:N); containing the solution of the *
\r
59 C the system of equations *
\r
60 C MARK : error parameter *
\r
61 C MARK=-1 : condition N > 3 is not satisfied *
\r
62 C MARK= 0 : numerically the matrix A is not strongly *
\r
64 C MARK= 1 : everything is o.k. *
\r
66 C NOTE: if MARK = 1, the determinant of A is given by: *
\r
67 C DET A = DM(1) * DM(2) * ... * DM(N) *
\r
69 C----------------------------------------------------------------*
\r
71 C subroutines required: FDIAGP, FDIAGS, MACHPD *
\r
73 C*****************************************************************
\r
75 C author : Gisela Engeln-Muellges *
\r
76 C date : 05.06.1988 *
\r
77 C source : FORTRAN 77 *
\r
80 C*****************************************************************
\r
83 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
\r
84 DOUBLE PRECISION DL1(1:N),DL2(1:N),DM(1:N)
\r
85 DOUBLE PRECISION DU1(1:N),DU2(1:N),RS(1:N),X(1:N)
\r
87 IF (N .LT. 4) RETURN
\r
89 C Factor the matrix A
\r
91 CALL FDIAGP(N,DL2,DL1,DM,DU1,DU2,MARK)
\r
93 C if MARK = 1, update and bachsubstitute
\r
95 IF (MARK .EQ. 1) THEN
\r
96 CALL FDIAGS(N,DL2,DL1,DM,DU1,DU2,RS,X)
\r
104 SUBROUTINE FDIAGP (N,DL2,DL1,DM,DU1,DU2,MARK)
\r
107 C*****************************************************************
\r
109 C Factor a five-diagonal, strongly nonsingular matrix A *
\r
110 C that is defined by the five N-vectors DL2, DL1, DM, DU1 *
\r
111 C and DU2, into its triangular factors L * R by applying *
\r
112 C Gaussian elimination specialized for five-diagonal matrices*
\r
113 C (without pivoting). *
\r
117 C INPUT PARAMETERS: *
\r
118 C ================= *
\r
119 C N : number of equations; N > 3 *
\r
120 C DL2 : N-vector DL2(1:N); second lower co-diagonal *
\r
121 C DL2(3), DL2(4), ... , DL2(N) *
\r
122 C DL1 : N-vector DL1(1:N); lower co-diagonal *
\r
123 C DL1(2), DL1(3), ... , DL1(N) *
\r
124 C DM : N-vector DM(1:N); main diagonal *
\r
125 C DM(1), DM(2), ... , DM(N) *
\r
126 C DU1 : N-vector DU1(1:N); upper co-diagonal *
\r
127 C DU1(1), DU1(2), ... , DU1(N-1) *
\r
128 C DU2 : N-vector DU2(1:N); second upper co-diagonal *
\r
129 C DU2(1), DU2(2), ... , DU2(N-2) *
\r
132 C OUTPUT PARAMETERS: *
\r
133 C ================== *
\r
134 C DL2 :) overwritten with auxiliary vectors that define *
\r
135 C DL1 :) the factors of the five-diagonal matrix A; *
\r
136 C DM :) the three co-diagonals of the lower triangular *
\r
137 C DU1 :) matrix L are stored in the vectors DL2, DL1 and *
\r
138 C DU2 :) DM. The two co-diagonals of the unit upper *
\r
139 C triangular matrix R are stored in the vectors DU1 *
\r
140 C and DU2, its diagonal elements each have the *
\r
142 C MARK : error parameter *
\r
143 C MARK=-1 : condition N > 3 is violated *
\r
144 C MARK= 0 : numerically the matrix is not strongly *
\r
146 C MARK= 1 : everything is o.k. *
\r
148 C----------------------------------------------------------------*
\r
150 C subroutines required: MACHPD *
\r
152 C*****************************************************************
\r
154 C author : Gisela Engeln-Muellges *
\r
155 C date : 05.06.1988 *
\r
156 C source : FORTRAN 77 *
\r
159 C*****************************************************************
\r
162 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
\r
163 DOUBLE PRECISION DL2(1:N),DL1(1:N),DM(1:N),DU1(1:N),DU2(1:N)
\r
165 C testing whether N > 3
\r
168 IF (N .LT. 4) RETURN
\r
170 C calculating the machine constant
\r
173 10 FMACHP = 0.5D0 * FMACHP
\r
174 IF (MACHPD(1.0D0+FMACHP) .EQ. 1) GOTO 10
\r
175 FMACHP = FMACHP * 2.0D0
\r
177 C determining relative error bounds
\r
179 EPS = 4.0D0 * FMACHP
\r
181 C initializing the undefined vector components
\r
190 C factoring the matrix A while checking for strong nonsingularity
\r
193 ROW = DABS(DM(1)) + DABS(DU1(1)) + DABS(DU2(1))
\r
194 IF (ROW .EQ. 0.0D0) THEN
\r
199 IF (DABS(DM(1))*D .LE. EPS) THEN
\r
203 DU1(1) = DU1(1)/DM(1)
\r
204 DU2(1) = DU2(1)/DM(1)
\r
205 ROW = DABS(DL1(2)) + DABS(DM(2)) + DABS(DU1(2)) + DABS(DU2(2))
\r
206 IF (ROW .EQ. 0.0D0) THEN
\r
211 DM(2) = DM(2)-DL1(2)*DU1(1)
\r
212 IF (DABS(DM(2))*D .LE. EPS) THEN
\r
216 DU1(2) = (DU1(2)-DL1(2)*DU2(1))/DM(2)
\r
217 DU2(2) = DU2(2)/DM(2)
\r
219 C factoring A while checking for strong nonsingularity of A
\r
222 ROW = DABS(DL2(I))+DABS(DL1(I))+DABS(DM(I))+
\r
223 + DABS(DU1(I))+DABS(DU2(I))
\r
224 IF (ROW .EQ. 0.0D0) THEN
\r
229 DL1(I) = DL1(I)-DL2(I)*DU1(I-2)
\r
230 DM(I) = DM(I)-DL2(I)*DU2(I-2)-DL1(I)*DU1(I-1)
\r
231 IF (DABS(DM(I))*D .LE. EPS) THEN
\r
236 DU1(I) = (DU1(I)-DL1(I)*DU2(I-1))/DM(I)
\r
238 IF (I .LT. (N-1)) THEN
\r
239 DU2(I) = DU2(I)/DM(I)
\r
249 SUBROUTINE FDIAGS (N,DL2,DL1,DM,DU1,DU2,RS,X)
\r
252 C*****************************************************************
\r
254 C Solving a linear system of equations *
\r
256 C for a five-diagonal, strongly nonsingular matrix A, once *
\r
257 C the factor matrices L * R have been calculated by *
\r
258 C SUBROUTINE FDIAGP. *
\r
260 C Here they are used as input arrays and *
\r
261 C they are stored in the five N-vectors DL2, DL1, DM, DU1 *
\r
265 C INPUT PARAMETERS: *
\r
266 C ================= *
\r
267 C N : number of equations; N > 3 *
\r
268 C DL2 : N-vector DL2(1:N); ) lower triangular matrix L *
\r
269 C DL1 : N-vector DL1(1:N); ) including the diagonal *
\r
270 C DM : N-vector DM(1:N); ) elements *
\r
272 C DU1 : N-vector DU1(1:N); ) unit upper triangular matrix *
\r
273 C DU2 : N-vector DU2(1:N); ) R without its unit diagonal *
\r
275 C RS : N-vector RS1(1:N); right side of the linear system *
\r
278 C OUTPUT PARAMETERS: *
\r
279 C ================== *
\r
280 C X : N-vector X(1:N); the solution of the linear system *
\r
282 C----------------------------------------------------------------*
\r
284 C subroutines required: none *
\r
286 C*****************************************************************
\r
288 C author : Gisela Engeln-Muellges *
\r
289 C date : 05.06.1988 *
\r
290 C source : FORTRAN 77 *
\r
293 C*****************************************************************
\r
296 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
\r
297 DOUBLE PRECISION DL2(1:N),DL1(1:N),DM(1:N)
\r
298 DOUBLE PRECISION DU1(1:N),DU2(1:N),RS(1:N),X(1:N)
\r
303 RS(2)=(RS(2)-DL1(2)*RS(1))/DM(2)
\r
305 RS(I)=(RS(I)-DL2(I)*RS(I-2)-DL1(I)*RS(I-1))/DM(I)
\r
311 X(N-1)=RS(N-1)-DU1(N-1)*X(N)
\r
313 X(I)=RS(I)-DU1(I)*X(I+1)-DU2(I)*X(I+2)
\r