source/unres/src-HCD-5D/fdiag.f

   1 C[BA*)\r
   2 C[LE*)\r
   3 C[LE*)\r
   4 C[LE*)\r
   5 C[FE{F 4.12.1}{Systems with Five-Diagonal Matrices}\r
   6 C[            {Systems with Five-Diagonal Matrices}*)\r
   7 C[LE*)\r
   8       SUBROUTINE FDIAG (N,DL2,DL1,DM,DU1,DU2,RS,X,MARK)\r
   9 C[IX{FDIAG}*)\r
  10 C\r
  11 C*****************************************************************\r
  12 C                                                                *\r
  13 C     Solving a system of linear equations                       *\r
  14 C                      A * X = RS                                *\r
  15 C     with a five-diagonal, strongly nonsingular matrix A via    *\r
  16 C     Gauss algorithm without pivoting.                          *\r
  17 C[BE*)\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
  20 C                                                                *\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
  23 C                                                                *\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
  27 C                                                                *\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
  31 C                                                                *\r
  32 C                                                                *\r
  33 C                                                                *\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
  48 C             linear system                                      *\r
  49 C                                                                *\r
  50 C                                                                *\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
  55 C     DM    :) matrix A                                          *\r
  56 C     DU1   :)                                                   *\r
  57 C     DU2   :)                                                   *\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
  63 C                       nonsingular                              *\r
  64 C             MARK= 1 : everything is o.k.                       *\r
  65 C                                                                *\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
  68 C                                                                *\r
  69 C----------------------------------------------------------------*\r
  70 C                                                                *\r
  71 C  subroutines required: FDIAGP, FDIAGS, MACHPD                  *\r
  72 C                                                                *\r
  73 C*****************************************************************\r
  74 C                                                                *\r
  75 C  author   : Gisela Engeln-Muellges                             *\r
  76 C  date     : 05.06.1988                                         *\r
  77 C  source   : FORTRAN 77                                         *\r
  78 C                                                                *\r
  79 C[BA*)\r
  80 C*****************************************************************\r
  81 C[BE*)\r
  82 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
  86       MARK = -1\r
  87       IF (N .LT. 4) RETURN\r
  88 C\r
  89 C  Factor the matrix A\r
  90 C\r
  91       CALL FDIAGP(N,DL2,DL1,DM,DU1,DU2,MARK)\r
  92 C\r
  93 C  if MARK = 1, update and bachsubstitute\r
  94 C\r
  95       IF (MARK .EQ. 1) THEN\r
  96            CALL FDIAGS(N,DL2,DL1,DM,DU1,DU2,RS,X)\r
  97       END IF\r
  98       RETURN\r
  99       END\r
 100 C\r
 101 C\r
 102 C[BA*)\r
 103 C[LE*)\r
 104       SUBROUTINE FDIAGP (N,DL2,DL1,DM,DU1,DU2,MARK)\r
 105 C[IX{FDIAGP}*)\r
 106 C\r
 107 C*****************************************************************\r
 108 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
 114 C[BE*)\r
 115 C                                                                *\r
 116 C                                                                *\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
 130 C                                                                *\r
 131 C                                                                *\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
 141 C              value  1.                                         *\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
 145 C                       nonsingular                              *\r
 146 C             MARK= 1 : everything is o.k.                       *\r
 147 C                                                                *\r
 148 C----------------------------------------------------------------*\r
 149 C                                                                *\r
 150 C  subroutines required: MACHPD                                  *\r
 151 C                                                                *\r
 152 C*****************************************************************\r
 153 C                                                                *\r
 154 C  author   : Gisela Engeln-Muellges                             *\r
 155 C  date     : 05.06.1988                                         *\r
 156 C  source   : FORTRAN 77                                         *\r
 157 C                                                                *\r
 158 C[BA*)\r
 159 C*****************************************************************\r
 160 C[BE*)\r
 161 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
 164 C\r
 165 C  testing whether N > 3\r
 166 C\r
 167       MARK = -1\r
 168       IF (N .LT. 4) RETURN\r
 169 C\r
 170 C  calculating the machine constant\r
 171 C\r
 172       FMACHP = 1.0D0\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
 176 C\r
 177 C  determining relative error bounds\r
 178 C\r
 179       EPS = 4.0D0 * FMACHP\r
 180 C\r
 181 C  initializing the undefined vector components\r
 182 C\r
 183       DL2(1) = 0.0D0\r
 184       DL2(2) = 0.0D0\r
 185       DL1(1) = 0.0D0\r
 186       DU1(N) = 0.0D0\r
 187       DU2(N-1) = 0.0D0\r
 188       DU2(N) = 0.0D0\r
 189 C\r
 190 C  factoring the matrix A while checking for strong nonsingularity\r
 191 C  for N=1, 2\r
 192 C\r
 193       ROW = DABS(DM(1)) + DABS(DU1(1)) + DABS(DU2(1))\r
 194       IF (ROW .EQ. 0.0D0) THEN\r
 195          MARK = 0\r
 196          RETURN\r
 197       ENDIF\r
 198       D = 1.0D0/ROW\r
 199       IF (DABS(DM(1))*D .LE. EPS) THEN\r
 200          MARK = 0\r
 201          RETURN\r
 202       ENDIF\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
 207          MARK = 0\r
 208          RETURN\r
 209       ENDIF\r
 210       D = 1.0D0/ROW\r
 211       DM(2) = DM(2)-DL1(2)*DU1(1)\r
 212       IF (DABS(DM(2))*D .LE. EPS) THEN\r
 213          MARK = 0\r
 214          RETURN\r
 215       ENDIF\r
 216       DU1(2) = (DU1(2)-DL1(2)*DU2(1))/DM(2)\r
 217       DU2(2) = DU2(2)/DM(2)\r
 218 C\r
 219 C  factoring A while checking for strong nonsingularity of A\r
 220 C\r
 221       DO 20 I=3,N,1\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
 225             MARK = 0\r
 226             RETURN\r
 227          ENDIF\r
 228          D = 1.0D0/ROW\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
 232             MARK = 0\r
 233             RETURN\r
 234          ENDIF\r
 235          IF (I .LT. N) THEN\r
 236             DU1(I) = (DU1(I)-DL1(I)*DU2(I-1))/DM(I)\r
 237          ENDIF\r
 238          IF (I .LT. (N-1)) THEN\r
 239             DU2(I) = DU2(I)/DM(I)\r
 240          ENDIF\r
 241    20 CONTINUE\r
 242       MARK = 1\r
 243       RETURN\r
 244       END\r
 245 C\r
 246 C\r
 247 C[BA*)\r
 248 C[LE*)\r
 249       SUBROUTINE FDIAGS (N,DL2,DL1,DM,DU1,DU2,RS,X)\r
 250 C[IX{FDIAGS}*)\r
 251 C\r
 252 C*****************************************************************\r
 253 C                                                                *\r
 254 C     Solving a linear system of equations                       *\r
 255 C                A * X = RS                                      *\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
 259 C[BE*)\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
 262 C     and DU2.                                                   *\r
 263 C                                                                *\r
 264 C                                                                *\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
 271 C                                                                *\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
 274 C                                   elements                     *\r
 275 C     RS    : N-vector RS1(1:N); right side of the linear system *\r
 276 C                                                                *\r
 277 C                                                                *\r
 278 C     OUTPUT PARAMETERS:                                         *\r
 279 C     ==================                                         *\r
 280 C     X     : N-vector X(1:N); the solution of the linear system *\r
 281 C                                                                *\r
 282 C----------------------------------------------------------------*\r
 283 C                                                                *\r
 284 C  subroutines required: none                                    *\r
 285 C                                                                *\r
 286 C*****************************************************************\r
 287 C                                                                *\r
 288 C  author   : Gisela Engeln-Muellges                             *\r
 289 C  date     : 05.06.1988                                         *\r
 290 C  source   : FORTRAN 77                                         *\r
 291 C                                                                *\r
 292 C[BA*)\r
 293 C*****************************************************************\r
 294 C[BE*)\r
 295 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
 299 C\r
 300 C  updating\r
 301 C\r
 302       RS(1)=RS(1)/DM(1)\r
 303       RS(2)=(RS(2)-DL1(2)*RS(1))/DM(2)\r
 304       DO 10 I=3,N\r
 305          RS(I)=(RS(I)-DL2(I)*RS(I-2)-DL1(I)*RS(I-1))/DM(I)\r
 306    10 CONTINUE\r
 307 C\r
 308 C  backsubstitution\r
 309 C\r
 310       X(N)=RS(N)\r
 311       X(N-1)=RS(N-1)-DU1(N-1)*X(N)\r
 312       DO 20 I=N-2,1,-1\r
 313          X(I)=RS(I)-DU1(I)*X(I+1)-DU2(I)*X(I+2)\r
 314    20 CONTINUE\r
 315       RETURN\r
 316       END\r