source/unres/fdisy.f90

   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
   4 !                                                                       \r
   5 !*****************************************************************      \r
   6 !                                                                *      \r
   7 !  Solving a system of linear equations                          *      \r
   8 !                     A * X = RS                                 *      \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
  12 !                                                                *      \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
  15 !                                                                *      \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
  19 !                                                                *      \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
  23 !                                                                *      \r
  24 !                                                                *      \r
  25 !                                                                *      \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
  36 !                                                                *      \r
  37 !                                                                *      \r
  38 !  OUTPUT PARAMETERS:                                            *      \r
  39 !  ==================                                            *      \r
  40 !  DM   :)                                                       *      \r
  41 !  DU1  :) overwritten with intermediate quantities              *      \r
  42 !  DU2  :)                                                       *      \r
  43 !  RS   :)                                                       *      \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
  48 !                   definite                                     *      \r
  49 !         MARK= 0 : numerically the matrix A is not strongly     *      \r
  50 !                   nonsingular                                  *      \r
  51 !         MARK= 1 : A is positive definite                       *      \r
  52 !                                                                *      \r
  53 !  NOTE: If MARK = +/- 1, then the determinant of A is:          *      \r
  54 !           DET A = DM(1) * DM(2) * ... * DM(N)                  *      \r
  55 !                                                                *      \r
  56 !----------------------------------------------------------------*      \r
  57 !                                                                *      \r
  58 !  subroutines required: FDISYP, FDISYS, MACHPD                  *      \r
  59 !                                                                *      \r
  60 !*****************************************************************      \r
  61 !                                                                *      \r
  62 !  authors  : Gisela Engeln-Muellges                             *      \r
  63 !  date     : 01.07.1992                                         *      \r
  64 !  source   : FORTRAN 77                                         *      \r
  65 !                                                                *      \r
  66 !*****************************************************************      \r
  67 !                                                                       \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
  70       MARK = - 2 \r
  71       IF (N.LT.4) RETURN \r
  72 !                                                                       \r
  73 !  Factorization of the matrix A                                        \r
  74 !                                                                       \r
  75       CALL FDISYP (N, DM, DU1, DU2, MARK) \r
  76 !                                                                       \r
  77 !  if MARK = +/- 1 , update and backsubstitute                          \r
  78 !                                                                       \r
  79       IF (MARK.EQ.1) THEN \r
  80          CALL FDISYS (N, DM, DU1, DU2, RS, X) \r
  81       ENDIF \r
  82       RETURN \r
  83       END SUBROUTINE FDISY                          \r
  84 !                                                                       \r
  85 !                                                                       \r
  86       SUBROUTINE FDISYP (N, DM, DU1, DU2, MARK) \r
  87 !                                                                       \r
  88 !*****************************************************************      \r
  89 !                                                                *      \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
  96 !                                                                *      \r
  97 !                                                                *      \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
 109 !                                                                *      \r
 110 !                                                                *      \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
 120 !                   definite                                     *      \r
 121 !         MARK= 0 : numerically the matrix is not strongly       *      \r
 122 !                   nonsingular                                  *      \r
 123 !         MARK= 1 : A is positive definite                       *      \r
 124 !                                                                *      \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
 129 !                                                                *      \r
 130 !----------------------------------------------------------------*      \r
 131 !                                                                *      \r
 132 !  subroutines required: MACHPD                                  *      \r
 133 !                                                                *      \r
 134 !*****************************************************************      \r
 135 !                                                                *      \r
 136 !  authors  : Gisela Engeln-Muellges                             *      \r
 137 !  date     : 01.07.1988                                         *      \r
 138 !  source   : FORTRAN 77                                         *      \r
 139 !                                                                *      \r
 140 !*****************************************************************      \r
 141 !                                                                       \r
 142       IMPLICIT DOUBLEPRECISION (A - H, O - Z) \r
 143       DOUBLEPRECISION DM (1:N), DU1 (1:N), DU2 (1:N) \r
 144 !                                                                       \r
 145 !   calculating the machine constant                                    \r
 146 !                                                                       \r
 147       FMACHP = 1.0D0 \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
 151 !                                                                       \r
 152 !   determining the relative error bound                                \r
 153 !                                                                       \r
 154       EPS = 4.0D0 * FMACHP \r
 155 !                                                                       \r
 156 !   checking for N > 3                                                  \r
 157 !                                                                       \r
 158       MARK = - 2 \r
 159       IF (N.LT.4) RETURN \r
 160       DU1 (N) = 0.0D0 \r
 161       DU2 (N) = 0.0D0 \r
 162       DU2 (N - 1) = 0.0D0 \r
 163 !                                                                       \r
 164 !   checking for strong nonsingularity of the matrix A for N=1          \r
 165 !                                                                       \r
 166       ROW = DABS (DM (1) ) + DABS (DU1 (1) ) + DABS (DU2 (1) ) \r
 167       IF (ROW.EQ.0.0D0) THEN \r
 168          MARK = 0 \r
 169          RETURN \r
 170       ENDIF \r
 171       D = 1.0D0 / ROW \r
 172       IF (DM (1) .LT.0.0D0) THEN \r
 173          MARK = - 1 \r
 174          RETURN \r
 175       ELSEIF (DABS (DM (1) ) * D.LE.EPS) THEN \r
 176          MARK = 0 \r
 177          RETURN \r
 178       ENDIF \r
 179 !                                                                       \r
 180 !   factoring A while checking for strong nonsingularity                \r
 181 !                                                                       \r
 182       DUMMY = DU1 (1) \r
 183       DU1 (1) = DU1 (1) / DM (1) \r
 184       DUMMY1 = DU2 (1) \r
 185       DU2 (1) = DU2 (1) / DM (1) \r
 186       ROW = DABS (DUMMY) + DABS (DM (2) ) + DABS (DU1 (2) ) + DABS (DU2 &\r
 187       (2) )                                                             \r
 188       IF (ROW.EQ.0.0D0) THEN\r
 189          MARK = 0\r
 190          RETURN \r
 191       ENDIF \r
 192       D = 1.0D0 / ROW \r
 193       DM (2) = DM (2) - DUMMY * DU1 (1) \r
 194       IF (DM (2) .LT.0.0D0) THEN \r
 195          MARK = - 1\r
 196          RETURN \r
 197       ELSEIF (DABS (DM (2) ) .LE.EPS) THEN \r
 198          MARK = 0\r
 199          RETURN \r
 200       ENDIF \r
 201       DUMMY = DU1 (2) \r
 202       DU1 (2) = (DU1 (2) - DUMMY1 * DU1 (1) ) / DM (2) \r
 203       DUMMY2 = DU2 (2) \r
 204       DU2 (2) = DU2 (2) / DM (2) \r
 205       DO 20 I = 3, N, 1 \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
 209             MARK = 0\r
 210             RETURN \r
 211          ENDIF \r
 212          D = 1.0D0 / ROW \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
 216             MARK = - 1\r
 217             RETURN \r
 218          ELSEIF (DABS (DM (I) ) * D.LE.EPS) THEN \r
 219             MARK = 0 \r
 220             RETURN \r
 221          ENDIF \r
 222          IF (I.LT.N) THEN \r
 223             DUMMY = DU1 (I) \r
 224             DU1 (I) = (DU1 (I) - DUMMY2 * DU1 (I - 1) ) / DM (I) \r
 225             DUMMY1 = DUMMY2 \r
 226          ENDIF \r
 227          IF (I.LT.N - 1) THEN \r
 228             DUMMY2 = DU2 (I) \r
 229             DU2 (I) = DU2 (I) / DM (I) \r
 230          ENDIF \r
 231    20 END DO \r
 232       MARK = 1 \r
 233       RETURN \r
 234       END SUBROUTINE FDISYP                         \r
 235 !                                                                       \r
 236 !                                                                       \r
 237       SUBROUTINE FDISYS (N, DM, DU1, DU2, RS, X) \r
 238 !                                                                       \r
 239 !*****************************************************************      \r
 240 !                                                                *      \r
 241 !  Solving a linear system of equations                          *      \r
 242 !               A * X = RS                                       *      \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
 247 !  and DU2.                                                      *      \r
 248 !                                                                *      \r
 249 !                                                                *      \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
 257 !                                                                *      \r
 258 !                                                                *      \r
 259 !  OUTPUT PARAMETER:                                             *      \r
 260 !  =================                                             *      \r
 261 !  X    : N-vector X(1:N) containing the solution vector         *      \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     : 29.04.1988                                         *      \r
 271 !  source   : FORTRAN 77                                         *      \r
 272 !                                                                *      \r
 273 !*****************************************************************      \r
 274 !                                                                       \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
 277 !                                                                       \r
 278 !  updating                                                             \r
 279 !                                                                       \r
 280       DUMMY1 = RS (1) \r
 281       RS (1) = DUMMY1 / DM (1) \r
 282       DUMMY2 = RS (2) - DU1 (1) * DUMMY1 \r
 283       RS (2) = DUMMY2 / DM (2) \r
 284       DO 10 I = 3, N, 1 \r
 285          DUMMY1 = RS (I) - DU1 (I - 1) * DUMMY2 - DU2 (I - 2) * DUMMY1 \r
 286          RS (I) = DUMMY1 / DM (I) \r
 287          DUMMY3 = DUMMY2 \r
 288          DUMMY2 = DUMMY1 \r
 289          DUMMY1 = DUMMY3 \r
 290    10 END DO \r
 291 !                                                                       \r
 292 !  backsubstitution                                                     \r
 293 !                                                                       \r
 294       X (N) = RS (N) \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
 298    20 END DO \r
 299       RETURN \r
 300       END SUBROUTINE FDISYS                         \r