Return to Buram.e CVS log

Up to [CSRG BSD Unix] / 41BSD / cmd / efl / efltest

Annotation of 41BSD/cmd/efl/efltest/Buram.e, revision 1.1

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.