SUBROUTINE BS(A,P,MM,M,N,PT) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(MM,N),P(M),PT(N) DO 3 IR=1,N S2=0.0 DO 11 I=IR,M 11 S2=S2+P(I)*A(I,IR) S2=-S2/(A(IR,IR)*PT(IR)) DO 12 I=IR,M 12 P(I)=P(I)-S2*A(I,IR) 3 CONTINUE NM=N-1 P(N)=P(N)/PT(N) DO 15 II=2,N I=N-II+1 DO 14 JJ=I,NM J=N-JJ+I 14 P(I)=P(I)-A(I,J)*P(J) 15 P(I)=P(I)/PT(I) RETURN END SUBROUTINE CFDIRZ(Z,A,B2,LEN,TERM,VALZ) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LEN),B2(LEN) COMPLEX*16 Z,TERM,VALZ VALZ=TERM DO 1 LEV=LEN,1,-1 VALZ=B2(LEV)/(Z-A(LEV)-VALZ) 1 CONTINUE RETURN END SUBROUTINE CFGEN(X,W,N,EPS,A,B2,LM1,WK) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(N),W(N),A(LM1),B2(LM1),WK(N,2) SA=0 SB=0 DO 1 I=1,N SB=SB+W(I) SA=SA+X(I)*W(I) WK(I,1)=0 WK(I,2)=1 1 CONTINUE B2(1)=SB A(1)=SA/SB IF(LM1.LT.2)RETURN DO 3 L=2,LM1 ANRM=SQRT(SB) SA=0 SB=0 DO 2 I=1,N DUM=WK(I,2) WK(I,2)=((X(I)-A(L-1))*WK(I,2)-B2(L-1)*WK(I,1))/ANRM WK(I,1)=DUM/ANRM DUM=WK(I,2)*WK(I,2)*W(I) SB=SB+DUM SA=SA+X(I)*DUM 2 CONTINUE B2(L)=SB IF(SB.LT.EPS)GOTO 4 3 A(L)=SA/SB L=LM1+1 4 LM1=L-1 RETURN END FUNCTION DENIN(E,A,B2,NA,NP,LL,ALP,EPS,WK,IWK,ICODE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(NA,NP),B2(NA,NP),WK(NA,2),IWK(NA,2) 1,POLY(2,3) ICODE=0 SUM=0 ICD=1 DO 2 II=1,NP IC=-1 NSTART=1 DUM=RECWT(E,A(1,II),B2(1,II),LL,EPS,IC,POLY,NSTART) SUM=SUM+ALP*DUM LEN=LL+IC NE=0 CALL RECRTS(A(1,II),B2(1,II),LEN,EPS,E-EPS,NE,WK,IWK,WK(1,2), 1 IWK(1,2)) IF (NE.LT.1) GOTO 2 DO 1 I=1,NE IC=0 IF ( ABS(IWK(I,1)) .NE. 1) ICD=-1 WK(I,2)=RECWT(WK(I,1),A(1,II),B2(1,II),LEN,EPS,IC,POLY,NSTART) SUM=SUM+ABS(FLOAT(IWK(I,1)))*WK(I,2) 1 CONTINUE 2 CONTINUE ICODE=(NE+1)*ICD WK(ICODE,1)=E WK(ICODE,2)=DUM DENIN=SUM RETURN END FUNCTION DENQD(E,EMX,A,B2,LL,ALP,EPS,TB,NT,NQ,NE,IWK) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(LL),B2(LL),TB(NT,4),IWK(LL) DOUBLE PRECISION P(2),PP(2),PPP(2),FI(2),FIP(2),XW,SC,DUM N=-1 DENQD=0 DUM=RECWT(E,A,B2,LL,EPS,N,TB,1) NN=LL+N NR=0 CALL RECRTS(A,B2,NN,EPS,EMX+EPS,NR,TB,IWK,TB(1,2),TB(1,3)) IC=1 NE=0 DO 5 NQ=1,NR IF ( IWK(NQ) .NE. 1) RETURN IF ( TB(NQ,1) .GT. EMX+EPS) GOTO 6 XV=TB(NQ,1) P(1)=1 P(2)=XV-A(1) PP(1)=0 PP(2)=1 PPP(1)=0 PPP(2)=0 FI(1)=0 FI(2)=B2(1) FIP(1)=0 FIP(2)=0 I=2 DO 4 J=2,NN XW=XV-A(J) I1=3-I P(I1)=XW*P(I)-B2(J)*P(I1) PP(I1)=XW*PP(I)-B2(J)*PP(I1)+P(I) PPP(I1)=XW*PPP(I)-B2(J)*PPP(I1)+2.0*PP(I) FI(I1)=XW*FI(I)-B2(J)*FI(I1) FIP(I1)=XW*FIP(I)-B2(J)*FIP(I1)+FI(I) SC=ABS(P(I1))+ABS(PP(I1))+ABS(PPP(I1))+ABS(FI(I1)) 1 +ABS(FIP(I1)) DO 3 K=1,2 P(K)=P(K)/SC PP(K)=PP(K)/SC PPP(K)=PPP(K)/SC FI(K)=FI(K)/SC FIP(K)=FIP(K)/SC 3 CONTINUE I=I1 4 CONTINUE I2=3-I TB(NQ,2)=FI(I)/PP(I) P2=PP(I)*PP(I) TB(NQ,4)=P(I2)/PP(I) TB(NQ,3)=(FI(I)*PP(I2)-PP(I)*FI(I2)+TB(NQ,4) 1 *(PP(I)*FIP(I)-FI(I)*PPP(I)))/P2 5 CONTINUE NQ=NR+1 6 NQ=NQ-1 DEV=ABS(TB(1,1)-E) DEN=0 I=2 IF ( NQ .EQ. 1 ) GOTO 8 DO 7 I=2,NQ WD=ABS(TB(I,1)-E) IF ( DEV .LT. WD ) GOTO 8 DEV=WD DEN=DEN+TB(I-1,3) 7 CONTINUE I=NQ+1 8 NE=I-1 DENQD=(DEN+TB(NE,3)*ALP)/TB(NE,4) IF ( DEV .GT. 10.0*EPS ) NE=-NE RETURN END SUBROUTINE DOSCAL(E,A,B,LEN,DOS,WK,NW) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LEN),B(LEN),WK(NW,4) PI=4*ATAN(1.0D0) WK(1,3)=1 WK(2,3)=(E-A(1))/B(2) WK(1,4)=0 WK(2,4)=1/B(2) DO 1 L = 3, LEN WK(L,3)= ((E-A(L-1))*WK(L-1,3) - B(L-1)*WK(L-2,3)) / B(L) WK(L,4)= ((E-A(L-1))*WK(L-1,4) - B(L-1)*WK(L-2,4)) / B(L) 1 CONTINUE WK(1,2)=0 WK(2,2)=1/((WK(1,3)*WK(2,4) - WK(2,3)*WK(1,4))**2) DO 3 L = 3, LEN SUM = 0 DO 2 I = 1, L-1 SUM = SUM + (WK(I,3)*WK(L,4) - WK(L,3)*WK(I,4))**2 2 CONTINUE WK(L,2) = WK(L-1,2) / (1 + WK(L-1,2)*SUM) 3 CONTINUE WK(1,1)=1/WK(1,3)**2 DO 4 I=2,LEN WK(I,1)=WK(I-1,1)/(1+WK(I-1,1)*WK(I,3)**2) 4 CONTINUE DOS=B(1)*B(1)*WK(LEN,1)/(PI*SQRT(WK(LEN,2))) RETURN END SUBROUTINE DYCFTL(LEN,LFIT,CFTL,EMAT,NE,EDIAG,DRATE,WK) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(LEN),B(LEN),CFTL(LEN),EMAT(NE,2),EDIAG(NE),WK(LFIT) LSTART=LEN-LFIT+1 C C fit the asymptotic parameters C DO 1 I=1,LFIT NVAL=I+LSTART-1 WK(I)= LOG(ABS(CFTL(NVAL))) 1 CONTINUE CALL BS(EMAT,WK,LFIT,LFIT,2,EDIAG) DRATE=WK(2) RETURN END SUBROUTINE GRASSZ(Z,A,B,LEN,GRVAL,WK,NW) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LEN),B(LEN) COMPLEX*16 Z,GRVAL(2),WK(NW,2) COMPLEX*16 AQ,BQ,CQ C C compute polynomials C WK(1,1)=(1.0D0,0.0D0)/B(1) WK(2,1)=(Z-A(1))/(B(2)*B(1)) WK(1,2)=(0.0D0,0.0D0) WK(2,2)=(1.0D0,0.0D0)*B(1)/B(2) DO 1 L = 3, LEN WK(L,1)= ((Z-A(L-1))*WK(L-1,1) - B(L-1)*WK(L-2,1)) / B(L) WK(L,2)= ((Z-A(L-1))*WK(L-1,2) - B(L-1)*WK(L-2,2)) / B(L) 1 CONTINUE C C compute Greenian C L = LEN AQ = WK(L-1,1)**2 - WK(L-2,1)*WK(L,1) BQ =WK(L,1)*WK(L-2,2) - 2*WK(L-1,1)*WK(L-1,2) + WK(L-2,1)*WK(L,2) CQ = WK(L-1,2)**2 - WK(L-2,2)*WK(L,2) GRVAL(1) = ( -BQ - (0.0D0,1.0D0)*SQRT(4*AQ*CQ-BQ*BQ)) / (2*AQ) GRVAL(2) = ( -BQ + (0.0D0,1.0D0)*SQRT(4*AQ*CQ-BQ*BQ)) / (2*AQ) RETURN END SUBROUTINE GRCALZ(Z,A,B,LEN,GRVAL,WK,WKR,NW) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LEN),B(LEN),WKR(NW,2) COMPLEX*16 Z,GRVAL,WK(NW,3) C C compute polynomials C WK(1,1)=(1.0D0,0.0D0)/B(1) WK(2,1)=(Z-A(1))/(B(2)*B(1)) WK(1,2)=(0.0D0,0.0D0) WK(2,2)=(1.0D0,0.0D0)*B(1)/B(2) DO 1 L = 3, LEN WK(L,1)= ((Z-A(L-1))*WK(L-1,1) - B(L-1)*WK(L-2,1)) / B(L) WK(L,2)= ((Z-A(L-1))*WK(L-1,2) - B(L-1)*WK(L-2,2)) / B(L) 1 CONTINUE C C compute Greenian C PINV=1/(WK(1,1)*CONJG(WK(1,1))+WK(2,1)*CONJG(WK(2,1))) QINV=1/(WK(2,2)*CONJG(WK(2,2))) T = REAL(WK(2,1)*CONJG(WK(2,2))) WKR(1,1) = B(1)*B(1) WKR(2,1) = PINV WKR(2,2) = QINV DO 2 L=3,LEN PINV = PINV/(1+WK(L,1)*CONJG(WK(L,1))*PINV) QINV = QINV/(1+WK(L,2)*CONJG(WK(L,2))*QINV) T = T + REAL(WK(L,1)*CONJG(WK(L,2))) DRT = T*PINV WKR(L,1)=PINV WKR(L,2)=QINV WK(L,3) = (DRT - (0.0D0,1.0D0)*SQRT(ABS(PINV/QINV - DRT*DRT))) 2 CONTINUE GRVAL=WK(LEN,3) RETURN END SUBROUTINE GRGENZ(Z,A,B,LEN,GRVAL,WK,WKR,NW) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LEN),B(LEN),WKR(NW,3) COMPLEX*16 Z,GRVAL,WK(NW,2) COMPLEX*16 TERM1,TERM2,TERM3,TERM4 C C compute polynomials C WK(1,1)=(1.0D0,0.0D0) WK(2,1)=(Z-A(1))/B(2) WK(1,2)=(0.0D0,0.0D0) WK(2,2)=(1.0D0,0.0D0)/B(2) DO 1 L = 3, LEN WK(L,1)= ((Z-A(L-1))*WK(L-1,1) - B(L-1)*WK(L-2,1)) / B(L) WK(L,2)= ((Z-A(L-1))*WK(L-1,2) - B(L-1)*WK(L-2,2)) / B(L) 1 CONTINUE C C compute inverse outer product C WKR(1,2) = 0.0D0 TERM1 = WK(2,1)*CONJG(WK(1,2)) - CONJG(WK(1,1))*WK(2,2) TERM2 = CONJG(WK(2,1))*CONJG(WK(1,2)) - 1 CONJG(WK(1,1))*CONJG(WK(2,2)) TERM3 = WK(1,1)*CONJG(WK(2,2)) - CONJG(WK(2,1))*WK(1,2) TERM4 = WK(2,1)*CONJG(WK(2,2)) - CONJG(WK(2,1))*WK(2,2) WKR(2,2) = 1 / (REAL(TERM1*CONJG(TERM1) + 2*TERM2*CONJG(TERM2) 1 + TERM3*CONJG(TERM3) + TERM4*CONJG(TERM4))) IF (LEN .GE. 3) THEN DO 3 L = 3, LEN TERM1 = (WK(L,1)*CONJG(WK(L,2)) - CONJG(WK(L,1))*WK(L,2)) WKR(L,2) = WKR(L-1,2) / (1 + WKR(L-1,2)*REAL(TERM1*CONJG(TERM1))) DO 2 I = 1, L-1 TERM1 = WK(L,1)*CONJG(WK(I,2)) - CONJG(WK(I,1))*WK(L,2) TERM3 = WK(I,1)*CONJG(WK(L,2)) - CONJG(WK(L,1))*WK(I,2) TERM2 = CONJG(WK(L,1))*CONJG(WK(I,2)) - 1 CONJG(WK(I,1))*CONJG(WK(L,2)) WKR(L,2) = WKR(L,2) / (1 + WKR(L,2)*REAL(TERM1*CONJG(TERM1) + 1 2*TERM2*CONJG(TERM2) + TERM3*CONJG(TERM3))) 2 CONTINUE 3 CONTINUE ENDIF DO 6 L=1,LEN WKR(L,2) = WKR(L,2)*4 6 CONTINUE C C compute Christoffel function C WKR(1,1)=1/REAL(WK(1,1)*CONJG(WK(1,1))) DO 4 I=2,LEN WKR(I,1)=WKR(I-1,1)/(1+WKR(I-1,1)*REAL(WK(I,1)*CONJG(WK(I,1)))) 4 CONTINUE C C compute real inner product C SUM=0.0D0 DO 5 I=1,LEN SUM = SUM + REAL(WK(I,1)*CONJG(WK(I,2))) WKR(I,3) = SUM 5 CONTINUE GRVAL=B(1)*B(1)* 1 WKR(LEN,1)*(WKR(LEN,3)-(0.0D0,1.0D0)/SQRT(WKR(LEN,2))) RETURN END SUBROUTINE HH(A,MM,M,N,PT,EMACH) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(MM,N),PT(N) DO 3 IR=1,N S2=0.0 DO 4 I=IR,M 4 S2=S2+A(I,IR)**2 S=SQRT(S2) U=A(IR,IR)+S IF (ABS(ABS(U)-ABS(A(IR,IR))-S)-EMACH) 5,5,6 5 V=1.0 GOTO 7 6 V=-1.0 7 A(IR,IR)=A(IR,IR)+V*S SK=A(IR,IR)*S*V IF (IR-N) 16,3,3 16 IRP=IR+1 DO 9 J=IRP,N PT(J)=0.0 DO 8 I=IR,M 8 PT(J)=PT(J)+A(I,IR)*A(I,J) 9 PT(J)=PT(J)/SK DO 10 I=IR,M DO 10 J=IRP,N 10 A(I,J)=A(I,J)-A(I,IR)*PT(J) 3 PT(IR)=-S*V RETURN END SUBROUTINE INFGEN(ICN,N,WT,EPS,RNS,WTS,NPTS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION RNS(N),WTS(N) C C a rewrite of INITAL (q.v.) which does not need a common block C and has no internal limit on the number of quadrature points. C PI2=8*ATAN(1.0D0) J=1 RM=4 IC=1 M=4 3 MM=M-1 DO 4 I=1,MM,IC RI=I THETA=RI*PI2/RM RNS(J)=RI/RM-SIN(THETA)/PI2 WTS(J)=1-COS(THETA) J=J+1 4 CONTINUE M=2*M RM=2*RM GOTO(5,6),IC 5 IC=2 GOTO 3 6 IF(M .LE. N+1) GOTO 3 NPTS=MM CON=WT/(NPTS+1) IF (ICN .EQ. 6) THEN C tansform the half infinite range DO 10 I=1,NPTS T = RNS(I) RNS(I)=T/(1-T) WTS(I)=CON*WTS(I)*EXP(-RNS(I))/(1-T)**2 10 CONTINUE ELSE C transform the doubly infnite range DO 7 I=1,NPTS T=RNS(I)*2-1 RNS(I)=T/SQRT(1-T*T) WTS(I)=WTS(I)*SQRT((1+RNS(I)*RNS(I))**3) 7 CONTINUE IF (ICN .EQ. 4) THEN C quad points for gaussian (exp(-x**2)) CON = 2 * CON / SQRT(PI2/2) DO 8 I=1,NPTS WTS(I) = CON * WTS(I) * EXP(-RNS(I)*RNS(I)) 8 CONTINUE ELSE C quad points for SECH (Meixner) PIB2=PI2/4 CON = CON*2 DO 9 I=1,NPTS E=EXP(-ABS(PIB2*RNS(I))) WTS(I) = CON * WTS(I) * E / (1 + E*E) 9 CONTINUE ENDIF ENDIF ICOUNT=1 E20 = EPS**20 DO 11 I=1,NPTS IF(WTS(I) .GE. E20)THEN WTS(ICOUNT) = WTS(I) RNS(ICOUNT) = RNS(I) ICOUNT = ICOUNT + 1 ENDIF 11 CONTINUE NPTS = ICOUNT - 1 RETURN END SUBROUTINE INITAL(MMAX) IMPLICIT REAL*8 (A-H,O-Z) COMMON /BLKI/XV(511),WV(511) PI=1 PI2=8*ATAN(PI) J=1 RM=4 IC=1 M=4 2 MM=M-1 DO 3 I=1,MM,IC RI=I THETA=RI*PI2/RM XV(J)=RI/RM-SIN(THETA)/PI2 WV(J)=1-COS(THETA) J=J+1 3 CONTINUE M=2*M RM=2*RM GOTO(4,6),IC 4 IC=2 GOTO 2 6 IF(M .LE. MMAX+1) GOTO 2 8 CONTINUE RETURN END FUNCTION NUMC(A,B2,ALAM,LM1) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LM1),B2(LM1) 7 NU=0 P0=0 P1=1 DO 6 I=1,LM1 P2=(A(I)-ALAM)*P1-B2(I)*P0 SC=ABS(P1)+ABS(P2) IF ( P2*P1 .GT. 0 ) GOTO 4 IF ( P2*P1 .LT. 0 ) GOTO 5 14 IF ( P0*P2 .GE. 0 ) GOTO 5 4 NU=NU+1 5 P0=P1/SC P1=P2/SC 6 CONTINUE NUMC=NU RETURN END SUBROUTINE NUMD(A,B2,ALAM,LM1,DE) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LM1),B2(LM1) P0=0 PP0=0 P1=1 PP1=0 DO 6 I=1,LM1 P2=(A(I)-ALAM)*P1-B2(I)*P0 PP2=(A(I)-ALAM)*PP1-B2(I)*PP0-P1 SC=ABS(P1)+ABS(P2) P0=P1/SC PP0=PP1/SC PP1=PP2/SC P1=P2/SC 6 CONTINUE DE=P1/PP1 RETURN END SUBROUTINE PREPAS(LEN,LFIT,EMACH,EMAT,NE,EDIAG) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION EMAT(NE,2),EDIAG(NE) C prepare matrices for L2 fittings to the log of a power law decay C using HH ready for BS to complete the solution. LSTART=LEN-LFIT+1 DO 1 I=1,LFIT EMAT(I,1)= 1 EMAT(I,2)= LOG(DFLOAT(I+LEN-LFIT)) 1 CONTINUE CALL HH(EMAT,LFIT,LFIT,2,EDIAG,EMACH) RETURN END SUBROUTINE QFINT(A,B,F,FINT,ERR,E,MMAX,IP) IMPLICIT REAL*8 (A-H,O-Z) COMMON /BLKI/XV(511),WV(511) IF(IP.GT.0)WRITE(IP,1) 1 FORMAT(27H NO. OF FN. CALLS INTEGRAL,15X,5HERROR) IC=1 M=4 J=1 ESM=0 2 MM=M-1 SIGMAM=0 DO 3 I=1,MM,IC ETA=(B-A)*XV(J)+A SIGMAM=SIGMAM+WV(J)*F(ETA) J=J+1 3 CONTINUE SIGMAM=SIGMAM*2*(B-A)/M E=ABS(ESM-SIGMAM)/2 ESM=(ESM+SIGMAM)/2 IF(IP.GT.0)WRITE(IP,5)MM,ESM,E 5 FORMAT(I12,4X,3E20.12) M=2*M GOTO (4,6),IC 4 IC=2 GOTO 2 6 IF ( E .LE. ERR ) GOTO 8 7 IF ( M .LE. MMAX+1) GOTO 2 8 FINT=ESM RETURN END SUBROUTINE QVINT(ENDS,N,F,FINT,ERR,E,MMAX,IP,WK) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ENDS(N,2),WK(N) COMMON /BLKI/XV(511),WV(511) IF(IP.GT.0)WRITE(IP,1) 1 FORMAT(27H NO. OF FN. CALLS INTEGRAL,15X,5HERROR) IC=1 M=4 J=1 ESM=0 RLEN=0 DO 11 IV=1,N RLEN = RLEN + (ENDS(IV,2) - ENDS(IV,1))**2 11 CONTINUE RLEN = SQRT(RLEN) 2 MM=M-1 SIGMAM=0 DO 3 I=1,MM,IC DO 21 IV=1,N WK(IV)=(ENDS(IV,2)-ENDS(IV,1))*XV(J)+ENDS(IV,1) 21 CONTINUE SIGMAM=SIGMAM+WV(J)*F(WK,N) J=J+1 3 CONTINUE SIGMAM=SIGMAM*2*RLEN/M E=ABS(ESM-SIGMAM)/2 ESM=(ESM+SIGMAM)/2 IF(IP.GT.0)WRITE(IP,5)MM,ESM,E 5 FORMAT(I12,4X,3E20.12) M=2*M GOTO (4,6),IC 4 IC=2 GOTO 2 6 IF ( E .LE. ERR ) GOTO 8 7 IF ( M .LE. MMAX+1) GOTO 2 8 FINT=ESM RETURN END SUBROUTINE QVINT1(CORN,N,F,FINT,ERR,E,MMAX,IP,AVAC,WK) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CORN(N,3),WK(N,2) COMMON /BLKI/XV(511),WV(511) EXTERNAL F IF(IP.GT.0)WRITE(IP,1) 1 FORMAT(27H NO. OF FN. CALLS INTEGRAL,15X,5HERROR) IC=1 M=4 J=1 ESM=0 ERR1=ERR*0.1 IC1=0 R1=0 R2=0 R3=0 DO 11 IV=1,N R1 = R1 + (CORN(IV,3) - CORN(IV,1))**2 R2 = R2 + (CORN(IV,3) - CORN(IV,1))*(CORN(IV,2)-CORN(IV,1)) R3 = R3 + (CORN(IV,2) - CORN(IV,1))**2 11 CONTINUE RLEN = SQRT(R1 - R2*R2/R3) 2 MM=M-1 SIGMAM=0 DO 3 I=1,MM,IC DO 21 IV=1,N WK(IV,1)=(CORN(IV,3)-CORN(IV,1))*XV(J)+CORN(IV,1) 21 CONTINUE SIGMAM=SIGMAM+WV(J)*FVINT1(WK,F,CORN,N,ERR1,MMAX,ACC1,IC1) J=J+1 3 CONTINUE SIGMAM=SIGMAM*2*RLEN/M E=ABS(ESM-SIGMAM)/2 ESM=(ESM+SIGMAM)/2 IF(IP.GT.0)WRITE(IP,5)MM,ESM,E 5 FORMAT(I12,4X,3E20.12) M=2*M GOTO (4,6),IC 4 IC=2 GOTO 2 6 IF ( E .LE. ERR ) GOTO 8 7 IF ( M .LE. MMAX+1) GOTO 2 8 FINT=ESM AVAC=ACC1/IC1 RETURN END FUNCTION FVINT1(ARG,F,CORN,N,ERR,MMAX,ACC1,IC1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ARG(N,2),CORN(N,4) EXTERNAL F DO 1 IV=1,N ARG(IV,2) = ARG(IV,1) + CORN(IV,2) - CORN(IV,1) 1 CONTINUE CALL QVINT(ARG,N,F,ANS,ERR,ACC,MMAX,0,WK) FVINT1=ANS ACC1 = ACC1 + ACC IC1=IC1+1 RETURN END SUBROUTINE QVINT2(CORN,N,F,FINT,ERR,E,MMAX,IP,AVAC,WK) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CORN(N,4),WK(N,3),AVAC(2) COMMON /BLKI/XV(511),WV(511) EXTERNAL F IF(IP.GT.0)WRITE(IP,1) 1 FORMAT(27H NO. OF FN. CALLS INTEGRAL,15X,5HERROR) IC=1 M=4 J=1 ESM=0 ERR1=ERR*0.1 ICNT=0 R1=0 R2=0 R3=0 R12=0 R13=0 R23=0 DO 11 IV=1,N R1 = R1 + (CORN(IV,2)-CORN(IV,1))**2 R2 = R2 + (CORN(IV,3)-CORN(IV,1))**2 R3 = R3 + (CORN(IV,4)-CORN(IV,1))**2 R12 = R12 + (CORN(IV,2)-CORN(IV,1))*(CORN(IV,3)-CORN(IV,1)) R23 = R23 + (CORN(IV,3)-CORN(IV,1))*(CORN(IV,4)-CORN(IV,1)) R13 = R13 + (CORN(IV,2)-CORN(IV,1))*(CORN(IV,4)-CORN(IV,1)) 11 CONTINUE RLEN=SQRT((R1*R2*R3-R1*R23*R23-R2*R13*R13-R3*R12*R12+ 1 2*R12*R23*R13)/(R1*R2-R12**2)) 2 MM=M-1 SIGMAM=0 DO 3 I=1,MM,IC DO 21 IV=1,N WK(IV,1)=(CORN(IV,4)-CORN(IV,1))*XV(J)+CORN(IV,1) 21 CONTINUE SIGMAM=SIGMAM+WV(J)*FVINT2(WK,F,CORN,N,ERR1,MMAX,AVAC,ICNT) J=J+1 3 CONTINUE SIGMAM=SIGMAM*2*RLEN/M E=ABS(ESM-SIGMAM)/2 ESM=(ESM+SIGMAM)/2 IF(IP.GT.0)WRITE(IP,5)MM,ESM,E 5 FORMAT(I12,4X,3E20.12) M=2*M GOTO (4,6),IC 4 IC=2 GOTO 2 6 IF ( E .LE. ERR ) GOTO 8 7 IF ( M .LE. MMAX+1) GOTO 2 8 FINT=ESM DO 9 I=1,2 AVAC(I)=AVAC(I)/ICNT 9 CONTINUE RETURN END FUNCTION FVINT2(ARG,F,CORN,N,ERR,MMAX,VACC,ICNT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ARG(N,3),CORN(N,4),VACC(2) EXTERNAL F DO 1 IV=1,N DO 1 ID=2,3 ARG(IV,ID) = ARG(IV,1) + CORN(IV,ID) - CORN(IV,1) 1 CONTINUE DO 2 I=1,2 VACC(I)=0 2 CONTINUE CALL QVINT1(ARG,N,F,ANS,ERR,ACC,MMAX,0,ACC1,WK) FVINT2=ANS VACC(1) = VACC(1) + ACC1 VACC(2) = VACC(2) + ACC ICNT=ICNT+1 RETURN END SUBROUTINE RECQD(A,B2,LM1,X,WT,M,EPS,WK,MU) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LM1),B2(LM1),X(LM1),WT(LM1),WK(LM1),MU(LM1) DIMENSION P(2,3) M=1 CALL RECRTS(A,B2,LM1,EPS,XLIM,M,X,MU,WK,WT) DO 7 II=1,M WT(II)=RECWT(X(II),A,B2,LM1,EPS,0,P,1) 7 CONTINUE RETURN END SUBROUTINE RECRTS(A,B2,LM1,EPS,XLIM,N,X,MULT,BI,NI) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LM1),B2(LM1),X(LM1),MULT(LM1),BI(LM1),NI(LM1) IF(LM1.GT.1)GOTO 1 IF(N.EQ.0.AND.A(1).GT.XLIM)RETURN X(1)=A(1) MULT(1)=1 N=1 RETURN 1 MULT(1)=0 C C ESTIMATE LIMITS OF ROOTS C CN=MIN(A(1)-B2(2),A(LM1)-1) NID=LM1-1 IF(NID.LT.2)GOTO 20 DO 2 I=2,NID CN=MIN(CN,A(I)-1-B2(I+1)) 2 CONTINUE 20 AA=CN IF(N.EQ.0)GOTO 4 CN=MAX(A(1)+B2(2),A(LM1)+1) IF(NID.LT.2)GOTO 21 DO 3 I=2,NID CN=MAX(CN,A(I)+1+B2(I+1)) 3 CONTINUE 21 B=CN 4 IF(N.EQ.0) B=XLIM NMAX=MAX(LOG(EPS/(B-AA))/LOG(0.5D0)*10,70.0D0) NA=NUMC(A,B2,AA,LM1) NB=NUMC(A,B2,B,LM1) 5 NL=NA NM=NB N=0 IF (NL.LE.NM) GOTO 19 NE=NL-NM DO 22 I=1,NE X(I)=0 MULT(I)=0 22 CONTINUE I=1 NI(1)=NB BI(1)=B NCOUNT=0 C C START BISECTION SEARCH C 9 C=(AA+B)/2 NC=NUMC(A,B2,C,LM1) NCOUNT=NCOUNT+1 IF(NCOUNT.GE.NMAX)GOTO 16 IF(NC.LT.NL) GOTO 11 AA=C NA=NC GOTO 15 11 IF(NC.LE.NM)GOTO 14 IF(NI(I).EQ.NC) GOTO 14 I=I+1 BI(I)=C NI(I)=NC 14 B=C NB=NC 15 IF( (NA-NB).GT.1 ) GOTO 36 C DO A FEW NEWTON-RAPHSON STEPS C=(AA+B)/2 DO 34 ITR=1,10 CALL NUMD(A,B2,C,LM1,DE) C=C-DE IF(C.GT.B .OR. C.LT.AA) GOTO 9 IF( ABS(DE).LE.EPS) GOTO 35 34 CONTINUE GOTO 9 35 XX=MAX(C-EPS,AA) NXX=NUMC(A,B2,XX,LM1) IF( NXX.NE.NA) GOTO 9 XX=MIN(C+EPS,B) NXX=NUMC(A,B2,XX,LM1) IF( NXX.NE.NB) GOTO 9 N=N+1 X(N)=C MULT(N)=1 GOTO 37 36 IF(ABS(AA-B).GT.EPS)GOTO 9 16 N=N+1 X(N)=(AA+B)/2 MULT(N)=NA-NB IF(NCOUNT.GE.NMAX)MULT(N)=-MULT(N) 37 IF(NB.LE.NM)RETURN IF(I.LE.1)GOTO 19 18 AA=BI(I) NA=NI(I) I=I-1 B=BI(I) NB=NI(I) NL=NL-ABS(MULT(N)) NCOUNT=0 GOTO 9 19 MULT(1)=-MULT(1) RETURN END SUBROUTINE RECSTD(A,B,IC,WT,N,EPS,RNS,WTS,WK,NWK) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION RNS(N),WTS(N),WK(NWK,4) IF (IC .EQ. 0)THEN WTS(1)=WT RNS(1)=A N=1 RETURN ENDIF PI=1 PI=4*ATAN(PI) IF (IC .GT. 1) GOTO 1 C constant weight: Legendre polynomial weight C generate recurrence coefficients for Legendre polys DO 11 L=1,N WK(L,1)=A AL2=(L-1)*(L-1) WK(L,2)=B*B*AL2/(4*AL2-1) 11 CONTINUE WK(1,2)=WT GOTO 50 1 IF (IC .GT. 2) GOTO 2 C 1/squareroot(1-x*x): Chebyshev weight of the first kind TH=PI/(2*N) SC=WT/N DO 21 I = 1, N WTS(I)=SC RNS(I) = A-B*COS((2*I-1)*TH) 21 CONTINUE RETURN 2 IF (IC .GT. 3) GOTO 3 C squareroot(1-x*x): Chebyshev weight of the second kind TH=PI/(N+1) CON = 2*WT/(N+1) DO 31 I=1,N RNS(I)=COS(TH*(N-I+1)) WTS(I)=CON*(1-RNS(I)*RNS(I)) RNS(I)=A+B*RNS(I) 31 CONTINUE RETURN C gaussian & Meixner and Laguerre: nodes generated by quad. method in QFINT 3 IF (IC .GT. 6) GOTO 4 CALL INFGEN(IC,N,WT,EPS,RNS,WTS,NACT) N=NACT DO 41 I=1,N RNS(I)=RNS(I)*B+A 41 CONTINUE RETURN 4 IF (IC .GT. 7) GOTO 5 C Pollaczek weight. Cfts from Chihara - CHA is Chihara's A, and CHB is C his B. Lambda is 1 and C is zero. CHB=A CHA=B DO 51 L=1,N WK(L,1) = -CHB/(L+CHA) WK(L,2) = L*(L-1)/(4*(L+CHA)*(L+CHA-1)) 51 CONTINUE WK(1,2)=WT GOTO 50 5 CHA=B C Charlier weights. Cfts from Chihara - CHA is Chihara's A DO 61 L=1,N WK(L,1) = L-1+CHA WK(L,2) = CHA*(L-1) 61 CONTINUE WK(1,2)=WT 50 CONTINUE C Generate nodes from recurrence CALL RECQD(WK(1,1),WK(1,2),N,RNS,WTS,LO,EPS,WK(1,3), 1 WK(1,4)) IF ( LO .NE. N) THEN IC = -1 RETURN ENDIF RETURN END FUNCTION RECWT(E,A,B2,LL,EPS,N,P,NS) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(LL),B2(LL),P(2,3) RECWT=B2(1) IF(LL.EQ.1)RETURN K=NS IF(K.GE.2)GOTO 1 P(1,1)=1 P(2,1)=E-A(1) P(1,2)=0 P(2,2)=1 P(1,3)=0 P(2,3)=B2(1) K=2 1 LLIM=LL-ABS(N) IF(LLIM.LT.K)GOTO 4 DO 3 L=K,LLIM SC=ABS(P(1,1))+ABS(P(1,3)) DO 2 J=1,3 DUM=P(2,J) P(2,J)=((E-A(L))*DUM-B2(L)*P(1,J))/SC P(1,J)=DUM/SC 2 CONTINUE P(2,2)=P(2,2)+P(1,1) 3 CONTINUE 4 IF ( N .NE. 0 ) GOTO 6 5 RECWT=P(2,3)/P(2,2) GOTO 7 6 RECWT=(P(1,1)*P(2,3)-P(2,1)*P(1,3))/(P(1,1)*P(2,2)-P(1,2) 1 *P(2,1)+P(2,1)*P(2,1)/B2(LL)) 7 IF(RECWT.LT.0)RECWT=0 IF(N.GE.0)RETURN IF(ABS(P(2,1)).LE.EPS*ABS(P(1,1)))RETURN N=0 A(LL)=E-B2(LL)*P(1,1)/P(2,1) RETURN END SUBROUTINE SIMBIS(FUN,A,B,EPS,NIT,X) IMPLICIT REAL*8(A-H,O-Z) FL=FUN(A) FR=FUN(B) XL=A XR=B IF( FL*FR .GT. 0)THEN NIT=0 RETURN ENDIF NIT = LOG((B-A)/EPS)/LOG(2.0)+1 DO 1 IT=1,NIT X=(XL+XR)/2 FX=FUN(X) IF( FX*FL .LT. 0) THEN XR=X FR=FX ELSE XL=X FL=FX ENDIF 1 CONTINUE X=(XL+XR)/2 RETURN END SUBROUTINE VECCF(Z,AV,BETA,NV,NVE,LEN,TERM,VAL) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Z(NVE),AV(NV,LEN),BETA(LEN),TERM(NVE),VAL(NVE) DO 1 I=1,NVE VAL(I)=TERM(I) 1 CONTINUE DO 4 LEV=LEN,1,-1 SUM=0 DO 2 I=1,NV VAL(I) = Z(I)-AV(I,LEV)-VAL(I) SUM = SUM+VAL(I)*VAL(I) 2 CONTINUE VAL(NVE) = -(Z(NVE)-VAL(NVE)) SUM = SUM + VAL(NVE)*VAL(NVE) C CALL CONJG(VAL,NVE) DO 3 I=1,NVE VAL(I) = VAL(I)*BETA(LEV)/SUM 3 CONTINUE 4 CONTINUE RETURN END SUBROUTINE VECPOL(AV,BETA,NV,LEN,ARG,PV,PSCAL,CFTL) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AV(NV,LEN),BETA(LEN),ARG(NV), 1 PV(NV,LEN),PSCAL(LEN),CFTL(LEN) SUM1=0 DO 1 I=1,NV PV(I,1)=0 PV(I,2)=(ARG(I)-AV(I,1))/(BETA(1)*BETA(2)) SUM1 = SUM1 + (ARG(I)-AV(I,1))*(ARG(I)-AV(I,1)) 1 CONTINUE PSCAL(1)=1/BETA(1) PSCAL(2)=SUM1*PSCAL(1)/BETA(2) CFTL(1)=1/PSCAL(1) CFTL(2)=1/(PSCAL(1)+PSCAL(2)) DO 3 LEV=3,LEN SUM1=0 SUM2=0 DO 2 I=1,NV PV(I,LEV) = ((ARG(I)-AV(I,LEV-1))*PSCAL(LEV-1) - 1 BETA(LEV-1)*PV(I,LEV-1)) / BETA(LEV) SUM1 = SUM1 + (ARG(I)-AV(I,LEV-1))*(ARG(I)-AV(I,LEV-1)) SUM2 = SUM2 + (ARG(I)-AV(I,LEV-1))*PV(I,LEV-1) 2 CONTINUE PSCAL(LEV) = ( SUM1*PSCAL(LEV-1) - 2*BETA(LEV-1)*SUM2 + 1 BETA(LEV-1)*PSCAL(LEV-2) ) / BETA(LEV) CFTL(LEV) = CFTL(LEV-1) / (1 + CFTL(LEV-1)*PSCAL(LEV)) 3 CONTINUE RETURN END