Add intro and pdf for lognormal
[maxima.git] / share / colnew / fortran / newmsh.f
blob03167f5216a256c33bdb220e6541eea7c7add6d1
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----------------------------------------------------------------------
7 SUBROUTINE NEWMSH (MODE, XI, XIOLD, Z, DMZ, VALSTR,
8 1 SLOPE, ACCUM, NFXPNT, FIXPNT)
10 C**********************************************************************
12 C purpose
13 C select a mesh on which a collocation solution is to be
14 C determined
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
33 C**********************************************************************
35 C variables
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
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**********************************************************************
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)
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
94 NFXP1 = NFXPNT +1
95 GO TO (180, 100, 50, 20, 10), MODE
97 C... mode=5 set mshlmt=1 so that no mesh selection is performed
99 10 MSHLMT = 1
101 C... mode=4 the user-specified initial mesh is already in place.
103 20 IF ( IGUESS .LT. 2 ) GO TO 40
105 C... iguess=2, 3 or 4.
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
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 .
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
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.
129 50 IF ( N .LT. NFXP1 ) N = NFXP1
130 NP1 = N + 1
131 XI(1) = ALEFT
132 ILEFT = 1
133 XLEFT = ALEFT
135 C... loop over the subregions between fixed points.
137 DO 90 J = 1, NFXP1
138 XRIGHT = ARIGHT
139 IRIGHT = NP1
140 IF ( J .EQ. NFXP1 ) GO TO 60
141 XRIGHT = FIXPNT(J)
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)
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
152 C... generate equally spaced points between the j-1st and the
153 C... j-th fixed points.
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
165 C... mode=2 halve the current mesh (i.e. double its size)
167 100 N2 = 2 * N
169 C... check that n does not exceed storage limitations
171 IF ( N2 .LE. NMAX ) GO TO 120
173 C... if possible, try with n=nmax. redistribute first.
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
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 .
188 120 IF ( MSHFLG .EQ. 0 ) GO TO 140
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.
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
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.
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
226 C... generate the halved mesh.
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
236 C... mode=1 we do mesh selection if it is deemed worthwhile
238 180 IF ( NOLD .EQ. 1 ) GO TO 100
239 IF ( NOLD .LE. 2*NFXPNT ) GO TO 100
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.)
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
263 C... go through the remaining intervals generating slope
264 C... and accum .
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
275 C... evaluate function to be equidistributed
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))
283 C... accumulate approximate integral of function to be
284 C... equidistributed
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)
293 C... naccum=expected n to achieve .1x user requested tolerances
295 NACCUM = ACCUM(NOLD+1) + 1.D0
296 IF ( IPRINT .LT. 0 ) WRITE(IOUT,350) DEGEQU, NACCUM
298 C... decide if mesh selection is worthwhile (otherwise, halve)
300 IF ( AVRG .LT. PRECIS ) GO TO 100
301 IF ( DEGEQU .GE. .5D0 ) GO TO 100
303 C... nmx assures mesh has at least half as many subintervals as the
304 C... previous mesh
306 NMX = MAX0 ( NOLD+1, NACCUM ) / 2
308 C... this assures that halving will be possible later (for error est)
310 NMAX2 = NMAX / 2
312 C... the mesh is at most halved
314 N = MIN0 ( NMAX2, NOLD, NMX )
315 220 NOLDP1 = NOLD + 1
316 IF ( N .LT. NFXP1 ) N = NFXP1
317 MSHNUM = MSHNUM + 1
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.
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
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.
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 )