+ subroutine buram(npts, mesh, fn, m, n, p, q, delk)
+ integer npts
+ integer m, n
+ real mesh(npts), fn(npts), p(1), q(1), delk
+ common /dfccom/ nitr
+ integer nitr
+ integer ier, maxitr, nerror, i, itol
+ real fnmin, fnmax
+ logical smonor
+c Buram is a double precision subroutine which finds a
+c a rational function which is the best approximation,
+c in the uniform or minimax sense, to a given discrete
+c function. The rational function is represented as
+c the quotient of two polynomials each expanded in terms
+c of Tchebychev polynomials. This routine is a shell
+c which in turn calls the routine Burm1 with certain
+c default values for the initial approximation and for
+c the stopping criteria.
+c Input:
+c npts - the number of mesh points.
+c mesh - the array of mesh points.
+c fn - the array of function values.
+c m - the degree of the desired numerator polynomial.
+c n - the degree of the desired denominator polynomial.
+c Output:
+c p - the array of coefficients for the numerator polynomial.
+c q - the array of coefficients for the denominator polynomial.
+c delk - the maximum error in the approximation.
+c Error States (asterisk indicates recoverable):
+c 1 - invalid degree
+c 2 - too few mesh points
+c 3 - mesh is not strictly monotone
+c 4* - approximation equals function
+c 5* - no improvement in approximation
+c 6* - reached 50 iterations
+ call enter(1)
+ if (m .lt. 0 .or. n .lt. 0) call seterr(
+ 1 23h Buram - invalid degree, 23, 1, 2)
+ if (npts .lt. m+n+2) call seterr(28h Buram - too few mesh points
+ 1 , 28, 2, 2)
+ if (.not. smonor(mesh, npts, 1)) call seterr(
+ 1 38h Buram - mesh is not strictly monotone, 38, 3, 2)
+c Initialize the numerator and demoninator polynomials.
+ fnmax = fn(1)
+ fnmin = fn(1)
+ do 3 i = 2, npts
+ if (fnmax .ge. fn(i)) goto 1
+ fnmax = fn(i)
+ goto 2
+ 1 if (fn(i) .lt. fnmin) fnmin = fn(i)
+ 2 continue
+ 3 continue
+ call setr(m+1, 0.0e0, p)
+ p(1) = 0.5e0*(fnmax+fnmin)
+ call setr(n+1, 0.0e0, q)
+ q(1) = 1.0e0
+ delk = fnmax-p(1)
+ nitr = 0
+ if (m .eq. 0 .and. n .eq. 0) goto 11
+ maxitr = 50
+ itol = 2
+ call burm1(npts, mesh, fn, maxitr, itol, m, n, p, q, delk)
+ if (nerror(ier) .eq. 0) goto 10
+ if (ier .ne. 7) goto 4
+ call newerr(38h Buram - approximation equals function,
+ 1 39, 4, 1)
+ goto 9
+ 4 if (ier .ne. 8) goto 5
+ call newerr(
+ 1 40h Buram - no improvement in approximation, 40, 5,
+ 2 1)
+ goto 8
+ 5 if (ier .ne. 9) goto 6
+ call newerr(30h Buram - reached 50 iterations, 30
+ 1 , 6, 1)
+ goto 7
+ 6 call eprint
+ 7 continue
+ 8 continue
+ 9 continue
+ 10 continue
+ 11 call leave
+ return
+ end
+ subroutine burm1(npts, mesh, fn, maxitr, itol, m, n, p, q,
+ 1 delk)
+ integer npts
+ integer maxitr, itol, m, n
+ real mesh(npts), fn(npts), p(1), q(1), delk
+ common /cstak/ dstak
+ double precision dstak(500)
+ integer istkgt, iexptr, j, idig, iflr, istak(1000)
+ integer enptr, qkptr, i1mach, npptr, nqptr
+ real abs, qlrg, float, ws(1000), r1mach
+ logical smonor
+ equivalence (dstak, istak)
+ equivalence (dstak, ws)
+c Burm1 is a double precision subroutine which finds a
+c a rational function which is the best approximation,
+c in the uniform or minimax sense, to a given discrete
+c function. The rational function is represented as
+c the quotient of two polynomials each expanded in terms
+c of Tchebychev polynomials. This routine starts from an
+c initial approximation and terminates for one of four
+c reasons: (1) the error curve equioscillates and the
+c alternating extrema match to ITOL digits, (2) the number
+c of iterations exceeds MAXITR, (3) the approximation
+c cannot be improved, or (4) the approximation is essentially
+c equal to the given discrete function.
+c Input:
+c npts - the number of mesh points.
+c mesh - the array of mesh points.
+c fn - the array of function values.
+c maxitr - the maximum number of iterations.
+c itol - the number of digits to which the extrema should match.
+c m - the degree of the desired numerator polynomial.
+c n - the degree of the desired denominator polynomial.
+c p - the array of coefficients for the initial numerator.
+c q - the array of coefficients for the initial denominator.
+c Output:
+c p - the array of coefficients for the numerator polynomial.
+c q - the array of coefficients for the denominator polynomial.
+c delk - the maximum error in the approximation.
+c Error States (asterisk indicates recoverable):
+c 1 - invalid degree
+c 2 - too few mesh points
+c 3 - mesh is not strictly monotone
+c 4 - maxitr .lt. 0
+c 5 - invalid accuracy request
+c 6 - denominator is nonpositive
+c 7* - approximation equals function
+c 8* - no improvement in approximation
+c 9* - reached maximum no. of iterations
+ call enter(1)
+ if (m .lt. 0 .or. n .lt. 0) call seterr(
+ 1 23h Burm1 - invalid degree, 23, 1, 2)
+ if (npts .lt. m+n+2) call seterr(28h Burm1 - too few mesh points
+ 1 , 28, 2, 2)
+ if (.not. smonor(mesh, npts, 1)) call seterr(
+ 1 38h Burm1 - mesh is not strictly monotone, 38, 3, 2)
+ if (maxitr .lt. 0) call seterr(22h Burm1 - maxitr .lt. 0, 22, 4, 2
+ 1 )
+ idig = iflr(r1mach(5)*float(i1mach(11)))
+ if (itol .lt. 1 .or. idig .lt. itol) call seterr(
+ 1 33h Burm1 - invalid accuracy request, 36, 5, 2)
+ qlrg = abs(q(1))
+ j = 2
+ goto 2
+ 1 j = j+1
+ 2 if (j .gt. n+1) goto 3
+ if (qlrg .lt. abs(q(j))) qlrg = abs(q(j))
+ goto 1
+ 3 if (qlrg .ne. 0.e0) goto 4
+ call seterr(35h Burm1 - denominator is nonpositive, 35, 6, 2)
+ goto 11
+ 4 j = 1
+ goto 6
+ 5 j = j+1
+ 6 if (j .gt. n+1) goto 7
+ q(j) = q(j)/qlrg
+ goto 5
+ 7 j = 1
+ goto 9
+ 8 j = j+1
+ 9 if (j .gt. m+1) goto 10
+ p(j) = p(j)/qlrg
+ goto 8
+ 10 continue
+ 11 npptr = istkgt(m+1, 3)
+ nqptr = istkgt(n+1, 3)
+ enptr = istkgt(npts, 3)
+ qkptr = istkgt(npts, 3)
+ iexptr = istkgt(npts, 2)
+ call b1rm1(npts, mesh, fn, maxitr, itol, m, n, p, q, delk, ws(
+ 1 npptr), ws(nqptr), ws(enptr), ws(qkptr), istak(iexptr))
+ call leave
+ return
+ end
+ subroutine b1rm1(npts, x, fn, maxitr, itol, m, n, p, q,
+ 1 delk, newp, newq, en, qk, iext)
+ integer npts
+ integer maxitr, itol, m, n, iext(npts)
+ real x(npts), fn(npts), p(1), q(1), delk, newp(1)
+ real newq(1), en(npts), qk(npts)
+ common /dfccom/ nitr
+ integer nitr
+ integer ier, nex, nerror, i, imin, imax
+ integer ilrg, lrgex
+ real bnd, abs, eps, delnew, r1mach
+ eps = r1mach(4)*10.0e0**itol
+ call extrmr(npts, fn, nex, iext, imax, imin, ilrg)
+ bnd = abs(fn(ilrg))*eps
+ call enqk(npts, x, fn, m, n, p, q, qk, en)
+ do 1 i = 1, npts
+ if (qk(i) .le. 0.0e0) call seterr(
+ 1 35h Burm1 - denominator is nonpositive, 35, 6, 2)
+ 1 continue
+ call extrmr(npts, en, nex, iext, imax, imin, ilrg)
+ delk = abs(en(ilrg))
+ delnew = delk
+ call movefr(m+1, p, newp)
+ call movefr(n+1, q, newq)
+ nitr = 0
+ goto 3
+ 2 nitr = nitr+1
+ 3 if (nitr .ge. maxitr) goto 6
+c call Outpt3 (x,npts,p,q,delk,m,n,en,iext,nex)
+ if (delk .gt. bnd) goto 4
+ call seterr(38h Burm1 - approximation equals function, 39, 7
+ 1 , 1)
+ return
+c Test for optimal solution.
+ 4 if (lrgex(npts, en, nex, iext, ilrg, itol) .ge. m+n+2) return
+ call lpstp(npts, x, fn, qk, delnew, m, n, newp, newq)
+ if (nerror(ier) .ne. 0) call erroff
+ call enqk(npts, x, fn, m, n, newp, newq, qk, en)
+ call extrmr(npts, en, nex, iext, imax, imin, ilrg)
+ delnew = abs(en(ilrg))
+ if (delk .gt. delnew) goto 5
+ call seterr(40h Burm1 - no improvement in approximation, 40,
+ 1 8, 1)
+ return
+ 5 call movefr(m+1, newp, p)
+ call movefr(n+1, newq, q)
+ delk = delnew
+ goto 2
+ 6 call seterr(42h Burm1 - reached maximum no. of iterations, 42, 9
+ 1 , 1)
+ return
+ end
+ subroutine enqk(npts, x, fn, m, n, p, q, qk, en)
+ integer npts
+ integer m, n
+ real x(npts), fn(npts), p(1), q(1), qk(npts), en(npts)
+ integer i
+ real pk, tchbp
+c
+c Subroutine Enqk computes en & Qk.
+c en=error values at mesh points.
+c Qk=value of denominator polynomial at mesh points.
+c
+ if (npts .le. 0 .or. m .lt. 0 .or. n .lt. 0) call seterr(
+ 1 22henQk-invalid dimension, 22, 1, 2)
+ do 1 i = 1, npts
+ qk(i) = tchbp(n, q, x(i), x(1), x(npts))
+ if (qk(i) .eq. 0.e0) call seterr(20henQk-divisor .eq. 0., 20, 2
+ 1 , 2)
+ pk = tchbp(m, p, x(i), x(1), x(npts))
+ en(i) = (fn(i)*qk(i)-pk)/qk(i)
+ 1 continue
+ return
+ end
+ integer function lrgex(npts, en, nex, iext, ilrg, tol)
+ integer nex, npts
+ integer iext(nex), ilrg, tol
+ real en(npts)
+ integer j, k, l
+ real abs, hold
+c
+c Function Lrgex finds the no. of error extrema with magnitudes
+c within tolerance of magnitude of largest error.
+c
+ if (npts .le. 0) call seterr(24hLrgex -invalid dimension, 24, 1, 2
+ 1 )
+ if (nex .le. 0 .or. ilrg .le. 0) call seterr(
+ 1 20hLrgex -invalid index, 20, 2, 2)
+ k = 0
+ do 1 j = 1, nex
+ l = iext(j)
+ hold = abs(en(ilrg))-abs(en(l))
+ if (hold .le. 10.**(-tol)*abs(en(ilrg))) k = k+1
+ 1 continue
+ lrgex = k
+ return
+ end
+ subroutine lpstp(npts, mesh, fn, qk, delk, m, n, p, q)
+ integer npts
+ integer m, n
+ real mesh(npts), fn(npts), qk(npts), delk, p(1), q(1)
+ common /cstak/ dstak
+ double precision dstak(500)
+ integer istkgt, aptr, bptr, cptr, xptr, istak(1000)
+ real ws(1000)
+ equivalence (dstak, istak)
+ equivalence (dstak, ws)
+c Lpstp defines the linear programming subproblem of the
+c Differential Correction algorithm. It also provides
+c the interface to the general purpose linear programming
+c package.
+c Input:
+c npts - the number of mesh points.
+c mesh - the array of mesh points.
+c fn - the array of function values.
+c Qk - the array of current denominator values.
+c delk - the current minimax error.
+c m - the degree of the numerator polynomial.
+c n - the degree of the denominator polynomial.
+c p - the current numerator polynomial.
+c q - the current denominator polynomial.
+c Output:
+c p - the array of coefficients for the numerator polynomial.
+c q - the array of coefficients for the denominator polynomial.
+c Error States (asterisk indicates fatal):
+c 1* - invalid degree
+c 2* - too few mesh points
+c 3* - nonpositive delk
+c 4 - no improvement in the lp subproblem
+ call enter(1)
+ if (m .lt. 0 .or. n .lt. 0) call seterr(
+ 1 23h Lpstp - invalid degree, 23, 1, 2)
+ if (npts .lt. m+n+2) call seterr(28h Lpstp - too few mesh points
+ 1 , 28, 2, 2)
+ aptr = istkgt(3*npts+1, 3)
+ bptr = istkgt(2*(npts+n+1), 3)
+ cptr = istkgt(m+n+3, 3)
+ xptr = istkgt(m+n+3, 3)
+ call l9stp(npts, mesh, fn, qk, delk, m, n, p, q, ws(aptr), ws(
+ 1 bptr), ws(cptr), ws(xptr))
+ call leave
+ return
+ end
+ subroutine l9stp(npts, mesh, fn, qk, delk, m, n, p, q, a, b,
+ 1 c, x)
+ integer npts
+ integer m, n
+ real mesh(npts), fn(npts), qk(npts), delk, p(1), q(1)
+ real a(1), b(1), c(1), x(1)
+ common /difcom/ nptsc, mc, nc, i1, i2, i3, i4
+ integer nptsc, mc, nc, i1, i2, i3
+ integer i4
+ external difmt
+ integer ier, nerror, i, j, ierr, mm
+ integer nn
+ real abs, ctx, ctxnew, qlrg, float, r1mach
+ nptsc = npts
+ mc = m
+ nc = n
+ i1 = npts
+ i2 = i1+npts
+ i3 = i2+n+1
+ i4 = i3+n+1
+ mm = i4
+ nn = m+n+3
+ call movefr(n+1, q, x)
+ call movefr(m+1, p, x(n+2))
+ x(nn) = 0.e0
+ call setr(i2, 0.0e0, b)
+ call setr(i4-i2, -1.0e0, b(i2+1))
+ call setr(nn, 0.0e0, c)
+ c(nn) = -1.0e0
+ call movefr(npts, mesh, a)
+ call movefr(npts, fn, a(npts+1))
+ call movefr(npts, qk, a(2*npts+1))
+ if (delk .le. 0.0e0) call seterr(25h Lpstp - nonpositive delk, 25,
+ 1 3, 2)
+ a(3*npts+1) = delk
+ ctx = 0.0e0
+c Solve the LP problem: max C(T)X subject to AX >= B.
+c The subroutine Difmt derives the matrix A from
+c the data stored in the array A.
+ call lpph2(a, mm, nn, difmt, b, c, x, 4*mm, ctxnew)
+ if (nerror(ier) .ne. 0) call erroff
+ if (ctx .ge. ctxnew) goto 10
+ qlrg = 0.0e0
+ j = 1
+ goto 2
+ 1 j = j+1
+ 2 if (j .gt. n+1) goto 3
+ if (qlrg .lt. abs(x(j))) qlrg = abs(x(j))
+ goto 1
+ 3 j = 1
+ goto 5
+ 4 j = j+1
+ 5 if (j .gt. n+1) goto 6
+ q(j) = x(j)/qlrg
+ goto 4
+ 6 i = 0
+ j = n+2
+ goto 8
+ 7 j = j+1
+ 8 if (j .gt. m+n+2) goto 9
+ i = i+1
+ p(i) = x(j)/qlrg
+ goto 7
+ 9 continue
+ goto 11
+ 10 call seterr(44h Lpstp - no improvement in the lp subproblem,
+ 1 44, 4, 1)
+ 11 return
+ end
+ subroutine difmt(inprod, a, mm, nn, irow, x, dinprd)
+ integer nn
+ integer mm, irow
+ real a(1), x(nn), dinprd
+ logical inprod
+ common /difcom/ npts, m, n, i1, i2, i3, i4
+ integer npts, m, n, i1, i2, i3
+ integer i4
+ integer max0, irm1, irm2, irm3, j, zptr
+ integer jp, maxmn, fnptr, qzptr
+ real fct, delk, z, fn, fdelk, tchbp
+ real qz
+c Difmt handles references by the LP routine to
+c the matrix for the linear programming subproblem.
+ call enter(1)
+ if (mm .ne. i4 .or. nn .ne. m+n+3) call seterr(
+ 1 26h Difmt - invalid dimension, 26, 1, 2)
+ if (irow .lt. 0 .or. mm .lt. irow) call seterr(
+ 1 22h Difmt - invalid index, 22, 2, 2)
+ irm1 = irow-i1
+ irm2 = irow-i2
+ irm3 = irow-i3
+ if ((.not. inprod) .or. i2 .ge. irow) goto 3
+ if (i3 .ge. irow) goto 1
+ dinprd = -x(irm3)
+ goto 2
+ 1 dinprd = x(irm2)
+ 2 continue
+ goto 18
+ 3 if (i2 .ge. irow) goto 6
+ call setr(nn, 0.0e0, x)
+ if (i3 .ge. irow) goto 4
+ x(irm3) = -1.0e0
+ goto 5
+ 4 x(irm2) = 1.0e0
+ 5 continue
+ goto 17
+ 6 if (i1 .ge. irow) goto 7
+ fct = -1.0e0
+ zptr = irm1
+ goto 8
+ 7 fct = 1.0e0
+ zptr = irow
+ 8 z = a(zptr)
+ fnptr = zptr+npts
+ fn = a(fnptr)
+ qzptr = fnptr+npts
+ qz = a(qzptr)
+ delk = a(3*npts+1)
+ fdelk = fct*fn+delk
+ if (.not. inprod) goto 9
+ dinprd = fdelk*tchbp(n, x, z, a(1), a(npts))-fct*tchbp(m,
+ 1 x(n+2), z, a(1), a(npts))+qz*x(nn)
+ goto 16
+ 9 maxmn = max0(m, n)
+ call tchcf(z, a(1), a(npts), maxmn, x)
+ j = m+1
+ goto 11
+ 10 j = j-1
+ 11 if (1 .gt. j) goto 12
+ jp = j+n+1
+ x(jp) = (-fct)*x(j)
+ goto 10
+ 12 j = 1
+ goto 14
+ 13 j = j+1
+ 14 if (j .gt. n+1) goto 15
+ x(j) = fdelk*x(j)
+ goto 13
+ 15 x(nn) = qz
+ 16 continue
+ 17 continue
+ 18 call leave
+ return
+ end
+ subroutine tchcf(x, a, b, deg, xx)
+ integer deg
+ real x, a, b, xx(1)
+ integer i
+ real twoxx
+c
+c Subroutine Tchcf computes the deg+1 Tchebycheff
+c coefficients of the point x.
+c
+ call enter(1)
+ if (deg .lt. 0) call seterr(21h Tchcf-invalid degree, 21, 1, 2)
+ xx(1) = 1.e0
+ if (deg .le. 0) goto 3
+ if (b .gt. a) goto 1
+ call seterr(23h Tchcf-invalid interval, 23, 2, 2)
+ goto 2
+ 1 xx(2) = 2.e0*(x-(a+b)/2.e0)/(b-a)
+cscale x to the interval (-1.e0,1.e0)
+ 2 continue
+ 3 if (deg .gt. 1) twoxx = 2.e0*xx(2)
+ i = 3
+ goto 5
+ 4 i = i+1
+ 5 if (i .gt. deg+1) goto 6
+ xx(i) = twoxx*xx(i-1)-xx(i-2)
+ goto 4
+ 6 call leave
+ return
+ end