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

   1 C[BA*)\r
   2 C[LE*)\r
   3 C[LE*)\r
   4 C[LE*)\r
   5 C[FE{F 4.12.2}\r
   6 C[  {Systems with Five-Diagonal Symmetric Matrices}\r
   7 C[  {Systems with Five-Diagonal Symmetric Matrices}*)\r
   8 C[LE*)\r
   9       SUBROUTINE FDISY (N,DM,DU1,DU2,RS,X,MARK)\r
  10 C[IX{FDISY}*)\r
  11 C\r
  12 C*****************************************************************\r
  13 C                                                                *\r
  14 C  Solving a system of linear equations                          *\r
  15 C                     A * X = RS                                 *\r
  16 C  for a five-diagonal, symmetric and strongly nonsingular       *\r
  17 C  matrix A.                                                     *\r
  18 C[BE*)\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
  21 C                                                                *\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
  24 C                                                                *\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
  28 C                                                                *\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
  32 C                                                                *\r
  33 C                                                                *\r
  34 C                                                                *\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
  45 C                                                                *\r
  46 C                                                                *\r
  47 C  OUTPUT PARAMETERS:                                            *\r
  48 C  ==================                                            *\r
  49 C  DM   :)                                                       *\r
  50 C  DU1  :) overwritten with intermediate quantities              *\r
  51 C  DU2  :)                                                       *\r
  52 C  RS   :)                                                       *\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
  57 C                   definite                                     *\r
  58 C         MARK= 0 : numerically the matrix A is not strongly     *\r
  59 C                   nonsingular                                  *\r
  60 C         MARK= 1 : A is positive definite                       *\r
  61 C                                                                *\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
  64 C                                                                *\r
  65 C----------------------------------------------------------------*\r
  66 C                                                                *\r
  67 C  subroutines required: FDISYP, FDISYS, MACHPD                  *\r
  68 C                                                                *\r
  69 C*****************************************************************\r
  70 C                                                                *\r
  71 C  authors  : Gisela Engeln-Muellges                             *\r
  72 C  date     : 01.07.1992                                         *\r
  73 C  source   : FORTRAN 77                                         *\r
  74 C                                                                *\r
  75 C[BA*)\r
  76 C*****************************************************************\r
  77 C[BE*)\r
  78 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
  81       MARK = -2\r
  82       IF (N .LT. 4) RETURN\r
  83 C\r
  84 C  Factorization of the matrix A\r
  85 C\r
  86       CALL FDISYP (N,DM,DU1,DU2,MARK)\r
  87 C\r
  88 C  if MARK = +/- 1 , update and backsubstitute\r
  89 C\r
  90       IF (MARK .EQ. 1) THEN\r
  91          CALL FDISYS (N,DM,DU1,DU2,RS,X)\r
  92       ENDIF\r
  93       RETURN\r
  94       END\r
  95 C\r
  96 C\r
  97 C[BA*)\r
  98 C[LE*)\r
  99       SUBROUTINE FDISYP (N,DM,DU1,DU2,MARK)\r
 100 C[IX{FDISYP}*)\r
 101 C\r
 102 C*****************************************************************\r
 103 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
 110 C[BE*)\r
 111 C                                                                *\r
 112 C                                                                *\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
 124 C                                                                *\r
 125 C                                                                *\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
 135 C                   definite                                     *\r
 136 C         MARK= 0 : numerically the matrix is not strongly       *\r
 137 C                   nonsingular                                  *\r
 138 C         MARK= 1 : A is positive definite                       *\r
 139 C                                                                *\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
 144 C                                                                *\r
 145 C----------------------------------------------------------------*\r
 146 C                                                                *\r
 147 C  subroutines required: MACHPD                                  *\r
 148 C                                                                *\r
 149 C*****************************************************************\r
 150 C                                                                *\r
 151 C  authors  : Gisela Engeln-Muellges                             *\r
 152 C  date     : 01.07.1988                                         *\r
 153 C  source   : FORTRAN 77                                         *\r
 154 C                                                                *\r
 155 C*****************************************************************\r
 156 C\r
 157       IMPLICIT DOUBLE PRECISION (A-H,O-Z)\r
 158       DOUBLE PRECISION DM(1:N),DU1(1:N),DU2(1:N)\r
 159 C\r
 160 C   calculating the machine constant\r
 161 C\r
 162       FMACHP = 1.0D0\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
 166 C\r
 167 C   determining the relative error bound\r
 168 C\r
 169       EPS = 4.0D0 * FMACHP\r
 170 C\r
 171 C   checking for N > 3\r
 172 C\r
 173       MARK = -2\r
 174       IF (N .LT. 4) RETURN\r
 175       DU1(N) = 0.0D0\r
 176       DU2(N) = 0.0D0\r
 177       DU2(N-1) = 0.0D0\r
 178 C\r
 179 C   checking for strong nonsingularity of the matrix A for N=1\r
 180 C\r
 181       ROW = DABS(DM(1)) + DABS(DU1(1)) + DABS(DU2(1))\r
 182       IF (ROW .EQ. 0.0D0) THEN\r
 183          MARK = 0\r
 184          RETURN\r
 185       ENDIF\r
 186       D = 1.0D0/ROW\r
 187       IF (DM(1) .LT. 0.0D0) THEN\r
 188          MARK =-1\r
 189          RETURN\r
 190       ELSEIF (DABS(DM(1))*D .LE. EPS) THEN\r
 191          MARK = 0\r
 192          RETURN\r
 193       ENDIF\r
 194 C\r
 195 C   factoring A while checking for strong nonsingularity\r
 196 C\r
 197       DUMMY = DU1(1)\r
 198       DU1(1) = DU1(1)/DM(1)\r
 199       DUMMY1 = DU2(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
 203          MARK = 0\r
 204          RETURN\r
 205       ENDIF\r
 206       D = 1.0D0/ROW\r
 207       DM(2) = DM(2) - DUMMY*DU1(1)\r
 208       IF (DM(2) .LT. 0.0D0) THEN\r
 209          MARK =-1\r
 210          RETURN\r
 211       ELSEIF (DABS(DM(2)) .LE. EPS) THEN\r
 212          MARK = 0\r
 213          RETURN\r
 214       ENDIF\r
 215       DUMMY = DU1(2)\r
 216       DU1(2) = (DU1(2)-DUMMY1*DU1(1))/DM(2)\r
 217       DUMMY2 = DU2(2)\r
 218       DU2(2) = DU2(2)/DM(2)\r
 219       DO 20 I=3,N,1\r
 220          ROW = DABS(DUMMY1)+DABS(DUMMY)+DABS(DM(I))+DABS(DU1(I))+\r
 221      +         DABS(DU2(I))\r
 222          IF (ROW .EQ. 0.0D0) THEN\r
 223             MARK = 0\r
 224             RETURN\r
 225          ENDIF\r
 226          D = 1.0D0/ROW\r
 227          DM(I) = DM(I) - DM(I-1) * DU1(I-1) * DU1(I-1)\r
 228      +           -DUMMY1*DU2(I-2)\r
 229          IF (DM(I) .LT. 0.0D0) THEN\r
 230             MARK = -1\r
 231             RETURN\r
 232          ELSEIF (DABS(DM(I))*D .LE. EPS) THEN\r
 233             MARK = 0\r
 234             RETURN\r
 235          ENDIF\r
 236          IF (I .LT. N) THEN\r
 237             DUMMY = DU1(I)\r
 238             DU1(I) = (DU1(I)-DUMMY2*DU1(I-1))/DM(I)\r
 239             DUMMY1 = DUMMY2\r
 240          ENDIF\r
 241          IF (I .LT. N-1) THEN\r
 242             DUMMY2 = DU2(I)\r
 243             DU2(I) = DU2(I)/DM(I)\r
 244          ENDIF\r
 245    20 CONTINUE\r
 246       MARK = 1\r
 247       RETURN\r
 248       END\r
 249 C\r
 250 C\r
 251 C[BA*)\r
 252 C[LE*)\r
 253       SUBROUTINE FDISYS (N,DM,DU1,DU2,RS,X)\r
 254 C[IX{FDISYS}*)\r
 255 C\r
 256 C*****************************************************************\r
 257 C                                                                *\r
 258 C  Solving a linear system of equations                          *\r
 259 C               A * X = RS                                       *\r
 260 C  for a five-diagonal, symmetric and strongly nonsingular       *\r
 261 C  matrix A.                                                     *\r
 262 C[BE*)\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
 266 C  and DU2.                                                      *\r
 267 C                                                                *\r
 268 C                                                                *\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
 276 C                                                                *\r
 277 C                                                                *\r
 278 C  OUTPUT PARAMETER:                                             *\r
 279 C  =================                                             *\r
 280 C  X    : N-vector X(1:N) containing the solution vector         *\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     : 29.04.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 DM(1:N),DU1(1:N),DU2(1:N),RS(1:N),X(1:N)\r
 298 C\r
 299 C  updating\r
 300 C\r
 301       DUMMY1 = RS(1)\r
 302       RS(1) = DUMMY1/DM(1)\r
 303       DUMMY2 = RS(2)-DU1(1)*DUMMY1\r
 304       RS(2) = DUMMY2/DM(2)\r
 305       DO 10 I=3,N,1\r
 306          DUMMY1 = RS(I)-DU1(I-1)*DUMMY2-DU2(I-2)*DUMMY1\r
 307          RS(I) = DUMMY1/DM(I)\r
 308          DUMMY3 = DUMMY2\r
 309          DUMMY2 = DUMMY1\r
 310          DUMMY1 = DUMMY3\r
 311    10 CONTINUE\r
 312 C\r
 313 C  backsubstitution\r
 314 C\r
 315       X(N) = RS(N)\r
 316       X(N-1) = RS(N-1)-DU1(N-1)*X(N)\r
 317       DO 20 I=N-2,1,-1\r
 318          X(I) = RS(I)-DU1(I)*X(I+1)-DU2(I)*X(I+2)\r
 319    20 CONTINUE\r
 320       RETURN\r
 321       END\r