![[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 |
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.