1 ![ {Systems with Five--Diagonal Symmetric Matrices}
\r
2 ![ {Systems with Five--Diagonal Symmetric Matrices}*)
\r
3 SUBROUTINE FDISY (N, DM, DU1, DU2, RS, X, MARK)
\r
5 !*****************************************************************
\r
7 ! Solving a system of linear equations *
\r
9 ! for a five-diagonal, symmetric and strongly nonsingular *
\r
10 ! matrix A. The matrix A is given by the three N-vectors DM, *
\r
11 ! DU1 and DU2. The system of equations has the form : *
\r
13 ! DM(1)*X(1) + DU1(1)*X(2) + DU2(1)*X(3) = RS(1) *
\r
14 ! DU1(1)*X(1) + DM(2)*X(2) + DU1(2)*X(3) + DU2(2)*X(4) = RS(2) *
\r
16 ! DU2(I-2)*X(I-2) + DU1(I-1)*X(I-1) + DM(I)*X(I) + *
\r
17 ! + DU1(I)*X(I+1) + DU2(I)*X(I+2) = RS(I) *
\r
18 ! for I = 3, ..., N - 2, and *
\r
20 ! DU2(N-3)*X(N-2) + DU1(N-2)*X(N-1) + DM(N-1)*X(N-1) + *
\r
21 ! + DU1(N-1)*X(N) = RS(N-1)*
\r
22 ! DU2(N-2)*X(N-2) + OD(N-1)*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 ! DM : N-vector DM(1:N); main diagonal of A *
\r
30 ! DM(1), DM(2), ... , DM(N) *
\r
31 ! DU1 : N-vector DU1(1:N); co-diagonal of A *
\r
32 ! DU1(1), DU1(2), ... , DU1(N-1) *
\r
33 ! DU2 : N-vector DU2(1:N); second co-diagonal of A *
\r
34 ! DU2(1), DU2(2), ... , DU2(N-2) *
\r
35 ! RS : N-vector RS(1:N); the right hand side *
\r
38 ! OUTPUT PARAMETERS: *
\r
39 ! ================== *
\r
41 ! DU1 :) overwritten with intermediate quantities *
\r
44 ! X : N-vector X(1:N) containing the solution vector *
\r
45 ! MARK : error parameter *
\r
46 ! MARK=-2 : condition N > 3 is not satisfied *
\r
47 ! MARK=-1 : A is strongly nonsingular, but not positive *
\r
49 ! MARK= 0 : numerically the matrix A is not strongly *
\r
51 ! MARK= 1 : A is positive definite *
\r
53 ! NOTE: If MARK = +/- 1, then the determinant of A is: *
\r
54 ! DET A = DM(1) * DM(2) * ... * DM(N) *
\r
56 !----------------------------------------------------------------*
\r
58 ! subroutines required: FDISYP, FDISYS, MACHPD *
\r
60 !*****************************************************************
\r
62 ! authors : Gisela Engeln-Muellges *
\r
63 ! date : 01.07.1992 *
\r
64 ! source : FORTRAN 77 *
\r
66 !*****************************************************************
\r
68 IMPLICIT DOUBLEPRECISION (A - H, O - Z)
\r
69 DOUBLEPRECISION DM (1:N), DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N)
\r
73 ! Factorization of the matrix A
\r
75 CALL FDISYP (N, DM, DU1, DU2, MARK)
\r
77 ! if MARK = +/- 1 , update and backsubstitute
\r
79 IF (MARK.EQ.1) THEN
\r
80 CALL FDISYS (N, DM, DU1, DU2, RS, X)
\r
83 END SUBROUTINE FDISY
\r
86 SUBROUTINE FDISYP (N, DM, DU1, DU2, MARK)
\r
88 !*****************************************************************
\r
90 ! Factor a five-diagonal, symmetric and strongly nonsingular *
\r
91 ! matrix A, that is given by the three N-vectors DM, DU1 and *
\r
92 ! DU2, into its Cholesky factors A = R(TRANSP) * D * R by *
\r
93 ! applying the root-free Cholesky method for five-diagonal *
\r
94 ! matrices. The form of the linear system is identical with *
\r
95 ! the one in SUBROUTINE FDISY. *
\r
98 ! INPUT PARAMETERS: *
\r
99 ! ================= *
\r
100 ! N : number of equations, N > 3 *
\r
101 ! DM : N-vector DM(1:N); main diagonal of A *
\r
102 ! DM(1), DM(2), ... , DM(N) *
\r
103 ! DU1 : N-vector DU1(1:N); upper co-diagonal of A *
\r
104 ! DU1(1), DU1(2), ... , DU1(N-1) *
\r
105 ! DU2 : N-vector DU2(1:N); second upper co-diagonal of A *
\r
106 ! DU2(1), DU2(2), ... , DU2(N-2); *
\r
107 ! due to symmetry the lower co-diagonals do not need to *
\r
108 ! be stored separately. *
\r
111 ! OUTPUT PARAMETERS: *
\r
112 ! ================== *
\r
113 ! DM :) overwritten with auxiliary vectors containing the *
\r
114 ! DU1 :) Cholesky factors of A. The co-diagonals of the unit *
\r
115 ! DU2 :) upper tridiagonal matrix R are stored in DU1 and DU2, *
\r
116 ! the diagonal matrix D in DM. *
\r
117 ! MARK : error parameter *
\r
118 ! MARK=-2 : condition N > 3 is not satisfied *
\r
119 ! MARK=-1 : A is strongly nonsingular, but not positive *
\r
121 ! MARK= 0 : numerically the matrix is not strongly *
\r
123 ! MARK= 1 : A is positive definite *
\r
125 ! NOTE : If MARK = +/-1, then the inertia of A, i. e., the *
\r
126 ! number of positive and negative eigenvalues of A, *
\r
127 ! is the same as the number of positive and negative *
\r
128 ! numbers among the components of DM. *
\r
130 !----------------------------------------------------------------*
\r
132 ! subroutines required: MACHPD *
\r
134 !*****************************************************************
\r
136 ! authors : Gisela Engeln-Muellges *
\r
137 ! date : 01.07.1988 *
\r
138 ! source : FORTRAN 77 *
\r
140 !*****************************************************************
\r
142 IMPLICIT DOUBLEPRECISION (A - H, O - Z)
\r
143 DOUBLEPRECISION DM (1:N), DU1 (1:N), DU2 (1:N)
\r
145 ! calculating the machine constant
\r
148 10 FMACHP = 0.5D0 * FMACHP
\r
149 IF (MACHPD (1.0D0 + FMACHP) .EQ.1) GOTO 10
\r
150 FMACHP = FMACHP * 2.0D0
\r
152 ! determining the relative error bound
\r
154 EPS = 4.0D0 * FMACHP
\r
156 ! checking for N > 3
\r
159 IF (N.LT.4) RETURN
\r
162 DU2 (N - 1) = 0.0D0
\r
164 ! checking for strong nonsingularity of the matrix A for N=1
\r
166 ROW = DABS (DM (1) ) + DABS (DU1 (1) ) + DABS (DU2 (1) )
\r
167 IF (ROW.EQ.0.0D0) THEN
\r
172 IF (DM (1) .LT.0.0D0) THEN
\r
175 ELSEIF (DABS (DM (1) ) * D.LE.EPS) THEN
\r
180 ! factoring A while checking for strong nonsingularity
\r
183 DU1 (1) = DU1 (1) / DM (1)
\r
185 DU2 (1) = DU2 (1) / DM (1)
\r
186 ROW = DABS (DUMMY) + DABS (DM (2) ) + DABS (DU1 (2) ) + DABS (DU2 &
\r
188 IF (ROW.EQ.0.0D0) THEN
\r
193 DM (2) = DM (2) - DUMMY * DU1 (1)
\r
194 IF (DM (2) .LT.0.0D0) THEN
\r
197 ELSEIF (DABS (DM (2) ) .LE.EPS) THEN
\r
202 DU1 (2) = (DU1 (2) - DUMMY1 * DU1 (1) ) / DM (2)
\r
204 DU2 (2) = DU2 (2) / DM (2)
\r
206 ROW = DABS (DUMMY1) + DABS (DUMMY) + DABS (DM (I) ) + DABS ( &
\r
207 DU1 (I) ) + DABS (DU2 (I) )
\r
208 IF (ROW.EQ.0.0D0) THEN
\r
213 DM (I) = DM (I) - DM (I - 1) * DU1 (I - 1) * DU1 (I - 1) &
\r
214 - DUMMY1 * DU2 (I - 2)
\r
215 IF (DM (I) .LT.0.0D0) THEN
\r
218 ELSEIF (DABS (DM (I) ) * D.LE.EPS) THEN
\r
224 DU1 (I) = (DU1 (I) - DUMMY2 * DU1 (I - 1) ) / DM (I)
\r
227 IF (I.LT.N - 1) THEN
\r
229 DU2 (I) = DU2 (I) / DM (I)
\r
234 END SUBROUTINE FDISYP
\r
237 SUBROUTINE FDISYS (N, DM, DU1, DU2, RS, X)
\r
239 !*****************************************************************
\r
241 ! Solving a linear system of equations *
\r
243 ! for a five-diagonal, symmetric and strongly nonsingular *
\r
244 ! matrix A, once its Cholesky factors have been calculated by *
\r
245 ! SUBROUTINE FDISYP. Here the factors of A are used as input *
\r
246 ! arrays and they are stored in the three N-vectors DM, DU1 *
\r
250 ! INPUT PARAMETER: *
\r
251 ! ================ *
\r
252 ! N : number of equations, N > 3 *
\r
253 ! DM : N-vector DM(1:N); diagonal matrix D *
\r
254 ! DU1 : N-vector DM(1:N); ) co-diagonals of the upper *
\r
255 ! DU2 : N-vector DM(1:N); ) triangular matrix R *
\r
256 ! RS : N-vector DM(1:N); the right hand side *
\r
259 ! OUTPUT PARAMETER: *
\r
260 ! ================= *
\r
261 ! X : N-vector X(1:N) containing the solution vector *
\r
263 !----------------------------------------------------------------*
\r
265 ! subroutines required: none *
\r
267 !*****************************************************************
\r
269 ! author : Gisela Engeln-Muellges *
\r
270 ! date : 29.04.1988 *
\r
271 ! source : FORTRAN 77 *
\r
273 !*****************************************************************
\r
275 IMPLICIT DOUBLEPRECISION (A - H, O - Z)
\r
276 DOUBLEPRECISION DM (1:N), DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N)
\r
281 RS (1) = DUMMY1 / DM (1)
\r
282 DUMMY2 = RS (2) - DU1 (1) * DUMMY1
\r
283 RS (2) = DUMMY2 / DM (2)
\r
285 DUMMY1 = RS (I) - DU1 (I - 1) * DUMMY2 - DU2 (I - 2) * DUMMY1
\r
286 RS (I) = DUMMY1 / DM (I)
\r
292 ! backsubstitution
\r
295 X (N - 1) = RS (N - 1) - DU1 (N - 1) * X (N)
\r
296 DO 20 I = N - 2, 1, - 1
\r
297 X (I) = RS (I) - DU1 (I) * X (I + 1) - DU2 (I) * X (I + 2)
\r
300 END SUBROUTINE FDISYS
\r