![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to Buram.out CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 41BSD / cmd / efl / efltest|
Return to Buram.out CVS log | Up to [CSRG BSD Unix] / 41BSD / cmd / efl / efltest |
1.1 root 1: subroutine buram(npts, mesh, fn, m, n, p, q, delk)
2: integer npts
3: integer m, n
4: real mesh(npts), fn(npts), p(1), q(1), delk
5: common /dfccom/ nitr
6: integer nitr
7: integer ier, maxitr, nerror, i, itol
8: real fnmin, fnmax
9: logical smonor
10: c Buram is a double precision subroutine which finds a
11: c a rational function which is the best approximation,
12: c in the uniform or minimax sense, to a given discrete
13: c function. The rational function is represented as
14: c the quotient of two polynomials each expanded in terms
15: c of Tchebychev polynomials. This routine is a shell
16: c which in turn calls the routine Burm1 with certain
17: c default values for the initial approximation and for
18: c the stopping criteria.
19: c Input:
20: c npts - the number of mesh points.
21: c mesh - the array of mesh points.
22: c fn - the array of function values.
23: c m - the degree of the desired numerator polynomial.
24: c n - the degree of the desired denominator polynomial.
25: c Output:
26: c p - the array of coefficients for the numerator polynomial.
27: c q - the array of coefficients for the denominator polynomial.
28: c delk - the maximum error in the approximation.
29: c Error States (asterisk indicates recoverable):
30: c 1 - invalid degree
31: c 2 - too few mesh points
32: c 3 - mesh is not strictly monotone
33: c 4* - approximation equals function
34: c 5* - no improvement in approximation
35: c 6* - reached 50 iterations
36: call enter(1)
37: if (m .lt. 0 .or. n .lt. 0) call seterr(
38: 1 23h Buram - invalid degree, 23, 1, 2)
39: if (npts .lt. m+n+2) call seterr(28h Buram - too few mesh points
40: 1 , 28, 2, 2)
41: if (.not. smonor(mesh, npts, 1)) call seterr(
42: 1 38h Buram - mesh is not strictly monotone, 38, 3, 2)
43: c Initialize the numerator and demoninator polynomials.
44: fnmax = fn(1)
45: fnmin = fn(1)
46: do 3 i = 2, npts
47: if (fnmax .ge. fn(i)) goto 1
48: fnmax = fn(i)
49: goto 2
50: 1 if (fn(i) .lt. fnmin) fnmin = fn(i)
51: 2 continue
52: 3 continue
53: call setr(m+1, 0.0e0, p)
54: p(1) = 0.5e0*(fnmax+fnmin)
55: call setr(n+1, 0.0e0, q)
56: q(1) = 1.0e0
57: delk = fnmax-p(1)
58: nitr = 0
59: if (m .eq. 0 .and. n .eq. 0) goto 11
60: maxitr = 50
61: itol = 2
62: call burm1(npts, mesh, fn, maxitr, itol, m, n, p, q, delk)
63: if (nerror(ier) .eq. 0) goto 10
64: if (ier .ne. 7) goto 4
65: call newerr(38h Buram - approximation equals function,
66: 1 39, 4, 1)
67: goto 9
68: 4 if (ier .ne. 8) goto 5
69: call newerr(
70: 1 40h Buram - no improvement in approximation, 40, 5,
71: 2 1)
72: goto 8
73: 5 if (ier .ne. 9) goto 6
74: call newerr(30h Buram - reached 50 iterations, 30
75: 1 , 6, 1)
76: goto 7
77: 6 call eprint
78: 7 continue
79: 8 continue
80: 9 continue
81: 10 continue
82: 11 call leave
83: return
84: end
85: subroutine burm1(npts, mesh, fn, maxitr, itol, m, n, p, q,
86: 1 delk)
87: integer npts
88: integer maxitr, itol, m, n
89: real mesh(npts), fn(npts), p(1), q(1), delk
90: common /cstak/ dstak
91: double precision dstak(500)
92: integer istkgt, iexptr, j, idig, iflr, istak(1000)
93: integer enptr, qkptr, i1mach, npptr, nqptr
94: real abs, qlrg, float, ws(1000), r1mach
95: logical smonor
96: equivalence (dstak, istak)
97: equivalence (dstak, ws)
98: c Burm1 is a double precision subroutine which finds a
99: c a rational function which is the best approximation,
100: c in the uniform or minimax sense, to a given discrete
101: c function. The rational function is represented as
102: c the quotient of two polynomials each expanded in terms
103: c of Tchebychev polynomials. This routine starts from an
104: c initial approximation and terminates for one of four
105: c reasons: (1) the error curve equioscillates and the
106: c alternating extrema match to ITOL digits, (2) the number
107: c of iterations exceeds MAXITR, (3) the approximation
108: c cannot be improved, or (4) the approximation is essentially
109: c equal to the given discrete function.
110: c Input:
111: c npts - the number of mesh points.
112: c mesh - the array of mesh points.
113: c fn - the array of function values.
114: c maxitr - the maximum number of iterations.
115: c itol - the number of digits to which the extrema should match.
116: c m - the degree of the desired numerator polynomial.
117: c n - the degree of the desired denominator polynomial.
118: c p - the array of coefficients for the initial numerator.
119: c q - the array of coefficients for the initial denominator.
120: c Output:
121: c p - the array of coefficients for the numerator polynomial.
122: c q - the array of coefficients for the denominator polynomial.
123: c delk - the maximum error in the approximation.
124: c Error States (asterisk indicates recoverable):
125: c 1 - invalid degree
126: c 2 - too few mesh points
127: c 3 - mesh is not strictly monotone
128: c 4 - maxitr .lt. 0
129: c 5 - invalid accuracy request
130: c 6 - denominator is nonpositive
131: c 7* - approximation equals function
132: c 8* - no improvement in approximation
133: c 9* - reached maximum no. of iterations
134: call enter(1)
135: if (m .lt. 0 .or. n .lt. 0) call seterr(
136: 1 23h Burm1 - invalid degree, 23, 1, 2)
137: if (npts .lt. m+n+2) call seterr(28h Burm1 - too few mesh points
138: 1 , 28, 2, 2)
139: if (.not. smonor(mesh, npts, 1)) call seterr(
140: 1 38h Burm1 - mesh is not strictly monotone, 38, 3, 2)
141: if (maxitr .lt. 0) call seterr(22h Burm1 - maxitr .lt. 0, 22, 4, 2
142: 1 )
143: idig = iflr(r1mach(5)*float(i1mach(11)))
144: if (itol .lt. 1 .or. idig .lt. itol) call seterr(
145: 1 33h Burm1 - invalid accuracy request, 36, 5, 2)
146: qlrg = abs(q(1))
147: j = 2
148: goto 2
149: 1 j = j+1
150: 2 if (j .gt. n+1) goto 3
151: if (qlrg .lt. abs(q(j))) qlrg = abs(q(j))
152: goto 1
153: 3 if (qlrg .ne. 0.e0) goto 4
154: call seterr(35h Burm1 - denominator is nonpositive, 35, 6, 2)
155: goto 11
156: 4 j = 1
157: goto 6
158: 5 j = j+1
159: 6 if (j .gt. n+1) goto 7
160: q(j) = q(j)/qlrg
161: goto 5
162: 7 j = 1
163: goto 9
164: 8 j = j+1
165: 9 if (j .gt. m+1) goto 10
166: p(j) = p(j)/qlrg
167: goto 8
168: 10 continue
169: 11 npptr = istkgt(m+1, 3)
170: nqptr = istkgt(n+1, 3)
171: enptr = istkgt(npts, 3)
172: qkptr = istkgt(npts, 3)
173: iexptr = istkgt(npts, 2)
174: call b1rm1(npts, mesh, fn, maxitr, itol, m, n, p, q, delk, ws(
175: 1 npptr), ws(nqptr), ws(enptr), ws(qkptr), istak(iexptr))
176: call leave
177: return
178: end
179: subroutine b1rm1(npts, x, fn, maxitr, itol, m, n, p, q,
180: 1 delk, newp, newq, en, qk, iext)
181: integer npts
182: integer maxitr, itol, m, n, iext(npts)
183: real x(npts), fn(npts), p(1), q(1), delk, newp(1)
184: real newq(1), en(npts), qk(npts)
185: common /dfccom/ nitr
186: integer nitr
187: integer ier, nex, nerror, i, imin, imax
188: integer ilrg, lrgex
189: real bnd, abs, eps, delnew, r1mach
190: eps = r1mach(4)*10.0e0**itol
191: call extrmr(npts, fn, nex, iext, imax, imin, ilrg)
192: bnd = abs(fn(ilrg))*eps
193: call enqk(npts, x, fn, m, n, p, q, qk, en)
194: do 1 i = 1, npts
195: if (qk(i) .le. 0.0e0) call seterr(
196: 1 35h Burm1 - denominator is nonpositive, 35, 6, 2)
197: 1 continue
198: call extrmr(npts, en, nex, iext, imax, imin, ilrg)
199: delk = abs(en(ilrg))
200: delnew = delk
201: call movefr(m+1, p, newp)
202: call movefr(n+1, q, newq)
203: nitr = 0
204: goto 3
205: 2 nitr = nitr+1
206: 3 if (nitr .ge. maxitr) goto 6
207: c call Outpt3 (x,npts,p,q,delk,m,n,en,iext,nex)
208: if (delk .gt. bnd) goto 4
209: call seterr(38h Burm1 - approximation equals function, 39, 7
210: 1 , 1)
211: return
212: c Test for optimal solution.
213: 4 if (lrgex(npts, en, nex, iext, ilrg, itol) .ge. m+n+2) return
214: call lpstp(npts, x, fn, qk, delnew, m, n, newp, newq)
215: if (nerror(ier) .ne. 0) call erroff
216: call enqk(npts, x, fn, m, n, newp, newq, qk, en)
217: call extrmr(npts, en, nex, iext, imax, imin, ilrg)
218: delnew = abs(en(ilrg))
219: if (delk .gt. delnew) goto 5
220: call seterr(40h Burm1 - no improvement in approximation, 40,
221: 1 8, 1)
222: return
223: 5 call movefr(m+1, newp, p)
224: call movefr(n+1, newq, q)
225: delk = delnew
226: goto 2
227: 6 call seterr(42h Burm1 - reached maximum no. of iterations, 42, 9
228: 1 , 1)
229: return
230: end
231: subroutine enqk(npts, x, fn, m, n, p, q, qk, en)
232: integer npts
233: integer m, n
234: real x(npts), fn(npts), p(1), q(1), qk(npts), en(npts)
235: integer i
236: real pk, tchbp
237: c
238: c Subroutine Enqk computes en & Qk.
239: c en=error values at mesh points.
240: c Qk=value of denominator polynomial at mesh points.
241: c
242: if (npts .le. 0 .or. m .lt. 0 .or. n .lt. 0) call seterr(
243: 1 22henQk-invalid dimension, 22, 1, 2)
244: do 1 i = 1, npts
245: qk(i) = tchbp(n, q, x(i), x(1), x(npts))
246: if (qk(i) .eq. 0.e0) call seterr(20henQk-divisor .eq. 0., 20, 2
247: 1 , 2)
248: pk = tchbp(m, p, x(i), x(1), x(npts))
249: en(i) = (fn(i)*qk(i)-pk)/qk(i)
250: 1 continue
251: return
252: end
253: integer function lrgex(npts, en, nex, iext, ilrg, tol)
254: integer nex, npts
255: integer iext(nex), ilrg, tol
256: real en(npts)
257: integer j, k, l
258: real abs, hold
259: c
260: c Function Lrgex finds the no. of error extrema with magnitudes
261: c within tolerance of magnitude of largest error.
262: c
263: if (npts .le. 0) call seterr(24hLrgex -invalid dimension, 24, 1, 2
264: 1 )
265: if (nex .le. 0 .or. ilrg .le. 0) call seterr(
266: 1 20hLrgex -invalid index, 20, 2, 2)
267: k = 0
268: do 1 j = 1, nex
269: l = iext(j)
270: hold = abs(en(ilrg))-abs(en(l))
271: if (hold .le. 10.**(-tol)*abs(en(ilrg))) k = k+1
272: 1 continue
273: lrgex = k
274: return
275: end
276: subroutine lpstp(npts, mesh, fn, qk, delk, m, n, p, q)
277: integer npts
278: integer m, n
279: real mesh(npts), fn(npts), qk(npts), delk, p(1), q(1)
280: common /cstak/ dstak
281: double precision dstak(500)
282: integer istkgt, aptr, bptr, cptr, xptr, istak(1000)
283: real ws(1000)
284: equivalence (dstak, istak)
285: equivalence (dstak, ws)
286: c Lpstp defines the linear programming subproblem of the
287: c Differential Correction algorithm. It also provides
288: c the interface to the general purpose linear programming
289: c package.
290: c Input:
291: c npts - the number of mesh points.
292: c mesh - the array of mesh points.
293: c fn - the array of function values.
294: c Qk - the array of current denominator values.
295: c delk - the current minimax error.
296: c m - the degree of the numerator polynomial.
297: c n - the degree of the denominator polynomial.
298: c p - the current numerator polynomial.
299: c q - the current denominator polynomial.
300: c Output:
301: c p - the array of coefficients for the numerator polynomial.
302: c q - the array of coefficients for the denominator polynomial.
303: c Error States (asterisk indicates fatal):
304: c 1* - invalid degree
305: c 2* - too few mesh points
306: c 3* - nonpositive delk
307: c 4 - no improvement in the lp subproblem
308: call enter(1)
309: if (m .lt. 0 .or. n .lt. 0) call seterr(
310: 1 23h Lpstp - invalid degree, 23, 1, 2)
311: if (npts .lt. m+n+2) call seterr(28h Lpstp - too few mesh points
312: 1 , 28, 2, 2)
313: aptr = istkgt(3*npts+1, 3)
314: bptr = istkgt(2*(npts+n+1), 3)
315: cptr = istkgt(m+n+3, 3)
316: xptr = istkgt(m+n+3, 3)
317: call l9stp(npts, mesh, fn, qk, delk, m, n, p, q, ws(aptr), ws(
318: 1 bptr), ws(cptr), ws(xptr))
319: call leave
320: return
321: end
322: subroutine l9stp(npts, mesh, fn, qk, delk, m, n, p, q, a, b,
323: 1 c, x)
324: integer npts
325: integer m, n
326: real mesh(npts), fn(npts), qk(npts), delk, p(1), q(1)
327: real a(1), b(1), c(1), x(1)
328: common /difcom/ nptsc, mc, nc, i1, i2, i3, i4
329: integer nptsc, mc, nc, i1, i2, i3
330: integer i4
331: external difmt
332: integer ier, nerror, i, j, ierr, mm
333: integer nn
334: real abs, ctx, ctxnew, qlrg, float, r1mach
335: nptsc = npts
336: mc = m
337: nc = n
338: i1 = npts
339: i2 = i1+npts
340: i3 = i2+n+1
341: i4 = i3+n+1
342: mm = i4
343: nn = m+n+3
344: call movefr(n+1, q, x)
345: call movefr(m+1, p, x(n+2))
346: x(nn) = 0.e0
347: call setr(i2, 0.0e0, b)
348: call setr(i4-i2, -1.0e0, b(i2+1))
349: call setr(nn, 0.0e0, c)
350: c(nn) = -1.0e0
351: call movefr(npts, mesh, a)
352: call movefr(npts, fn, a(npts+1))
353: call movefr(npts, qk, a(2*npts+1))
354: if (delk .le. 0.0e0) call seterr(25h Lpstp - nonpositive delk, 25,
355: 1 3, 2)
356: a(3*npts+1) = delk
357: ctx = 0.0e0
358: c Solve the LP problem: max C(T)X subject to AX >= B.
359: c The subroutine Difmt derives the matrix A from
360: c the data stored in the array A.
361: call lpph2(a, mm, nn, difmt, b, c, x, 4*mm, ctxnew)
362: if (nerror(ier) .ne. 0) call erroff
363: if (ctx .ge. ctxnew) goto 10
364: qlrg = 0.0e0
365: j = 1
366: goto 2
367: 1 j = j+1
368: 2 if (j .gt. n+1) goto 3
369: if (qlrg .lt. abs(x(j))) qlrg = abs(x(j))
370: goto 1
371: 3 j = 1
372: goto 5
373: 4 j = j+1
374: 5 if (j .gt. n+1) goto 6
375: q(j) = x(j)/qlrg
376: goto 4
377: 6 i = 0
378: j = n+2
379: goto 8
380: 7 j = j+1
381: 8 if (j .gt. m+n+2) goto 9
382: i = i+1
383: p(i) = x(j)/qlrg
384: goto 7
385: 9 continue
386: goto 11
387: 10 call seterr(44h Lpstp - no improvement in the lp subproblem,
388: 1 44, 4, 1)
389: 11 return
390: end
391: subroutine difmt(inprod, a, mm, nn, irow, x, dinprd)
392: integer nn
393: integer mm, irow
394: real a(1), x(nn), dinprd
395: logical inprod
396: common /difcom/ npts, m, n, i1, i2, i3, i4
397: integer npts, m, n, i1, i2, i3
398: integer i4
399: integer max0, irm1, irm2, irm3, j, zptr
400: integer jp, maxmn, fnptr, qzptr
401: real fct, delk, z, fn, fdelk, tchbp
402: real qz
403: c Difmt handles references by the LP routine to
404: c the matrix for the linear programming subproblem.
405: call enter(1)
406: if (mm .ne. i4 .or. nn .ne. m+n+3) call seterr(
407: 1 26h Difmt - invalid dimension, 26, 1, 2)
408: if (irow .lt. 0 .or. mm .lt. irow) call seterr(
409: 1 22h Difmt - invalid index, 22, 2, 2)
410: irm1 = irow-i1
411: irm2 = irow-i2
412: irm3 = irow-i3
413: if ((.not. inprod) .or. i2 .ge. irow) goto 3
414: if (i3 .ge. irow) goto 1
415: dinprd = -x(irm3)
416: goto 2
417: 1 dinprd = x(irm2)
418: 2 continue
419: goto 18
420: 3 if (i2 .ge. irow) goto 6
421: call setr(nn, 0.0e0, x)
422: if (i3 .ge. irow) goto 4
423: x(irm3) = -1.0e0
424: goto 5
425: 4 x(irm2) = 1.0e0
426: 5 continue
427: goto 17
428: 6 if (i1 .ge. irow) goto 7
429: fct = -1.0e0
430: zptr = irm1
431: goto 8
432: 7 fct = 1.0e0
433: zptr = irow
434: 8 z = a(zptr)
435: fnptr = zptr+npts
436: fn = a(fnptr)
437: qzptr = fnptr+npts
438: qz = a(qzptr)
439: delk = a(3*npts+1)
440: fdelk = fct*fn+delk
441: if (.not. inprod) goto 9
442: dinprd = fdelk*tchbp(n, x, z, a(1), a(npts))-fct*tchbp(m,
443: 1 x(n+2), z, a(1), a(npts))+qz*x(nn)
444: goto 16
445: 9 maxmn = max0(m, n)
446: call tchcf(z, a(1), a(npts), maxmn, x)
447: j = m+1
448: goto 11
449: 10 j = j-1
450: 11 if (1 .gt. j) goto 12
451: jp = j+n+1
452: x(jp) = (-fct)*x(j)
453: goto 10
454: 12 j = 1
455: goto 14
456: 13 j = j+1
457: 14 if (j .gt. n+1) goto 15
458: x(j) = fdelk*x(j)
459: goto 13
460: 15 x(nn) = qz
461: 16 continue
462: 17 continue
463: 18 call leave
464: return
465: end
466: subroutine tchcf(x, a, b, deg, xx)
467: integer deg
468: real x, a, b, xx(1)
469: integer i
470: real twoxx
471: c
472: c Subroutine Tchcf computes the deg+1 Tchebycheff
473: c coefficients of the point x.
474: c
475: call enter(1)
476: if (deg .lt. 0) call seterr(21h Tchcf-invalid degree, 21, 1, 2)
477: xx(1) = 1.e0
478: if (deg .le. 0) goto 3
479: if (b .gt. a) goto 1
480: call seterr(23h Tchcf-invalid interval, 23, 2, 2)
481: goto 2
482: 1 xx(2) = 2.e0*(x-(a+b)/2.e0)/(b-a)
483: cscale x to the interval (-1.e0,1.e0)
484: 2 continue
485: 3 if (deg .gt. 1) twoxx = 2.e0*xx(2)
486: i = 3
487: goto 5
488: 4 i = i+1
489: 5 if (i .gt. deg+1) goto 6
490: xx(i) = twoxx*xx(i-1)-xx(i-2)
491: goto 4
492: 6 call leave
493: return
494: end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.