source/unres/fdiag.F90

   1 ![            {Systems with Five--Diagonal Matrices}*)                  \r
   2       SUBROUTINE FDIAG (N, DL2, DL1, DM, DU1, DU2, RS, X, MARK) \r
   3 !                                                                       \r
   4 !*****************************************************************      \r
   5 !                                                                *      \r
   6 !     Solving a system of linear equations                       *      \r
   7 !                      A * X = RS                                *      \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
  12 !                                                                *      \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
  15 !                                                                *      \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
  19 !                                                                *      \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
  23 !                                                                *      \r
  24 !                                                                *      \r
  25 !                                                                *      \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
  40 !             linear system                                      *      \r
  41 !                                                                *      \r
  42 !                                                                *      \r
  43 !     OUTPUT PARAMETERS:                                         *      \r
  44 !     ==================                                         *      \r
  45 !     DL2   :) overwritten with auxiliary vectors defining the   *      \r
  46 !     DL1   :) factorization of the cyclically tridiagonal       *      \r
  47 !     DM    :) matrix A                                          *      \r
  48 !     DU1   :)                                                   *      \r
  49 !     DU2   :)                                                   *      \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
  55 !                       nonsingular                              *      \r
  56 !             MARK= 1 : everything is o.k.                       *      \r
  57 !                                                                *      \r
  58 !     NOTE: if MARK = 1, the determinant of A is given by:       *      \r
  59 !                DET A = DM(1) * DM(2) * ... * DM(N)             *      \r
  60 !                                                                *      \r
  61 !----------------------------------------------------------------*      \r
  62 !                                                                *      \r
  63 !  subroutines required: FDIAGP, FDIAGS, MACHPD                  *      \r
  64 !                                                                *      \r
  65 !*****************************************************************      \r
  66 !                                                                *      \r
  67 !  author   : Gisela Engeln-Muellges                             *      \r
  68 !  date     : 05.06.1988                                         *      \r
  69 !  source   : FORTRAN 77                                         *      \r
  70 !                                                                *      \r
  71 !*****************************************************************      \r
  72 !                                                                       \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
  76       MARK = - 1 \r
  77       IF (N.LT.4) RETURN \r
  78 !                                                                       \r
  79 !  Factor the matrix A                                                  \r
  80 !                                                                       \r
  81       CALL FDIAGP (N, DL2, DL1, DM, DU1, DU2, MARK) \r
  82 !                                                                       \r
  83 !  if MARK = 1, update and bachsubstitute                               \r
  84 !                                                                       \r
  85       IF (MARK.EQ.1) THEN \r
  86          CALL FDIAGS (N, DL2, DL1, DM, DU1, DU2, RS, X) \r
  87       ENDIF \r
  88       RETURN \r
  89       END SUBROUTINE FDIAG                          \r
  90 !                                                                       \r
  91 !                                                                       \r
  92       SUBROUTINE FDIAGP (N, DL2, DL1, DM, DU1, DU2, MARK) \r
  93 !                                                                       \r
  94 !*****************************************************************      \r
  95 !                                                                *      \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
 101 !                                                                *      \r
 102 !                                                                *      \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
 116 !                                                                *      \r
 117 !                                                                *      \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
 127 !              value  1.                                         *      \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
 131 !                       nonsingular                              *      \r
 132 !             MARK= 1 : everything is o.k.                       *      \r
 133 !                                                                *      \r
 134 !----------------------------------------------------------------*      \r
 135 !                                                                *      \r
 136 !  subroutines required: MACHPD                                  *      \r
 137 !                                                                *      \r
 138 !*****************************************************************      \r
 139 !                                                                *      \r
 140 !  author   : Gisela Engeln-Muellges                             *      \r
 141 !  date     : 05.06.1988                                         *      \r
 142 !  source   : FORTRAN 77                                         *      \r
 143 !                                                                *      \r
 144 !*****************************************************************      \r
 145 !                                                                       \r
 146       IMPLICIT DOUBLEPRECISION (A - H, O - Z) \r
 147       DOUBLEPRECISION DL2 (1:N), DL1 (1:N), DM (1:N), DU1 (1:N),        &\r
 148       DU2 (1:N)                                                         \r
 149 !                                                                       \r
 150 !  testing whether N > 3                                                \r
 151 !                                                                       \r
 152       MARK = - 1 \r
 153       IF (N.LT.4) RETURN \r
 154 !                                                                       \r
 155 !  calculating the machine constant                                     \r
 156 !                                                                       \r
 157       FMACHP = 1.0D0 \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
 161 !                                                                       \r
 162 !  determining relative error bounds                                    \r
 163 !                                                                       \r
 164       EPS = 4.0D0 * FMACHP \r
 165 !                                                                       \r
 166 !  initializing the undefined vector components                         \r
 167 !                                                                       \r
 168       DL2 (1) = 0.0D0 \r
 169       DL2 (2) = 0.0D0 \r
 170       DL1 (1) = 0.0D0 \r
 171       DU1 (N) = 0.0D0 \r
 172       DU2 (N - 1) = 0.0D0 \r
 173       DU2 (N) = 0.0D0 \r
 174 !                                                                       \r
 175 !  factoring the matrix A while checking for strong nonsingularity      \r
 176 !  for N=1, 2                                                           \r
 177 !                                                                       \r
 178       ROW = DABS (DM (1) ) + DABS (DU1 (1) ) + DABS (DU2 (1) ) \r
 179       IF (ROW.EQ.0.0D0) THEN \r
 180          MARK = 0 \r
 181          RETURN \r
 182       ENDIF \r
 183       D = 1.0D0 / ROW \r
 184       IF (DABS (DM (1) ) * D.LE.EPS) THEN \r
 185          MARK = 0 \r
 186          RETURN \r
 187       ENDIF \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
 191       DU2 (2) )                                                         \r
 192       IF (ROW.EQ.0.0D0) THEN \r
 193          MARK = 0 \r
 194          RETURN \r
 195       ENDIF \r
 196       D = 1.0D0 / ROW \r
 197       DM (2) = DM (2) - DL1 (2) * DU1 (1) \r
 198       IF (DABS (DM (2) ) * D.LE.EPS) THEN \r
 199          MARK = 0 \r
 200          RETURN \r
 201       ENDIF \r
 202       DU1 (2) = (DU1 (2) - DL1 (2) * DU2 (1) ) / DM (2) \r
 203       DU2 (2) = DU2 (2) / DM (2) \r
 204 !                                                                       \r
 205 !  factoring A while checking for strong nonsingularity of A            \r
 206 !                                                                       \r
 207       DO 20 I = 3, N, 1 \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
 211             MARK = 0 \r
 212             RETURN \r
 213          ENDIF \r
 214          D = 1.0D0 / ROW \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
 218             MARK = 0 \r
 219             RETURN \r
 220          ENDIF \r
 221          IF (I.LT.N) THEN \r
 222             DU1 (I) = (DU1 (I) - DL1 (I) * DU2 (I - 1) ) / DM (I) \r
 223          ENDIF \r
 224          IF (I.LT. (N - 1) ) THEN \r
 225             DU2 (I) = DU2 (I) / DM (I) \r
 226          ENDIF \r
 227    20 END DO \r
 228       MARK = 1 \r
 229       RETURN \r
 230       END SUBROUTINE FDIAGP                         \r
 231 !                                                                       \r
 232 !                                                                       \r
 233       SUBROUTINE FDIAGS (N, DL2, DL1, DM, DU1, DU2, RS, X) \r
 234 !                                                                       \r
 235 !*****************************************************************      \r
 236 !                                                                *      \r
 237 !     Solving a linear system of equations                       *      \r
 238 !                A * X = RS                                      *      \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
 243 !     and DU2.                                                   *      \r
 244 !                                                                *      \r
 245 !                                                                *      \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
 252 !                                                                *      \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
 255 !                                   elements                     *      \r
 256 !     RS    : N-vector RS1(1:N); right side of the linear system *      \r
 257 !                                                                *      \r
 258 !                                                                *      \r
 259 !     OUTPUT PARAMETERS:                                         *      \r
 260 !     ==================                                         *      \r
 261 !     X     : N-vector X(1:N); the solution of the linear system *      \r
 262 !                                                                *      \r
 263 !----------------------------------------------------------------*      \r
 264 !                                                                *      \r
 265 !  subroutines required: none                                    *      \r
 266 !                                                                *      \r
 267 !*****************************************************************      \r
 268 !                                                                *      \r
 269 !  author   : Gisela Engeln-Muellges                             *      \r
 270 !  date     : 05.06.1988                                         *      \r
 271 !  source   : FORTRAN 77                                         *      \r
 272 !                                                                *      \r
 273 !*****************************************************************      \r
 274 !                                                                       \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
 278 !                                                                       \r
 279 !  updating                                                             \r
 280 !                                                                       \r
 281       RS (1) = RS (1) / DM (1) \r
 282       RS (2) = (RS (2) - DL1 (2) * RS (1) ) / DM (2) \r
 283       DO 10 I = 3, N \r
 284          RS (I) = (RS (I) - DL2 (I) * RS (I - 2) - DL1 (I) * RS (I - 1) &\r
 285          ) / DM (I)                                                     \r
 286    10 END DO \r
 287 !                                                                       \r
 288 !  backsubstitution                                                     \r
 289 !                                                                       \r
 290       X (N) = RS (N) \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
 294    20 END DO \r
 295       RETURN \r
 296       END SUBROUTINE FDIAGS                         \r