6 C[ {Systems with Five-Diagonal Symmetric Matrices}
\r
7 C[ {Systems with Five-Diagonal Symmetric Matrices}*)
\r
9 SUBROUTINE FDISY (N,DM,DU1,DU2,RS,X,MARK)
\r
12 C*****************************************************************
\r
14 C Solving a system of linear equations *
\r
16 C for a five-diagonal, symmetric and strongly nonsingular *
\r
19 C The matrix A is given by the three N-vectors DM, *
\r
20 C DU1 and DU2. The system of equations has the form : *
\r
22 C DM(1)*X(1) + DU1(1)*X(2) + DU2(1)*X(3) = RS(1) *
\r
23 C DU1(1)*X(1) + DM(2)*X(2) + DU1(2)*X(3) + DU2(2)*X(4) = RS(2) *
\r
25 C DU2(I-2)*X(I-2) + DU1(I-1)*X(I-1) + DM(I)*X(I) + *
\r
26 C + DU1(I)*X(I+1) + DU2(I)*X(I+2) = RS(I) *
\r
27 C for I = 3, ..., N - 2, and *
\r
29 C DU2(N-3)*X(N-2) + DU1(N-2)*X(N-1) + DM(N-1)*X(N-1) + *
\r
30 C + DU1(N-1)*X(N) = RS(N-1)*
\r
31 C DU2(N-2)*X(N-2) + OD(N-1)*X(N-1) + DM(N)*X(N) = RS(N) *
\r
35 C INPUT PARAMETERS: *
\r
36 C ================= *
\r
37 C N : number of equations, N > 3 *
\r
38 C DM : N-vector DM(1:N); main diagonal of A *
\r
39 C DM(1), DM(2), ... , DM(N) *
\r
40 C DU1 : N-vector DU1(1:N); co-diagonal of A *
\r
41 C DU1(1), DU1(2), ... , DU1(N-1) *
\r
42 C DU2 : N-vector DU2(1:N); second co-diagonal of A *
\r
43 C DU2(1), DU2(2), ... , DU2(N-2) *
\r
44 C RS : N-vector RS(1:N); the right hand side *
\r
47 C OUTPUT PARAMETERS: *
\r
48 C ================== *
\r
50 C DU1 :) overwritten with intermediate quantities *
\r
53 C X : N-vector X(1:N) containing the solution vector *
\r
54 C MARK : error parameter *
\r
55 C MARK=-2 : condition N > 3 is not satisfied *
\r
56 C MARK=-1 : A is strongly nonsingular, but not positive *
\r
58 C MARK= 0 : numerically the matrix A is not strongly *
\r
60 C MARK= 1 : A is positive definite *
\r
62 C NOTE: If MARK = +/- 1, then the determinant of A is: *
\r
63 C DET A = DM(1) * DM(2) * ... * DM(N) *
\r
65 C----------------------------------------------------------------*
\r
67 C subroutines required: FDISYP, FDISYS, MACHPD *
\r
69 C*****************************************************************
\r
71 C authors : Gisela Engeln-Muellges *
\r
72 C date : 01.07.1992 *
\r
73 C source : FORTRAN 77 *
\r
76 C*****************************************************************
\r
79 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
\r
80 DOUBLE PRECISION DM(1:N),DU1(1:N),DU2(1:N),RS(1:N),X(1:N)
\r
82 IF (N .LT. 4) RETURN
\r
84 C Factorization of the matrix A
\r
86 CALL FDISYP (N,DM,DU1,DU2,MARK)
\r
88 C if MARK = +/- 1 , update and backsubstitute
\r
90 IF (MARK .EQ. 1) THEN
\r
91 CALL FDISYS (N,DM,DU1,DU2,RS,X)
\r
99 SUBROUTINE FDISYP (N,DM,DU1,DU2,MARK)
\r
102 C*****************************************************************
\r
104 C Factor a five-diagonal, symmetric and strongly nonsingular *
\r
105 C matrix A, that is given by the three N-vectors DM, DU1 and *
\r
106 C DU2, into its Cholesky factors A = R(TRANSP) * D * R by *
\r
107 C applying the root-free Cholesky method for five-diagonal *
\r
108 C matrices. The form of the linear system is identical with *
\r
109 C the one in SUBROUTINE FDISY. *
\r
113 C INPUT PARAMETERS: *
\r
114 C ================= *
\r
115 C N : number of equations, N > 3 *
\r
116 C DM : N-vector DM(1:N); main diagonal of A *
\r
117 C DM(1), DM(2), ... , DM(N) *
\r
118 C DU1 : N-vector DU1(1:N); upper co-diagonal of A *
\r
119 C DU1(1), DU1(2), ... , DU1(N-1) *
\r
120 C DU2 : N-vector DU2(1:N); second upper co-diagonal of A *
\r
121 C DU2(1), DU2(2), ... , DU2(N-2); *
\r
122 C due to symmetry the lower co-diagonals do not need to *
\r
123 C be stored separately. *
\r
126 C OUTPUT PARAMETERS: *
\r
127 C ================== *
\r
128 C DM :) overwritten with auxiliary vectors containing the *
\r
129 C DU1 :) Cholesky factors of A. The co-diagonals of the unit *
\r
130 C DU2 :) upper tridiagonal matrix R are stored in DU1 and DU2, *
\r
131 C the diagonal matrix D in DM. *
\r
132 C MARK : error parameter *
\r
133 C MARK=-2 : condition N > 3 is not satisfied *
\r
134 C MARK=-1 : A is strongly nonsingular, but not positive *
\r
136 C MARK= 0 : numerically the matrix is not strongly *
\r
138 C MARK= 1 : A is positive definite *
\r
140 C NOTE : If MARK = +/-1, then the inertia of A, i. e., the *
\r
141 C number of positive and negative eigenvalues of A, *
\r
142 C is the same as the number of positive and negative *
\r
143 C numbers among the components of DM. *
\r
145 C----------------------------------------------------------------*
\r
147 C subroutines required: MACHPD *
\r
149 C*****************************************************************
\r
151 C authors : Gisela Engeln-Muellges *
\r
152 C date : 01.07.1988 *
\r
153 C source : FORTRAN 77 *
\r
155 C*****************************************************************
\r
157 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
\r
158 DOUBLE PRECISION DM(1:N),DU1(1:N),DU2(1:N)
\r
160 C calculating the machine constant
\r
163 10 FMACHP = 0.5D0 * FMACHP
\r
164 IF (MACHPD(1.0D0+FMACHP) .EQ. 1) GOTO 10
\r
165 FMACHP = FMACHP * 2.0D0
\r
167 C determining the relative error bound
\r
169 EPS = 4.0D0 * FMACHP
\r
171 C checking for N > 3
\r
174 IF (N .LT. 4) RETURN
\r
179 C checking for strong nonsingularity of the matrix A for N=1
\r
181 ROW = DABS(DM(1)) + DABS(DU1(1)) + DABS(DU2(1))
\r
182 IF (ROW .EQ. 0.0D0) THEN
\r
187 IF (DM(1) .LT. 0.0D0) THEN
\r
190 ELSEIF (DABS(DM(1))*D .LE. EPS) THEN
\r
195 C factoring A while checking for strong nonsingularity
\r
198 DU1(1) = DU1(1)/DM(1)
\r
200 DU2(1) = DU2(1)/DM(1)
\r
201 ROW = DABS(DUMMY) + DABS(DM(2)) + DABS(DU1(2)) + DABS(DU2(2))
\r
202 IF (ROW .EQ. 0.0D0) THEN
\r
207 DM(2) = DM(2) - DUMMY*DU1(1)
\r
208 IF (DM(2) .LT. 0.0D0) THEN
\r
211 ELSEIF (DABS(DM(2)) .LE. EPS) THEN
\r
216 DU1(2) = (DU1(2)-DUMMY1*DU1(1))/DM(2)
\r
218 DU2(2) = DU2(2)/DM(2)
\r
220 ROW = DABS(DUMMY1)+DABS(DUMMY)+DABS(DM(I))+DABS(DU1(I))+
\r
222 IF (ROW .EQ. 0.0D0) THEN
\r
227 DM(I) = DM(I) - DM(I-1) * DU1(I-1) * DU1(I-1)
\r
229 IF (DM(I) .LT. 0.0D0) THEN
\r
232 ELSEIF (DABS(DM(I))*D .LE. EPS) THEN
\r
238 DU1(I) = (DU1(I)-DUMMY2*DU1(I-1))/DM(I)
\r
241 IF (I .LT. N-1) THEN
\r
243 DU2(I) = DU2(I)/DM(I)
\r
253 SUBROUTINE FDISYS (N,DM,DU1,DU2,RS,X)
\r
256 C*****************************************************************
\r
258 C Solving a linear system of equations *
\r
260 C for a five-diagonal, symmetric and strongly nonsingular *
\r
263 C Before this its Cholesky must factors have been calculated by *
\r
264 C SUBROUTINE FDISYP. Here the factors of A are used as input *
\r
265 C arrays and they are stored in the three N-vectors DM, DU1 *
\r
269 C INPUT PARAMETER: *
\r
270 C ================ *
\r
271 C N : number of equations, N > 3 *
\r
272 C DM : N-vector DM(1:N); diagonal matrix D *
\r
273 C DU1 : N-vector DM(1:N); ) co-diagonals of the upper *
\r
274 C DU2 : N-vector DM(1:N); ) triangular matrix R *
\r
275 C RS : N-vector DM(1:N); the right hand side *
\r
278 C OUTPUT PARAMETER: *
\r
279 C ================= *
\r
280 C X : N-vector X(1:N) containing the solution vector *
\r
282 C----------------------------------------------------------------*
\r
284 C subroutines required: none *
\r
286 C*****************************************************************
\r
288 C author : Gisela Engeln-Muellges *
\r
289 C date : 29.04.1988 *
\r
290 C source : FORTRAN 77 *
\r
293 C*****************************************************************
\r
296 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
\r
297 DOUBLE PRECISION DM(1:N),DU1(1:N),DU2(1:N),RS(1:N),X(1:N)
\r
302 RS(1) = DUMMY1/DM(1)
\r
303 DUMMY2 = RS(2)-DU1(1)*DUMMY1
\r
304 RS(2) = DUMMY2/DM(2)
\r
306 DUMMY1 = RS(I)-DU1(I-1)*DUMMY2-DU2(I-2)*DUMMY1
\r
307 RS(I) = DUMMY1/DM(I)
\r
316 X(N-1) = RS(N-1)-DU1(N-1)*X(N)
\r
318 X(I) = RS(I)-DU1(I)*X(I+1)-DU2(I)*X(I+2)
\r