C.IRATCU........D....................................................... C C FUNCTION - RATIONAL CHEBYCHEV APPROXIMATION TO C THE FUNCTION F C USAGE - CALL IRATCU (F,PHI,G,AA,BB,L,M,P,Q,WK,IER) C PARAMETERS F - USER SUPPLIED EXTERNAL FUNCTION C SUBPROGRAM F(X). F IS TO BE APPROXIMATED C OVER THE INTERVAL (AA,BB) C PHI - VARIABLE TRANSFORMATION WHICH MUST BE C CONTINUOUS AND MONOTONIC. IT IS SPECIFIED C BY A USER SUPPLIED EXTERNAL FUNCTION C SUBPROGRAM C G - WEIGHT FUNCTION WHICH MUST BE CONTINUOUS C AND NON VANISHING IN (AA,BB). C IT IS SPECIFIED BY A USER SUPPLIED EXTERNAL C FUNCTION SUBPROGRAM C AA - LOWER END POINT OF INTERVAL (AA,BB) C BB - UPPER END POINT OF INTERVAL (AA,BB) C L - DEGREE OF NUMERATOR OF RATIONAL APPROXIMATION C M - DEGREE OF DENOMINATOR OF RATIONAL C APPROXIMATION. M MUST NOT BE GREATER THAN L C P - OUTPUT ARRAY P(1),P(2),...,P(L+1) OF C NUMERATOR COEFFICIENTS C Q - OUTPUT ARRAY Q(1),Q(2),...,Q(M+1) OF C DENOMINATOR COEFFICIENTS C WK - WORK AREA OF DIMENSION .GE. (NP1+6)*NP1 C WHERE NP1=L+M+2 C IER - ERROR PARAMETER C TERMINAL ERROR = 128+N C N = 1 INDICATES A POLE IN THE TRIAL C RATIONAL APPROXIMANT DUE TO A SEQ. OF C CRITICAL POINTS WHICH IS NOT MONOTONIC C N = 2 INDICATES A POLE IN THE TRIAL C RATIONAL APPROXIMANT DUE TO A LARGE C DEL(X(I)) WHERE X(I) IS A TRIAL CRITICAL C POINT C WARNING ERROR = 32+N C N = 1 INDICATES ALGORITHM DID NOT CONVERGE C WITHIN 10 ITERATIONS C PRECISION - DOUBLE C REQ'D IMSL ROUTINES - UERTST C AUTHOR/IMPLEMENTOR - O.G. JOHNSON/C.L. SMITH C LANGUAGE - FORTRAN C....................................................................... C LATEST REVISION - MARCH 29, 1971 C SUBROUTINE IRATCU (F,PHI,G,AA,BB,L,M,P,Q,WK,IER) C DIMENSION P(1),Q(1),WK(1) REAL*8 P,Q,A,WK,AA,BB,XLB,DN,DNP1,XXI REAL*8 BB1,XXM1,ZZ,TEST,TEMPL,U,Z1,Y2,H1,Y3,PZ,T1,T2 REAL*8 T3,T4,DEL,Z2,Z,Y,TEMP1,TEMPL1,SUM,PHI,F,G INTEGER V DATA PI,EPS,EQUIB/3.1415926535897932,.0001,1.E-12/ IER = 0 XLB = 0. ITNO = 1 LP1 = L+1 LPM = L+M N = LPM+1 DN = N DNP1 = DN+1.0 RTEST=0. NP1 = N+1 MP1 = M+1 C SET UP INDEXING FOR WORK AREA IB1 = NP1 IB2 = IB1+NP1 IB3 = IB2 + NP1 IB4 = IB3 + NP1 IB5 = IB4 + NP1 IB6 = IB5 + NP1 DO 5 I=1,LP1 P(I) = 0. 5 CONTINUE DO 10 I=1,MP1 Q(I) = 0.0 10 CONTINUE WK(1) = 1.0 DO 15 I=1,NP1 WK(I+1) = -WK(I) 15 CONTINUE Q(1) = 1.0 C COMPUTE INITIAL GUESSES C OF CRITICAL POINTS XXI = 0.0 KK = N/2 BB1 = (BB-AA)*.5 XXM1 = (BB+AA)*.5 DO 20 I=1,KK XXI = XXI + 1. ZZ = DCOS(PI*(1.-XXI/N))*BB1 IJK4 = IB4+I+1 WK(IJK4) = ZZ+XXM1 J = N-I+1 IJK4 = IB4+J WK(IJK4) = XXM1-ZZ 20 CONTINUE WK(IB4+1) = AA WK(IB5) = BB 25 DO 30 I=1,NP1 IJK4 = IB4+I IJK5 = IB5+I WK(IJK5) = PHI(WK(IJK4)) 30 CONTINUE C SET UP AND SOLVE LINEAR EQUATIONS KRET = 1 GO TO 45 35 XLB = (WK(IB2)+XLB*DN)/DNP1 IF (M.EQ.0) GO TO 135 DO 120 I=1,3 KRET = 2 GO TO 45 40 TEST = DABS(XLB-WK(IB2)) XLB = (WK(IB2)+XLB*DN)/DNP1 IF (TEST .LE. .01*DABS(XLB)) GO TO 125 GO TO 120 45 KK = L+2 C SET UP LINEAR EQUATIONS DO 70 V=1,NP1 IJK = IB6+V WK(IJK) = 1.0 IJK5 = IB5+V IF (L.LE.0) GO TO 55 DO 50 J=2,LP1 IC1 = IB6+(J-1)*NP1+V IC2 = IC1-NP1 WK(IC1)=WK(IC2)*WK(IJK5) 50 CONTINUE 55 IJK1 = IB1 + V IJK4 = IB4+V WK(IJK1) = F(WK(IJK4)) IF (M.LE.0) GO TO 65 IJK0 = V+1 TEMPL = XLB*G(WK(IJK4))*WK(IJK0)-WK(IJK1) IC1 = IB6+(KK-1)*NP1+V WK(IC1) = WK(IJK5)*TEMPL DO 60 J=KK,LPM IC1 = IB6+J*NP1+V IC2 = IC1-NP1 WK(IC1) = WK(IC2)*WK(IJK5) 60 CONTINUE IC1 = IB6+NP1*(NP1-1)+V WK(IC1) = WK(IJK0)*G(WK(IJK4)) 65 CONTINUE 70 CONTINUE C SOLVE LINEAR EQUATIONS USING GAUSIAN C ELIMINATION WITH PARTIAL PIVOTING DO 105 NN=1,NP1 IC3 = IB6+(NN-1)*NP1+NN TEMP1 = DABS(WK(IC3)) II = NN IJK1 = IB1+NN IF (NN.EQ.NP1) GO TO 85 NNP1 = NN+1 DO 75 KK=NNP1,NP1 IC1 = IB6+(NN-1)*NP1+KK TEMPL1 = DABS(WK(IC1)) IF (TEMP1.GE.TEMPL1) GO TO 75 II = KK TEMP1 = TEMPL1 75 CONTINUE IF (II.EQ.NN) GO TO 85 DO 80 MM=NN,NP1 IC1 = IB6+(MM-1)*NP1+NN IC2 = IC1-NN+II TEMP1 = WK(IC1) WK(IC1) = WK(IC2) WK(IC2) = TEMP1 80 CONTINUE TEMP1 = WK(IJK1) IJK2 = IB1+II WK(IJK1) = WK(IJK2) WK(IJK2) = TEMP1 85 WK(IJK1) = WK(IJK1)/WK(IC3) DO 90 JJ=NN,NP1 JK = NP1-JJ+NN IC1 = IB6+(JK-1)*NP1+NN WK(IC1) = WK(IC1)/WK(IC3) 90 CONTINUE IF (NN.EQ.NP1) GO TO 105 DO 100 II=NNP1,NP1 IC1 = IB6+(NN-1)*NP1+II DO 95 JJ=NNP1,NP1 IC2 = IB6+(JJ-1)*NP1+II IC4 = IC2-II+NN WK(IC2) = WK(IC2)-WK(IC1)*WK(IC4) 95 CONTINUE IJKN = IB1+II WK(IJKN) = WK(IJKN)-WK(IC1)*WK(IJK1) 100 CONTINUE 105 CONTINUE DO 115 II=1,NP1 JJ = NP1-II+1 IJK1 = IB1+JJ IC1 = IB6+(JJ-1)*NP1+JJ WK(IJK1) = WK(IJK1)/WK(IC1) JJM1 = JJ-1 IF (JJ.EQ.1) GO TO 115 DO 110 KK=1,JJM1 KJ = JJM1-KK+1 IJKN = IB1+KJ IC2 = IC1-JJ+KJ WK(IJKN) = WK(IJKN)-WK(IC2)*WK(IJK1) 110 CONTINUE 115 CONTINUE GO TO (35,40), KRET 120 CONTINUE 125 DO 130 I=1,M IJK1 = IB1+LP1+I Q(I+1) = WK(IJK1) 130 CONTINUE 135 CONTINUE DO 140 I=1,LP1 IJK1 = IB1+I P(I) = WK(IJK1) 140 CONTINUE 145 U = DSIGN(1.D0,XLB) Z1 = 0.0 C SEARCH FOR NEW CRITICAL POINTS IF (N.GT.1) GO TO 150 WK(IB3+1) = .015D0*(WK(IB4+2)-WK(IB4+1)) WK(IB3+2) = -WK(IB3+1) GO TO 160 150 DO 155 I=2,N IJK3 = IB3+I IJK4 = IB4+I+1 IJKN = IB4+I-1 WK(IJK3) = .015D0*(WK(IJK4)-WK(IJKN)) 155 CONTINUE WK(IB3+1) = WK(IB3+2)*.5 IJK3 = IB3+N WK(IB4) = -WK(IJK3)*.5 160 DO 240 I=1,NP1 IJK4 = IB4+I Y2 = WK(IJK4) IJK3 = IB3+I H1 = WK(IJK3) Y3 = Y2+H1 PZ = Y2 KRET = 1 C COMPUTE APPROXIMATION ERROR AT PZ 165 T1 = F(PZ) T2 = PHI(PZ) C ESTIMATE MAXIMUM VALUE OF F T5=DABS(T1) IF(T5.GT.RTEST) RTEST=T5 IJK1 = IB1+LP1 T3 = WK(IJK1) DO 170 II=1,L KD = LP1-II IJK1 = IB1+KD T3 = T3*T2+WK(IJK1) 170 CONTINUE T4 = 0. DO 175 II=1,M KD = NP1-II IJK1 = IB1+KD T4 = (T4+WK(IJK1))*T2 175 CONTINUE T4 = T4+1.0 T2 = T3/T4 DEL = (T2-T1)/G(PZ) C DEL IS THE APPROXIMATION ERROR AT PZ C RETURN TO THE CALLING SEGEMENT GO TO (180,185,205,225,250,275),KRET C DIRECT SEARCH FOR ADJUSTED CRITICAL C POINTS 180 Z2 = DEL*U PZ = Y3 KRET = 2 GO TO 165 185 Z3 = DEL*U IF (Z3.GT.Z2) GO TO 190 H1 = -H1 Z = Z3 Z3 = Z2 Z2 = Z Y = Y3 Y3 = Y2 Y2 = Y 190 Y = Y3+H1 IF (Y.GE.AA) GO TO 195 Y = AA GO TO 220 195 IF (Y.LE.BB) GO TO 200 Y = BB GO TO 220 200 PZ = Y KRET = 3 GO TO 165 205 Z = DEL*U IF (Z.LE.Z3) GO TO 210 Y2 = Y3 Y3 = Y Z2 = Z3 Z3 = Z GO TO 190 210 Y = -Z3-Z3+Z+Z2 IF(Y.EQ.0) GO TO 215 Y = (Y2+Y3)*.5+H1*(Z2-Z3)/Y GO TO 220 215 Y = Y3 220 WK(IJK4) = Y PZ = Y KRET = 4 GO TO 165 225 IJK2 = IB2+I WK(IJK2) = DEL U = -U C CHECK FOR PROPER ORDERING OF C ADJUSTED CRITICAL POINTS IF (.NOT.(I.GT.1.AND.WK(IJK4).LE.WK(IJK4-1))) GO TO 230 IER = 129 GO TO 9000 230 IJK2 = IB2+I Z = DABS(WK(IJK2)) C CHECK FOR POLES IN TRIAL C RATIONAL APPROXIMANT IF (Z.LT.10.) GO TO 235 IER = 130 GO TO 9000 235 Y = DABS(XLB) Z = DABS(Z-Y)/Y IF (Z1.LT.Z) Z1=Z 240 CONTINUE C ADJUST END CRITICAL POINTS SO THAT C THEY DO NOT LIE OUTSIDE THE C INTERVAL (AA,BB) IJK4 = IB4+1 IF (WK(IJK4).LE.AA) GO TO 265 H1 = (WK(IJK4)-AA)*.0625 U = -DSIGN(1.D0,XLB) Z3 = 0.0 Y = AA I = 1 245 PZ = Y KRET = 5 GO TO 165 250 Z = DEL*U IF (Z.LE.Z3) GO TO 255 Z3 = Z Z2 = Y 255 Y = Y+H1 I = I+1 IF (I.LE.16) GO TO 245 Z = DABS(XLB) IF (Z3.LE.Z) GO TO 265 DO 260 I=2,NP1 J = NP1-I+2 IJK2 = IB2+J WK(IJK2) = WK(IJK2-1) IJK4 = IB4+J WK(IJK4) = WK(IJK4-1) 260 CONTINUE IJK4 = IB4+1 WK(IJK4) = Z2 IJK2 = IB2+1 WK(IJK2) = Z3*U GO TO 290 265 IF (BB.LE.WK(IB5)) GO TO 295 H1 = .0625*(BB-WK(IB5)) Z3 = 0.0 Y = BB U = -DSIGN(1.D0,WK(IB3)) I = 1 270 PZ = Y KRET = 6 GO TO 165 275 Z = DEL*U IF (Z3.GE.Z) GO TO 280 Z3 = Z Z2 = Y 280 Y = Y-H1 I = I+1 IF (I.LE.16) GO TO 270 Z = DABS(XLB) IF (Z.GE.Z3) GO TO 295 DO 285 I=1,N IJK2 = IB2+I WK(IJK2) = WK(IJK2+1) IJK4 = IB4+I WK(IJK4) = WK(IJK4+1) 285 CONTINUE WK(IB5) = Z2 WK(IB3) = Z3*U 290 XLB = -XLB Z = DABS(Z3-Z)/Z IF (Z1.LT.Z) Z1=Z 295 CONTINUE C CHECK FOR NO CONVERGENCE IF (ITNO.GE.10) GO TO 315 C CHECK FOR CONVERGENCE IF (Z1.LE.EPS) GO TO 9005 C TEST FOR SIGNIFICIENT DIGITS IN Z1 IF(DABS(XLB).LT. RTEST*EQUIB) GO TO 9005 SUM=0. DO 310 I=1,NP1 IJK2 = IB2+I SUM = SUM+WK(I)*WK(IJK2) 310 CONTINUE XLB = SUM/DNP1 ITNO = ITNO+1 GO TO 25 315 IER = 33 9000 CONTINUE CALL UERTST(IER,'IRATCU') 9005 RETURN END