![[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] / 40BSD / cmd / efl / efltest|
Return to Buram.out CVS log | Up to [CSRG BSD Unix] / 40BSD / cmd / efl / efltest |
BSD 4.0
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
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.