C[BA*) C[LE*) C[LE*) C[LE*) C[FE{F 4.12.2} C[ {Systems with Five-Diagonal Symmetric Matrices} C[ {Systems with Five-Diagonal Symmetric Matrices}*) C[LE*) SUBROUTINE FDISY (N,DM,DU1,DU2,RS,X,MARK) C[IX{FDISY}*) C C***************************************************************** C * C Solving a system of linear equations * C A * X = RS * C for a five-diagonal, symmetric and strongly nonsingular * C matrix A. * C[BE*) C The matrix A is given by the three N-vectors DM, * C DU1 and DU2. The system of equations has the form : * C * C DM(1)*X(1) + DU1(1)*X(2) + DU2(1)*X(3) = RS(1) * C DU1(1)*X(1) + DM(2)*X(2) + DU1(2)*X(3) + DU2(2)*X(4) = RS(2) * C * C DU2(I-2)*X(I-2) + DU1(I-1)*X(I-1) + DM(I)*X(I) + * C + DU1(I)*X(I+1) + DU2(I)*X(I+2) = RS(I) * C for I = 3, ..., N - 2, and * C * C DU2(N-3)*X(N-2) + DU1(N-2)*X(N-1) + DM(N-1)*X(N-1) + * C + DU1(N-1)*X(N) = RS(N-1)* C DU2(N-2)*X(N-2) + OD(N-1)*X(N-1) + DM(N)*X(N) = RS(N) * C * C * C * C INPUT PARAMETERS: * C ================= * C N : number of equations, N > 3 * C DM : N-vector DM(1:N); main diagonal of A * C DM(1), DM(2), ... , DM(N) * C DU1 : N-vector DU1(1:N); co-diagonal of A * C DU1(1), DU1(2), ... , DU1(N-1) * C DU2 : N-vector DU2(1:N); second co-diagonal of A * C DU2(1), DU2(2), ... , DU2(N-2) * C RS : N-vector RS(1:N); the right hand side * C * C * C OUTPUT PARAMETERS: * C ================== * C DM :) * C DU1 :) overwritten with intermediate quantities * C DU2 :) * C RS :) * C X : N-vector X(1:N) containing the solution vector * C MARK : error parameter * C MARK=-2 : condition N > 3 is not satisfied * C MARK=-1 : A is strongly nonsingular, but not positive * C definite * C MARK= 0 : numerically the matrix A is not strongly * C nonsingular * C MARK= 1 : A is positive definite * C * C NOTE: If MARK = +/- 1, then the determinant of A is: * C DET A = DM(1) * DM(2) * ... * DM(N) * C * C----------------------------------------------------------------* C * C subroutines required: FDISYP, FDISYS, MACHPD * C * C***************************************************************** C * C authors : Gisela Engeln-Muellges * C date : 01.07.1992 * C source : FORTRAN 77 * C * C[BA*) C***************************************************************** C[BE*) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION DM(1:N),DU1(1:N),DU2(1:N),RS(1:N),X(1:N) MARK = -2 IF (N .LT. 4) RETURN C C Factorization of the matrix A C CALL FDISYP (N,DM,DU1,DU2,MARK) C C if MARK = +/- 1 , update and backsubstitute C IF (MARK .EQ. 1) THEN CALL FDISYS (N,DM,DU1,DU2,RS,X) ENDIF RETURN END C C C[BA*) C[LE*) SUBROUTINE FDISYP (N,DM,DU1,DU2,MARK) C[IX{FDISYP}*) C C***************************************************************** C * C Factor a five-diagonal, symmetric and strongly nonsingular * C matrix A, that is given by the three N-vectors DM, DU1 and * C DU2, into its Cholesky factors A = R(TRANSP) * D * R by * C applying the root-free Cholesky method for five-diagonal * C matrices. The form of the linear system is identical with * C the one in SUBROUTINE FDISY. * C[BE*) C * C * C INPUT PARAMETERS: * C ================= * C N : number of equations, N > 3 * C DM : N-vector DM(1:N); main diagonal of A * C DM(1), DM(2), ... , DM(N) * C DU1 : N-vector DU1(1:N); upper co-diagonal of A * C DU1(1), DU1(2), ... , DU1(N-1) * C DU2 : N-vector DU2(1:N); second upper co-diagonal of A * C DU2(1), DU2(2), ... , DU2(N-2); * C due to symmetry the lower co-diagonals do not need to * C be stored separately. * C * C * C OUTPUT PARAMETERS: * C ================== * C DM :) overwritten with auxiliary vectors containing the * C DU1 :) Cholesky factors of A. The co-diagonals of the unit * C DU2 :) upper tridiagonal matrix R are stored in DU1 and DU2, * C the diagonal matrix D in DM. * C MARK : error parameter * C MARK=-2 : condition N > 3 is not satisfied * C MARK=-1 : A is strongly nonsingular, but not positive * C definite * C MARK= 0 : numerically the matrix is not strongly * C nonsingular * C MARK= 1 : A is positive definite * C * C NOTE : If MARK = +/-1, then the inertia of A, i. e., the * C number of positive and negative eigenvalues of A, * C is the same as the number of positive and negative * C numbers among the components of DM. * C * C----------------------------------------------------------------* C * C subroutines required: MACHPD * C * C***************************************************************** C * C authors : Gisela Engeln-Muellges * C date : 01.07.1988 * C source : FORTRAN 77 * C * C***************************************************************** C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION DM(1:N),DU1(1:N),DU2(1:N) C C calculating the machine constant C FMACHP = 1.0D0 10 FMACHP = 0.5D0 * FMACHP IF (MACHPD(1.0D0+FMACHP) .EQ. 1) GOTO 10 FMACHP = FMACHP * 2.0D0 C C determining the relative error bound C EPS = 4.0D0 * FMACHP C C checking for N > 3 C MARK = -2 IF (N .LT. 4) RETURN DU1(N) = 0.0D0 DU2(N) = 0.0D0 DU2(N-1) = 0.0D0 C C checking for strong nonsingularity of the matrix A for N=1 C ROW = DABS(DM(1)) + DABS(DU1(1)) + DABS(DU2(1)) IF (ROW .EQ. 0.0D0) THEN MARK = 0 RETURN ENDIF D = 1.0D0/ROW IF (DM(1) .LT. 0.0D0) THEN MARK =-1 RETURN ELSEIF (DABS(DM(1))*D .LE. EPS) THEN MARK = 0 RETURN ENDIF C C factoring A while checking for strong nonsingularity C DUMMY = DU1(1) DU1(1) = DU1(1)/DM(1) DUMMY1 = DU2(1) DU2(1) = DU2(1)/DM(1) ROW = DABS(DUMMY) + DABS(DM(2)) + DABS(DU1(2)) + DABS(DU2(2)) IF (ROW .EQ. 0.0D0) THEN MARK = 0 RETURN ENDIF D = 1.0D0/ROW DM(2) = DM(2) - DUMMY*DU1(1) IF (DM(2) .LT. 0.0D0) THEN MARK =-1 RETURN ELSEIF (DABS(DM(2)) .LE. EPS) THEN MARK = 0 RETURN ENDIF DUMMY = DU1(2) DU1(2) = (DU1(2)-DUMMY1*DU1(1))/DM(2) DUMMY2 = DU2(2) DU2(2) = DU2(2)/DM(2) DO 20 I=3,N,1 ROW = DABS(DUMMY1)+DABS(DUMMY)+DABS(DM(I))+DABS(DU1(I))+ + DABS(DU2(I)) IF (ROW .EQ. 0.0D0) THEN MARK = 0 RETURN ENDIF D = 1.0D0/ROW DM(I) = DM(I) - DM(I-1) * DU1(I-1) * DU1(I-1) + -DUMMY1*DU2(I-2) IF (DM(I) .LT. 0.0D0) THEN MARK = -1 RETURN ELSEIF (DABS(DM(I))*D .LE. EPS) THEN MARK = 0 RETURN ENDIF IF (I .LT. N) THEN DUMMY = DU1(I) DU1(I) = (DU1(I)-DUMMY2*DU1(I-1))/DM(I) DUMMY1 = DUMMY2 ENDIF IF (I .LT. N-1) THEN DUMMY2 = DU2(I) DU2(I) = DU2(I)/DM(I) ENDIF 20 CONTINUE MARK = 1 RETURN END C C C[BA*) C[LE*) SUBROUTINE FDISYS (N,DM,DU1,DU2,RS,X) C[IX{FDISYS}*) C C***************************************************************** C * C Solving a linear system of equations * C A * X = RS * C for a five-diagonal, symmetric and strongly nonsingular * C matrix A. * C[BE*) C Before this its Cholesky must factors have been calculated by * C SUBROUTINE FDISYP. Here the factors of A are used as input * C arrays and they are stored in the three N-vectors DM, DU1 * C and DU2. * C * C * C INPUT PARAMETER: * C ================ * C N : number of equations, N > 3 * C DM : N-vector DM(1:N); diagonal matrix D * C DU1 : N-vector DM(1:N); ) co-diagonals of the upper * C DU2 : N-vector DM(1:N); ) triangular matrix R * C RS : N-vector DM(1:N); the right hand side * C * C * C OUTPUT PARAMETER: * C ================= * C X : N-vector X(1:N) containing the solution vector * C * C----------------------------------------------------------------* C * C subroutines required: none * C * C***************************************************************** C * C author : Gisela Engeln-Muellges * C date : 29.04.1988 * C source : FORTRAN 77 * C * C[BA*) C***************************************************************** C[BE*) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION DM(1:N),DU1(1:N),DU2(1:N),RS(1:N),X(1:N) C C updating C DUMMY1 = RS(1) RS(1) = DUMMY1/DM(1) DUMMY2 = RS(2)-DU1(1)*DUMMY1 RS(2) = DUMMY2/DM(2) DO 10 I=3,N,1 DUMMY1 = RS(I)-DU1(I-1)*DUMMY2-DU2(I-2)*DUMMY1 RS(I) = DUMMY1/DM(I) DUMMY3 = DUMMY2 DUMMY2 = DUMMY1 DUMMY1 = DUMMY3 10 CONTINUE C C backsubstitution C X(N) = RS(N) X(N-1) = RS(N-1)-DU1(N-1)*X(N) DO 20 I=N-2,1,-1 X(I) = RS(I)-DU1(I)*X(I+1)-DU2(I)*X(I+2) 20 CONTINUE RETURN END