share/colnew/fortran/newmsh.f

   1 C----------------------------------------------------------------------
   2 C                            p a r t  2
   3 C          mesh selection, error estimation, (and related
   4 C          constant assignment) routines -- see [3], [4]
   5 C----------------------------------------------------------------------
   6 C
   7       SUBROUTINE NEWMSH (MODE, XI, XIOLD, Z, DMZ, VALSTR,
   8      1                   SLOPE, ACCUM, NFXPNT, FIXPNT)
   9 C
  10 C**********************************************************************
  11 C
  12 C   purpose
  13 C            select a mesh on which a collocation solution is to be
  14 C            determined
  15 C
  16 C                           there are 5 possible modes of action:
  17 C            mode = 5,4,3 - deal mainly with definition of an initial
  18 C                           mesh for the current boundary value problem
  19 C                 = 2,1   - deal with definition of a new mesh, either
  20 C                           by simple mesh halving or by mesh selection
  21 C            more specifically, for
  22 C            mode = 5  an initial (generally nonuniform) mesh is
  23 C                      defined by the user and no mesh selection is to
  24 C                      be performed
  25 C                 = 4  an initial (generally nonuniform) mesh is
  26 C                      defined by the user
  27 C                 = 3  a simple uniform mesh (except possibly for some
  28 C                      fixed points) is defined; n= no. of subintervals
  29 C                 = 1  the automatic mesh selection procedure is used
  30 C                      (see [3] for details)
  31 C                 = 2  a simple mesh halving is performed
  32 C
  33 C**********************************************************************
  34 C
  35 C   variables
  36 C
  37 C            n      = number of mesh subintervals
  38 C            nold   = number of subintervals for former mesh
  39 C            xi     - mesh point array
  40 C            xiold  - former mesh point array
  41 C            mshlmt - maximum no. of mesh selections which are permitted
  42 C                     for a given n before mesh halving
  43 C            mshnum - no. of mesh selections which have actually been
  44 C                     performed for the given n
  45 C            mshalt - no. of consecutive times ( plus 1 ) the mesh
  46 C                     selection has alternately halved and doubled n.
  47 C                     if mshalt .ge. mshlmt then  contrl  requires
  48 C                     that the current mesh be halved.
  49 C            mshflg = 1  the mesh is a halving of its former mesh
  50 C                       (so an error estimate has been calculated)
  51 C                   = 0  otherwise
  52 C            iguess - ipar(9) in subroutine colnew.  it is used
  53 C                     here only for mode=5 and 4, where
  54 C                   = 2 the subroutine sets xi=xiold.  this is
  55 C                       used e.g. if continuation is being per-
  56 C                       formed, and a mesh for the old differen-
  57 C                       tial equation is being used
  58 C                   = 3 same as for =2, except xi uses every other
  59 C                       point of xiold (so mesh xiold is mesh xi
  60 C                       halved)
  61 C                   = 4 xi has been defined by the user, and an old
  62 C                       mesh xiold is also available
  63 C                       otherwise, xi has been defined by the user
  64 C                       and we set xiold=xi in this subroutine
  65 C            slope  - an approximate quantity to be equidistributed for
  66 C                     mesh selection (see [3]), viz,
  67 C                             .                        (k+mj)
  68 C                     slope(i)=     max   (weight(l) *u      (xi(i)))
  69 C                               1.le.l.le.ntol         j
  70 C
  71 C                     where j=jtol(l)
  72 C            slphmx - maximum of slope(i)*(xiold(i+1)-xiold(i)) for
  73 C                     i = 1 ,..., nold.
  74 C            accum  - accum(i) is the integral of  slope  from  aleft
  75 C                     to  xiold(i).
  76 C            valstr - is assigned values needed in  errchk  for the
  77 C                     error estimate.
  78 C**********************************************************************
  79 C
  80       IMPLICIT REAL*8 (A-H,O-Z)
  81       DIMENSION D1(40), D2(40), SLOPE(1), ACCUM(1), VALSTR(1)
  82       DIMENSION XI(1), XIOLD(1), Z(1), DMZ(1), FIXPNT(1), DUMMY(1)
  83 C
  84       COMMON /COLOUT/ PRECIS, IOUT, IPRINT
  85       COMMON /COLORD/ K, NCOMP, MSTAR, KD, MMAX, M(20)
  86       COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ
  87       COMMON /COLMSH/ MSHFLG, MSHNUM, MSHLMT, MSHALT
  88       COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS
  89       COMMON /COLSID/  ZETA(40), ALEFT, ARIGHT, IZETA, IDUM
  90       COMMON /COLBAS/ B(28), ACOL(28,7), ASAVE(28,4)
  91       COMMON /COLEST/ TOL(40), WGTMSH(40), WGTERR(40), TOLIN(40),
  92      1                ROOT(40), JTOL(40), LTOL(40), NTOL
  93 C
  94       NFXP1 = NFXPNT +1
  95       GO TO (180, 100, 50, 20, 10), MODE
  96 C
  97 C...  mode=5   set mshlmt=1 so that no mesh selection is performed
  98 C
  99    10 MSHLMT = 1
 100 C
 101 C...  mode=4   the user-specified initial mesh is already in place.
 102 C
 103    20 IF ( IGUESS .LT. 2 )                          GO TO 40
 104 C
 105 C...  iguess=2, 3 or 4.
 106 C
 107       NOLDP1 = NOLD + 1
 108       IF (IPRINT .LT. 1)  WRITE(IOUT,360) NOLD, (XIOLD(I), I=1,NOLDP1)
 109       IF ( IGUESS .NE. 3 )                          GO TO 40
 110 C
 111 C...  if iread ( ipar(8) ) .ge. 1 and iguess ( ipar(9) ) .eq. 3
 112 C...  then the first mesh is every second point of the
 113 C...  mesh in  xiold .
 114 C
 115       N = NOLD /2
 116       I = 0
 117       DO 30 J = 1, NOLD, 2
 118            I = I + 1
 119    30 XI(I) = XIOLD(J)
 120    40 CONTINUE
 121       NP1 = N + 1
 122       XI(1) = ALEFT
 123       XI(NP1) = ARIGHT
 124       GO TO 320
 125 C
 126 C...  mode=3   generate a (piecewise) uniform mesh. if there are
 127 C...  fixed points then ensure that the n being used is large enough.
 128 C
 129    50 IF ( N .LT. NFXP1 )  N = NFXP1
 130       NP1 = N + 1
 131       XI(1) = ALEFT
 132       ILEFT = 1
 133       XLEFT = ALEFT
 134 C
 135 C...  loop over the subregions between fixed points.
 136 C
 137       DO 90 J = 1, NFXP1
 138            XRIGHT = ARIGHT
 139            IRIGHT = NP1
 140            IF ( J .EQ. NFXP1 )                      GO TO 60
 141            XRIGHT = FIXPNT(J)
 142 C
 143 C...       determine where the j-th fixed point should fall in the
 144 C...       new mesh - this is xi(iright) and the (j-1)st fixed
 145 C...       point is in xi(ileft)
 146 C
 147            NMIN = (XRIGHT-ALEFT) / (ARIGHT-ALEFT) * DFLOAT(N) + 1.5D0
 148            IF (NMIN .GT. N-NFXPNT+J)  NMIN = N - NFXPNT + J
 149            IRIGHT = MAX0 (ILEFT+1, NMIN)
 150    60      XI(IRIGHT) = XRIGHT
 151 C
 152 C...       generate equally spaced points between the j-1st and the
 153 C...       j-th fixed points.
 154 C
 155            NREGN = IRIGHT - ILEFT - 1
 156            IF ( NREGN .EQ. 0 )                      GO TO 80
 157            DX = (XRIGHT - XLEFT) / DFLOAT(NREGN+1)
 158            DO 70 I = 1, NREGN
 159    70      XI(ILEFT+I) = XLEFT  +  DFLOAT(I) * DX
 160    80      ILEFT = IRIGHT
 161            XLEFT = XRIGHT
 162    90 CONTINUE
 163       GO TO 320
 164 C
 165 C...  mode=2  halve the current mesh (i.e. double its size)
 166 C
 167   100 N2 = 2 * N
 168 C
 169 C...  check that n does not exceed storage limitations
 170 C
 171       IF ( N2 .LE. NMAX )                           GO TO 120
 172 C
 173 C...  if possible, try with n=nmax. redistribute first.
 174 C
 175       IF ( MODE .EQ. 2 )                            GO TO 110
 176       N = NMAX / 2
 177       GO TO 220
 178   110 IF ( IPRINT .LT. 1 )  WRITE(IOUT,370)
 179       N = N2
 180       RETURN
 181 C
 182 C...  calculate the old approximate solution values at
 183 C...  points to be used in  errchk  for error estimates.
 184 C...  if  mshflg  =1 an error estimate was obtained for
 185 C...  for the old approximation so half the needed values
 186 C...  will already be in  valstr .
 187 C
 188   120 IF ( MSHFLG .EQ. 0 )                          GO TO 140
 189 C
 190 C...  save in  valstr  the values of the old solution
 191 C...  at the relative positions 1/6 and 5/6 in each subinterval.
 192 C
 193       KSTORE = 1
 194       DO 130 I = 1, NOLD
 195           HD6 = (XIOLD(I+1) - XIOLD(I)) / 6.D0
 196           X = XIOLD(I) + HD6
 197           CALL APPROX (I, X, VALSTR(KSTORE), ASAVE(1,1), DUMMY, XIOLD,
 198      1         NOLD, Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 4, DUMMY, 0)
 199           X = X + 4.D0 * HD6
 200           KSTORE = KSTORE  +  3 * MSTAR
 201           CALL APPROX (I, X, VALSTR(KSTORE), ASAVE(1,4), DUMMY, XIOLD,
 202      1         NOLD, Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 4, DUMMY, 0)
 203           KSTORE = KSTORE  +  MSTAR
 204   130 CONTINUE
 205       GO TO 160
 206 C
 207 C...  save in  valstr  the values of the old solution
 208 C...  at the relative positions 1/6, 2/6, 4/6 and 5/6 in
 209 C...  each subinterval.
 210 C
 211   140 KSTORE = 1
 212       DO 150 I = 1, N
 213          X = XI(I)
 214          HD6 = (XI(I+1) - XI(I)) / 6.D0
 215          DO 150 J = 1, 4
 216            X = X + HD6
 217            IF ( J.EQ.3 )  X = X + HD6
 218            CALL APPROX (I, X, VALSTR(KSTORE), ASAVE(1,J), DUMMY, XIOLD,
 219      1          NOLD, Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 4, DUMMY, 0)
 220            KSTORE = KSTORE  +  MSTAR
 221   150 CONTINUE
 222   160 MSHFLG = 0
 223       MSHNUM = 1
 224       MODE = 2
 225 C
 226 C...  generate the halved mesh.
 227 C
 228       J = 2
 229       DO 170 I = 1, N
 230            XI(J) = (XIOLD(I) + XIOLD(I+1)) / 2.D0
 231            XI(J+1) = XIOLD(I+1)
 232   170 J = J + 2
 233       N = N2
 234       GO TO 320
 235 C
 236 C...  mode=1  we do mesh selection if it is deemed worthwhile
 237 C
 238   180 IF ( NOLD .EQ. 1 )                            GO TO 100
 239       IF ( NOLD .LE. 2*NFXPNT )                     GO TO 100
 240 C
 241 C...  the first interval has to be treated separately from the
 242 C...  other intervals (generally the solution on the (i-1)st and ith
 243 C...  intervals will be used to approximate the needed derivative, but
 244 C...  here the 1st and second intervals are used.)
 245 C
 246       I = 1
 247       HIOLD = XIOLD(2) - XIOLD(1)
 248       CALL HORDER (1, D1, HIOLD, DMZ, NCOMP, K)
 249       HIOLD = XIOLD(3) - XIOLD(2)
 250       CALL HORDER (2, D2, HIOLD, DMZ, NCOMP, K)
 251       ACCUM(1) = 0.D0
 252       SLOPE(1) = 0.D0
 253       ONEOVH = 2.D0 / ( XIOLD(3) - XIOLD(1) )
 254       DO 190 J = 1, NTOL
 255            JJ = JTOL(J)
 256            JZ = LTOL(J)
 257   190 SLOPE(1) = DMAX1(SLOPE(1),(DABS(D2(JJ)-D1(JJ))*WGTMSH(J)*
 258      1           ONEOVH / (1.D0 + DABS(Z(JZ)))) **ROOT(J))
 259       SLPHMX = SLOPE(1) * (XIOLD(2) - XIOLD(1))
 260       ACCUM(2) = SLPHMX
 261       IFLIP = 1
 262 C
 263 C...  go through the remaining intervals generating  slope
 264 C...  and  accum .
 265 C
 266       DO 210 I = 2, NOLD
 267            HIOLD = XIOLD(I+1) - XIOLD(I)
 268            IF ( IFLIP .EQ. -1 )
 269      1                   CALL HORDER ( I, D1, HIOLD, DMZ, NCOMP, K)
 270            IF ( IFLIP .EQ. 1 )
 271      1                   CALL HORDER ( I, D2, HIOLD, DMZ, NCOMP, K)
 272            ONEOVH = 2.D0 / ( XIOLD(I+1) - XIOLD(I-1) )
 273            SLOPE(I) = 0.D0
 274 C
 275 C...       evaluate function to be equidistributed
 276 C
 277            DO 200 J = 1, NTOL
 278              JJ = JTOL(J)
 279              JZ = LTOL(J)  +  (I-1) * MSTAR
 280   200      SLOPE(I) = DMAX1(SLOPE(I),(DABS(D2(JJ)-D1(JJ))*WGTMSH(J)*
 281      1                ONEOVH / (1.D0 + DABS(Z(JZ)))) **ROOT(J))
 282 C
 283 C...       accumulate approximate integral of function to be
 284 C...       equidistributed
 285 C
 286            TEMP = SLOPE(I) * (XIOLD(I+1)-XIOLD(I))
 287            SLPHMX = DMAX1(SLPHMX,TEMP)
 288            ACCUM(I+1) = ACCUM(I) + TEMP
 289   210 IFLIP = - IFLIP
 290       AVRG = ACCUM(NOLD+1) / DFLOAT(NOLD)
 291       DEGEQU = AVRG / DMAX1(SLPHMX,PRECIS)
 292 C
 293 C...  naccum=expected n to achieve .1x user requested tolerances
 294 C
 295       NACCUM = ACCUM(NOLD+1) + 1.D0
 296       IF ( IPRINT .LT. 0 )  WRITE(IOUT,350) DEGEQU, NACCUM
 297 C
 298 C...  decide if mesh selection is worthwhile (otherwise, halve)
 299 C
 300       IF ( AVRG .LT. PRECIS )                       GO TO 100
 301       IF ( DEGEQU .GE. .5D0 )                       GO TO 100
 302 C
 303 C...  nmx assures mesh has at least half as many subintervals as the
 304 C...  previous mesh
 305 C
 306       NMX = MAX0 ( NOLD+1, NACCUM ) / 2
 307 C
 308 C...  this assures that halving will be possible later (for error est)
 309 C
 310       NMAX2 = NMAX / 2
 311 C
 312 C...  the mesh is at most halved
 313 C
 314       N = MIN0 ( NMAX2, NOLD, NMX )
 315   220 NOLDP1 = NOLD + 1
 316       IF ( N .LT. NFXP1 )  N = NFXP1
 317       MSHNUM = MSHNUM + 1
 318 C
 319 C...  if the new mesh is smaller than the old mesh set mshnum
 320 C...  so that the next call to  newmsh  will produce a halved
 321 C...  mesh. if n .eq. nold / 2 increment mshalt so there can not
 322 C...  be an infinite loop alternating between n and n/2 points.
 323 C
 324       IF ( N .LT. NOLD )  MSHNUM = MSHLMT
 325       IF ( N .GT. NOLD/2 )  MSHALT = 1
 326       IF ( N .EQ. NOLD/2 )  MSHALT = MSHALT + 1
 327       MSHFLG = 0
 328 C
 329 C...  having decided to generate a new mesh with n subintervals we now
 330 C...  do so, taking into account that the nfxpnt points in the array
 331 C...  fixpnt must be included in the new mesh.
 332 C
 333       IN = 1
 334       ACCL = 0.D0
 335       LOLD = 2
 336       XI(1) = ALEFT
 337       XI(N+1) = ARIGHT
 338       DO 310 I = 1, NFXP1
 339            IF ( I .EQ. NFXP1 )                      GO TO 250
 340            DO 230 J = LOLD, NOLDP1
 341              LNEW = J
 342              IF ( FIXPNT(I) .LE. XIOLD(J) )         GO TO 240
 343   230      CONTINUE
 344   240      CONTINUE
 345            ACCR = ACCUM(LNEW) + (FIXPNT(I)-XIOLD(LNEW))*SLOPE(LNEW-1)
 346            NREGN = (ACCR-ACCL) / ACCUM(NOLDP1) * DFLOAT(N) - .5D0
 347            NREGN = MIN0(NREGN, N - IN - NFXP1 + I)
 348            XI(IN + NREGN + 1) = FIXPNT(I)
 349            GO TO 260
 350   250      ACCR = ACCUM(NOLDP1)
 351            LNEW = NOLDP1
 352            NREGN = N - IN
 353   260      IF ( NREGN .EQ. 0 )                      GO TO 300
 354            TEMP = ACCL
 355            TSUM = (ACCR - ACCL) / DFLOAT(NREGN+1)
 356            DO 290 J = 1, NREGN
 357              IN = IN + 1
 358              TEMP = TEMP + TSUM
 359              DO 270 L = LOLD, LNEW
 360                LCARRY = L
 361                IF ( TEMP .LE. ACCUM(L) )            GO TO 280
 362   270        CONTINUE
 363   280        CONTINUE
 364              LOLD = LCARRY
 365   290      XI(IN) = XIOLD(LOLD-1) + (TEMP - ACCUM(LOLD-1)) /
 366      1     SLOPE(LOLD-1)
 367   300      IN = IN + 1
 368            ACCL = ACCR
 369            LOLD = LNEW
 370   310 CONTINUE
 371       MODE = 1
 372   320 CONTINUE
 373       NP1 = N + 1
 374       IF ( IPRINT .LT. 1 )  WRITE(IOUT,340) N, (XI(I),I=1,NP1)
 375       NZ   = MSTAR * (N + 1)
 376       NDMZ = KD * N
 377       RETURN
 378 C----------------------------------------------------------------
 379   340 FORMAT(/17H THE NEW MESH (OF,I5,16H SUBINTERVALS), ,100(/8F12.6))
 380   350 FORMAT(/21H MESH SELECTION INFO,/30H DEGREE OF EQUIDISTRIBUTION =
 381      1       , F8.5, 28H PREDICTION FOR REQUIRED N = , I8)
 382   360 FORMAT(/20H THE FORMER MESH (OF,I5,15H SUBINTERVALS),,
 383      1       100(/8F12.6))
 384   370 FORMAT (/23H  EXPECTED N TOO LARGE  )
 385       END