share/colnew/fortran/lsyslv.f

   1 C---------------------------------------------------------------------
   2 C                            p a r t  3
   3 C          collocation system setup routines
   4 C---------------------------------------------------------------------
   5 C
   6       SUBROUTINE LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ,
   7      1           G, W, V, RHS, DMZO, INTEGS, IPVTG, IPVTW, RNORM,
   8      2           MODE, FSUB, DFSUB, GSUB, DGSUB, GUESS )
   9 C*********************************************************************
  10 C
  11 C   purpose
  12 C         this routine controls the set up and solution of a linear
  13 C      system of collocation equations.
  14 C         the matrix  g  is cast into an almost block diagonal
  15 C      form by an appropriate ordering of the columns and solved
  16 C      using the package of de boor-weiss [5]. the matrix is composed
  17 C      of n blocks. the i-th block has the size
  18 C                  integs(1,i) * integs(2,i).
  19 C      it contains in its last rows the linearized collocation
  20 C      equations, condensed as described in [2],
  21 C      and the linearized side conditions corresponding to
  22 C      the i-th subinterval.  integs(3,i)  steps of gaussian
  23 C      elimination are applied to it to achieve a  partial plu
  24 C      decomposition.  the right hand side vector is put into  rhs
  25 C      and the solution vector is returned in  delz and deldmz.
  26 C
  27 C         lsyslv operates according to one of 5 modes:
  28 C      mode = 0 - set up the collocation matrices  v , w , g
  29 C                 and the right hand side  rhs ,  and solve.
  30 C                 (for linear problems only.)
  31 C      mode = 1 - set up the collocation matrices  v , w , g
  32 C                 and the right hand sides  rhs  and  dmzo ,
  33 C                 and solve. also set up  integs .
  34 C                 (first iteration of nonlinear problems only).
  35 C      mode = 2 - set up  rhs  only and compute its norm.
  36 C      mode = 3 - set up  v, w, g  only and solve system.
  37 C      mode = 4 - perform forward and backward substitution only
  38 C                 (do not set up the matrices nor form the rhs).
  39 C
  40 C   variables
  41 C
  42 C      ig,izeta  - pointers to g,zeta respectively
  43 C                       (necessary to keep track of blocks of g
  44 C                       during matrix manipulations)
  45 C      idmz,irhs,iv,iw - pointers to  rhs,v,w rspectively
  46 C      df    - partial derivatives of f from dfsub
  47 C      rnorm - euclidean norm of rhs
  48 C      lside - number of side conditions in current and previous blocks
  49 C      iguess = 1 when current soln is user specified via  guess
  50 C             = 0 otherwise
  51 C
  52 C*********************************************************************
  53       IMPLICIT REAL*8 (A-H,O-Z)
  54       DIMENSION  Z(1), DMZ(1), DELZ(1), DELDMZ(1), XI(1), XIOLD(1)
  55       DIMENSION  G(1), W(1), V(1),  RHS(1) , DMZO(1), DUMMY(1)
  56       DIMENSION  INTEGS(3,1), IPVTG(1), IPVTW(1)
  57       DIMENSION  ZVAL(40), F(40), DGZ(40), DMVAL(20), DF(800), AT(28)
  58 C
  59       COMMON /COLOUT/ PRECIS, IOUT, IPRINT
  60       COMMON /COLLOC/ RHO(7), COEF(49)
  61       COMMON /COLORD/ K, NCOMP, MSTAR, KD,  MMAX, M(20)
  62       COMMON /COLSID/ ZETA(40), ALEFT, ARIGHT, IZETA, IZSAVE
  63       COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ
  64       COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS
  65       COMMON /COLBAS/ B(28), ACOL(28,7), ASAVE(28,4)
  66 C
  67       EXTERNAL DFSUB, DGSUB
  68 C
  69       M1 = MODE + 1
  70       GO TO (10, 30, 30, 30, 310), M1
  71 C
  72 C...  linear problem initialization
  73 C
  74    10 DO 20 I=1,MSTAR
  75    20 ZVAL(I) = 0.D0
  76 C
  77 C...  initialization
  78 C
  79    30 IDMZ = 1
  80       IDMZO = 1
  81       IRHS = 1
  82       IG = 1
  83       IW = 1
  84       IV = 1
  85       IZETA = 1
  86       LSIDE = 0
  87       IOLD = 1
  88       NCOL = 2 * MSTAR
  89       RNORM = 0.D0
  90       IF ( MODE .GT. 1 )                            GO TO 80
  91 C
  92 C...  build integs (describing block structure of matrix)
  93 C
  94       DO 70 I = 1,N
  95            INTEGS(2,I) = NCOL
  96            IF (I .LT. N)                            GO TO 40
  97            INTEGS(3,N) = NCOL
  98            LSIDE = MSTAR
  99            GO TO 60
 100    40      INTEGS(3,I) = MSTAR
 101    50      IF( LSIDE .EQ. MSTAR )                   GO TO 60
 102            IF ( ZETA(LSIDE+1) .GE. XI(I)+PRECIS )   GO TO 60
 103            LSIDE = LSIDE + 1
 104            GO TO 50
 105    60      NROW = MSTAR + LSIDE
 106    70 INTEGS(1,I) = NROW
 107    80 CONTINUE
 108       IF ( MODE .EQ. 2 )                            GO TO 90
 109 C
 110 C...  zero the matrices to be computed
 111 C
 112       LW = KD * KD * N
 113       DO 84 L = 1, LW
 114    84   W(L) = 0.D0
 115 C
 116 C...  the do loop 290 sets up the linear system of equations.
 117 C
 118   90  CONTINUE
 119       DO 290 I=1, N
 120 C
 121 C...       construct a block of  a  and a corresponding piece of  rhs.
 122 C
 123            XII = XI(I)
 124            H = XI(I+1) - XI(I)
 125            NROW = INTEGS(1,I)
 126 C
 127 C...       go thru the ncomp collocation equations and side conditions
 128 C...       in the i-th subinterval
 129 C
 130   100      IF ( IZETA .GT. MSTAR )                  GO TO 140
 131            IF ( ZETA(IZETA) .GT. XII + PRECIS )      GO TO 140
 132 C
 133 C...       build equation for a side condition.
 134 C
 135            IF ( MODE .EQ. 0 )                       GO TO 110
 136            IF ( IGUESS .NE. 1 )                     GO TO 102
 137 C
 138 C...       case where user provided current approximation
 139 C
 140            CALL GUESS (XII, ZVAL, DMVAL)
 141            GO TO 110
 142 C
 143 C...       other nonlinear case
 144 C
 145   102      IF ( MODE .NE. 1 )                       GO TO 106
 146            CALL APPROX (IOLD, XII, ZVAL, AT, COEF, XIOLD, NOLD,
 147      1          Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 2, DUMMY, 0)
 148            GO TO 110
 149   106      CALL APPROX (I, XII, ZVAL, AT, DUMMY, XI, N, Z, DMZ,
 150      1                  K, NCOMP, MMAX, M, MSTAR, 1, DUMMY, 0)
 151   108      IF ( MODE .EQ. 3 )                       GO TO 120
 152 C
 153 C...       find  rhs  boundary value.
 154 C
 155   110      CALL GSUB (IZETA, ZVAL, GVAL)
 156            RHS(NDMZ+IZETA) = -GVAL
 157            RNORM = RNORM + GVAL**2
 158            IF ( MODE .EQ. 2 )                       GO TO 130
 159 C
 160 C...       build a row of  a  corresponding to a boundary point
 161 C
 162   120      CALL GDERIV (G(IG), NROW, IZETA, ZVAL, DGZ, 1, DGSUB)
 163   130      IZETA = IZETA + 1
 164            GO TO 100
 165 C
 166 C...       assemble collocation equations
 167 C
 168   140      DO 220 J = 1, K
 169              HRHO = H * RHO(J)
 170              XCOL = XII + HRHO
 171 C
 172 C...         this value corresponds to a collocation (interior)
 173 C...         point. build the corresponding  ncomp  equations.
 174 C
 175              IF ( MODE .EQ. 0 )                     GO TO 200
 176              IF ( IGUESS .NE. 1 )                   GO TO 160
 177 C
 178 C...         use initial approximation provided by the user.
 179 C
 180              CALL GUESS (XCOL, ZVAL, DMZO(IRHS) )
 181              GO TO 170
 182 C
 183 C...         find  rhs  values
 184 C
 185   160        IF ( MODE .NE. 1 )                     GO TO 190
 186              CALL APPROX (IOLD, XCOL, ZVAL, AT, COEF, XIOLD, NOLD,
 187      1            Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 2, DMZO(IRHS), 1)
 188 C
 189   170        CALL FSUB (XCOL, ZVAL, F)
 190              DO 180 JJ = 1, NCOMP
 191                VALUE = DMZO(IRHS) - F(JJ)
 192                RHS(IRHS) = - VALUE
 193                RNORM = RNORM + VALUE**2
 194                IRHS = IRHS + 1
 195   180        CONTINUE
 196              GO TO 210
 197 C
 198 C...         evaluate former collocation solution
 199 C
 200   190        CALL APPROX (I, XCOL, ZVAL, ACOL(1,J), COEF, XI, N,
 201      1            Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 4, DUMMY, 0)
 202              IF ( MODE .EQ. 3 )                     GO TO 210
 203 C
 204 C...         fill in  rhs  values (and accumulate its norm).
 205 C
 206              CALL FSUB (XCOL, ZVAL, F)
 207              DO 195 JJ = 1, NCOMP
 208                VALUE = DMZ(IRHS) - F(JJ)
 209                RHS(IRHS) = - VALUE
 210                RNORM = RNORM + VALUE**2
 211                IRHS = IRHS + 1
 212   195        CONTINUE
 213              GO TO 220
 214 C
 215 C...         the linear case
 216 C
 217   200        CALL FSUB (XCOL, ZVAL, RHS(IRHS))
 218              IRHS = IRHS + NCOMP
 219 C
 220 C...         fill in ncomp rows of  w and v
 221 C
 222   210        CALL VWBLOK (XCOL, HRHO, J, W(IW), V(IV), IPVTW(IDMZ),
 223      1       KD, ZVAL, DF, ACOL(1,J), DMZO(IDMZO), NCOMP, DFSUB, MSING)
 224              IF ( MSING .NE. 0 )                    RETURN
 225   220      CONTINUE
 226 C
 227 C...       build global bvp matrix  g
 228 C
 229            IF ( MODE .NE. 2 )
 230      1      CALL GBLOCK (H, G(IG), NROW, IZETA, W(IW), V(IV), KD,
 231      2                  DUMMY, DELDMZ(IDMZ), IPVTW(IDMZ), 1 )
 232            IF ( I .LT. N )                          GO TO 280
 233            IZSAVE = IZETA
 234   240      IF ( IZETA .GT. MSTAR )                  GO TO 290
 235 C
 236 C...       build equation for a side condition.
 237 C
 238            IF ( MODE .EQ. 0 )                       GO TO 250
 239            IF ( IGUESS .NE. 1 )                     GO TO 245
 240 C
 241 C...       case where user provided current approximation
 242 C
 243            CALL GUESS (ARIGHT, ZVAL, DMVAL)
 244            GO TO 250
 245 C
 246 C...       other nonlinear case
 247 C
 248   245      IF ( MODE .NE. 1 )                       GO TO 246
 249            CALL APPROX (NOLD+1, ARIGHT, ZVAL, AT, COEF, XIOLD, NOLD,
 250      1          Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 1, DUMMY, 0)
 251            GO TO 250
 252   246      CALL APPROX (N+1, ARIGHT, ZVAL, AT, COEF, XI, N,
 253      1          Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 1, DUMMY, 0)
 254   248      IF ( MODE .EQ. 3 )                       GO TO 260
 255 C
 256 C...       find  rhs  boundary value.
 257 C
 258   250      CALL GSUB (IZETA, ZVAL, GVAL)
 259            RHS(NDMZ+IZETA) = - GVAL
 260            RNORM = RNORM + GVAL**2
 261            IF ( MODE .EQ. 2 )                       GO TO 270
 262 C
 263 C...       build a row of  a  corresponding to a boundary point
 264 C
 265   260      CALL GDERIV (G(IG), NROW, IZETA+MSTAR, ZVAL, DGZ, 2, DGSUB)
 266   270      IZETA = IZETA + 1
 267            GO TO 240
 268 C
 269 C...       update counters -- i-th block completed
 270 C
 271   280      IG = IG + NROW * NCOL
 272            IV = IV + KD * MSTAR
 273            IW = IW + KD * KD
 274            IDMZ = IDMZ + KD
 275            IF ( MODE .EQ. 1 )  IDMZO = IDMZO + KD
 276   290 CONTINUE
 277 C
 278 C...       assembly process completed
 279 C
 280       IF ( MODE .EQ. 0 .OR. MODE .EQ. 3 )           GO TO 300
 281       RNORM = DSQRT(RNORM / DFLOAT(NZ+NDMZ))
 282       IF ( MODE .NE. 2 )                            GO TO 300
 283       RETURN
 284 C
 285 C...  solve the linear system.
 286 C
 287 C...  matrix decomposition
 288 C
 289   300 CALL FCBLOK (G, INTEGS, N, IPVTG, DF, MSING)
 290 C
 291 C...  check for singular matrix
 292 C
 293       MSING = - MSING
 294       IF( MSING .NE. 0 )                            RETURN
 295 C
 296 C...  perform forward and backward substitution for mode=4 only.
 297 C
 298   310 CONTINUE
 299       DO 311 L = 1, NDMZ
 300         DELDMZ(L) = RHS(L)
 301   311 CONTINUE
 302       IZ = 1
 303       IDMZ = 1
 304       IW = 1
 305       IZET = 1
 306       DO 320 I=1, N
 307          NROW = INTEGS(1,I)
 308          IZETA = NROW + 1 - MSTAR
 309          IF ( I .EQ. N ) IZETA = IZSAVE
 310   322    IF ( IZET .EQ. IZETA )                     GO TO 324
 311            DELZ(IZ-1+IZET) = RHS(NDMZ+IZET)
 312            IZET = IZET + 1
 313          GO TO 322
 314   324    H = XI(I+1) - XI(I)
 315          CALL GBLOCK (H, G(1), NROW, IZETA, W(IW), V(1), KD,
 316      1                DELZ(IZ), DELDMZ(IDMZ), IPVTW(IDMZ), 2 )
 317          IZ = IZ + MSTAR
 318          IDMZ = IDMZ + KD
 319          IW = IW + KD * KD
 320          IF ( I .LT. N )                            GO TO 320
 321   326    IF ( IZET .GT. MSTAR )                     GO TO 320
 322            DELZ(IZ-1+IZET) = RHS(NDMZ+IZET)
 323            IZET = IZET + 1
 324          GO TO 326
 325   320 CONTINUE
 326 C
 327 C...  perform forward and backward substitution for mode=0,2, or 3.
 328 C
 329       CALL SBBLOK (G, INTEGS, N, IPVTG, DELZ)
 330 C
 331 C...  finally find deldmz
 332 C
 333       CALL DMZSOL (KD, MSTAR, N, V, DELZ, DELDMZ)
 334 C
 335       IF ( MODE .NE. 1 )                            RETURN
 336       DO 321 L = 1, NDMZ
 337         DMZ(L) = DMZO(L)
 338   321 CONTINUE
 339       IZ = 1
 340       IDMZ = 1
 341       IW = 1
 342       IZET = 1
 343       DO 350 I=1, N
 344          NROW = INTEGS(1,I)
 345          IZETA = NROW + 1 - MSTAR
 346          IF ( I .EQ. N ) IZETA = IZSAVE
 347   330    IF ( IZET .EQ. IZETA )                     GO TO 340
 348            Z(IZ-1+IZET) = DGZ(IZET)
 349            IZET = IZET + 1
 350          GO TO 330
 351   340    H = XI(I+1) - XI(I)
 352          CALL GBLOCK (H, G(1), NROW, IZETA, W(IW), DF, KD,
 353      1                Z(IZ), DMZ(IDMZ), IPVTW(IDMZ), 2 )
 354          IZ = IZ + MSTAR
 355          IDMZ = IDMZ + KD
 356          IW = IW + KD * KD
 357          IF ( I .LT. N )                            GO TO 350
 358   342    IF ( IZET .GT. MSTAR )                     GO TO 350
 359             Z(IZ-1+IZET) = DGZ(IZET)
 360             IZET = IZET + 1
 361          GO TO 342
 362   350 CONTINUE
 363       CALL SBBLOK (G, INTEGS, N, IPVTG, Z)
 364 C
 365 C...  finally find dmz
 366 C
 367       CALL DMZSOL (KD, MSTAR, N, V, Z, DMZ)
 368 C
 369       RETURN
 370       END