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