C ALGORITHM 793, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 25,NO. 2, June, 1999, P. 213--239. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc/ # Doc/Makefile_dp # Doc/Makefile_sp # Doc/PackingList # Fortran77/ # Fortran77/Drivers/ # Fortran77/Drivers/Dp/ # Fortran77/Drivers/Dp/driver1.f # Fortran77/Drivers/Dp/driver2.f # Fortran77/Drivers/Dp/driver3.f # Fortran77/Drivers/Dp/driver5.f # Fortran77/Drivers/Dp/driver6.f # Fortran77/Drivers/Dp/driver7.f # Fortran77/Drivers/Dp/res1 # Fortran77/Drivers/Dp/res2 # Fortran77/Drivers/Dp/res3 # Fortran77/Drivers/Dp/res5 # Fortran77/Drivers/Dp/res6 # Fortran77/Drivers/Dp/res7 # Fortran77/Drivers/Sp/ # Fortran77/Drivers/Sp/driver1.f # Fortran77/Drivers/Sp/driver2.f # Fortran77/Drivers/Sp/driver3.f # Fortran77/Drivers/Sp/driver4.f # Fortran77/Drivers/Sp/driver5.f # Fortran77/Drivers/Sp/driver6.f # Fortran77/Drivers/Sp/driver7.f # Fortran77/Drivers/Sp/res1 # Fortran77/Drivers/Sp/res2 # Fortran77/Drivers/Sp/res3 # Fortran77/Drivers/Sp/res4 # Fortran77/Drivers/Sp/res5 # Fortran77/Drivers/Sp/res6 # Fortran77/Drivers/Sp/res7 # Fortran77/Src/ # Fortran77/Src/Dp/ # Fortran77/Src/Dp/src.f # Fortran77/Src/Sp/ # Fortran77/Src/Sp/src.f # This archive created: Fri Oct 22 09:59:36 1999 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'Makefile_dp' then echo shar: will not over-write existing file "'Makefile_dp'" else cat << SHAR_EOF > 'Makefile_dp' ######################################################################## # # # Makefile for GQRAT # # # # Usage: # # # # Type make to compile and run all test programs # # Type make resN to compile and run the test program # # number N # # Type make clean to remove all but the generated output # # # ######################################################################## #---------------------------------------------------------------------- # This section specifies the name of the linker and compilation options # for various systems. All but one are commented out. #---------------------------------------------------------------------- # ... For the Sun LD = f77 FFLAGS = # ... For the Cray #LD = cf77 #FFLAGS = # ... For the Convex #LD = fc #FFLAGS = -O2 # ... For the RS/6000 #LD = xlf #FFLAGS = -O # ... For the HP 9000 series #LD = f77 #FFLAGS = -O -K -w # ... For Silicon Graphics #LD = f77 #FFLAGS = -O2 #------------------------------------------------------------------- # This section specifies which object files are needed for each test # program. #------------------------------------------------------------------- TEST1 = driver1.o d1mach.o src.f TEST2 = driver2.o d1mach.o src.f TEST3 = driver3.o d1mach.o src.f TEST5 = driver5.o d1mach.o src.f TEST6 = driver6.o d1mach.o src.f TEST7 = driver7.o d1mach.o src.f #-------------------------------------------------------------- # This is the default target. It will cause all the tests to be # generated and run, saving their outputs in files. #-------------------------------------------------------------- all : res1 res2 res3 res5 res6 res7 #---------------------------------------------------------------------- # This section shows how to run each test program. The output for driverN # is put in resN. #---------------------------------------------------------------------- res1 : test1 ; test1 > $@ res2 : test2 ; test2 > $@ res3 : test3 ; test3 > $@ res5 : test5 ; test5 > $@ res6 : test6 ; test6 > $@ res7 : test7 ; test7 > $@ #-------------------------------------------------- # This section shows how to link each test program. #-------------------------------------------------- test1 : $(TEST1) ; $(LD) -o $@ $(TEST1) test2 : $(TEST2) ; $(LD) -o $@ $(TEST2) test3 : $(TEST3) ; $(LD) -o $@ $(TEST3) test5 : $(TEST5) ; $(LD) -o $@ $(TEST5) test6 : $(TEST6) ; $(LD) -o $@ $(TEST6) test7 : $(TEST7) ; $(LD) -o $@ $(TEST7) #---------------------------------------------- # This section shows how to clean up afterward. #---------------------------------------------- clean : rm -f core test1 test2 test3 test5 test6 test7 SHAR_EOF fi # end of overwriting check if test -f 'Makefile_sp' then echo shar: will not over-write existing file "'Makefile_sp'" else cat << SHAR_EOF > 'Makefile_sp' ######################################################################## # # # Makefile for GQRAT # # # # Usage: # # # # Type make to compile and run all test programs # # Type make resN to compile and run the test program # # number N # # Type make clean to remove all but the generated output # # # ######################################################################## #---------------------------------------------------------------------- # This section specifies the name of the linker and compilation options # for various systems. All but one are commented out. #---------------------------------------------------------------------- # ... For the Sun LD = f77 FFLAGS = # ... For the Cray #LD = cf77 #FFLAGS = # ... For the Convex #LD = fc #FFLAGS = -O2 # ... For the RS/6000 #LD = xlf #FFLAGS = -O # ... For the HP 9000 series #LD = f77 #FFLAGS = -O -K -w # ... For Silicon Graphics #LD = f77 #FFLAGS = -O2 #------------------------------------------------------------------- # This section specifies which object files are needed for each test # program. #------------------------------------------------------------------- TEST1 = driver1.o r1mach.o src.f TEST2 = driver2.o r1mach.o src.f TEST3 = driver3.o r1mach.o src.f TEST4 = driver4.o r1mach.o src.f TEST5 = driver5.o r1mach.o src.f TEST6 = driver6.o r1mach.o src.f TEST7 = driver7.o r1mach.o src.f #-------------------------------------------------------------- # This is the default target. It will cause all the tests to be # generated and run, saving their outputs in files. #-------------------------------------------------------------- all : res1 res2 res3 res4 res5 res6 res7 #---------------------------------------------------------------------- # This section shows how to run each test program. The output for driverN # is put in resN #---------------------------------------------------------------------- res1 : test1 ; test1 > $@ res2 : test2 ; test2 > $@ res3 : test3 ; test3 > $@ res4 : test4 ; test4 > $@ res5 : test5 ; test5 > $@ res6 : test6 ; test6 > $@ res7 : test7 ; test7 > $@ #-------------------------------------------------- # This section shows how to link each test program. #-------------------------------------------------- test1 : $(TEST1) ; $(LD) -o $@ $(TEST1) test2 : $(TEST2) ; $(LD) -o $@ $(TEST2) test3 : $(TEST3) ; $(LD) -o $@ $(TEST3) test4 : $(TEST4) ; $(LD) -o $@ $(TEST4) test5 : $(TEST5) ; $(LD) -o $@ $(TEST5) test6 : $(TEST6) ; $(LD) -o $@ $(TEST6) test7 : $(TEST7) ; $(LD) -o $@ $(TEST7) #---------------------------------------------- # This section shows how to clean up afterward. #---------------------------------------------- clean : rm -f core test1 test2 test3 test4 test5 test6 test7 SHAR_EOF fi # end of overwriting check if test -f 'PackingList' then echo shar: will not over-write existing file "'PackingList'" else cat << SHAR_EOF > 'PackingList' Algorithm:793 Language:Fortran77 Precision:Dp Src:src.f SrcLibs:port Driver_0:driver1.f @Src Res_0:stdout=res1 Driver_1:driver2.f @Src Res_1:stdout=res2 Driver_2:driver3.f @Src Res_2:stdout=res3 Driver_3:driver5.f @Src Res_3:stdout=res5 Driver_4:driver6.f @Src Res_4:stdout=res6 Driver_5:driver7.f @Src Res_5:stdout=res7 Precision:Sp Src:src.f SrcLibs:port Driver_6:driver1.f @Src Res_6:stdout=res1 Driver_7:driver2.f @Src Res_7:stdout=res2 Driver_8:driver3.f @Src Res_8:stdout=res3 Driver_9:driver4.f @Src Res_9:stdout=res4 Driver_10:driver5.f @Src Res_10:stdout=res5 Driver_11:driver6.f @Src Res_11:stdout=res6 Driver_12:driver7.f @Src Res_12:stdout=res7 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'driver1.f' then echo shar: will not over-write existing file "'driver1.f'" else cat << SHAR_EOF > 'driver1.f' c c PROGRAM DTEST1 c c This is a double-precision version of test1 c C .. Parameters .. INTEGER NMAX,NCAPM PARAMETER (NMAX=10,NCAPM=100) DOUBLE PRECISION OM PARAMETER (OM=1.1D0) C .. C .. Local Scalars .. DOUBLE PRECISION CONST,EPS,EXACT,SGN,SUM REAL ERR INTEGER IERAB,IERGQ,IERR,IROUT,K,KOUNT,M,MU,N,NCAP C .. C .. Local Arrays .. DOUBLE PRECISION A(NCAPM),ALPHA(NMAX+1),B(NCAPM),BE(NMAX+1), + BETA(NMAX+1),E(NCAPM),P0(NCAPM),P1(NCAPM), + P2(NCAPM),W(NCAPM),WG(NMAX),X(NCAPM),XII(2*NMAX), + XIR(2*NMAX),ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. DOUBLE PRECISION D1MACH,F EXTERNAL D1MACH,F C .. C .. External Subroutines .. EXTERNAL DABMOD,DGQRAT,DRECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,REAL integer i1mach, nout nout = i1mach(2) C .. IROUT = 1 c irout=0 EPS = D1MACH(3)*1.D2 CALL DRECUR(NCAPM,1,0.D0,0.D0,A,B,IERR) WRITE(NOUT,FMT=9000) DO 30 N = 1,NMAX M = 2*N c m=n c m=2 c m=0 SGN = 1.D0 DO 10 MU = 1,M SGN = -SGN XIR(MU) = SGN/ (OM*DBLE((MU+1)/2)) XII(MU) = 0.D0 IS(MU) = 1 10 CONTINUE CALL DABMOD(N+1,NCAPM,M,EPS,IROUT,A,B,XIR,XII,IS,ALPHA,BETA, + NCAP,KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9010) IERAB CALL DGQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9020) IERGQ SUM = 0.D0 DO 20 K = 1,N SUM = SUM + WG(K)*F(ZG(K),OM) 20 CONTINUE EXACT = 4.467773646387766D0 ERR = ABS(REAL((SUM-EXACT)/EXACT)) WRITE(NOUT,FMT=9030) N,M,SUM,ERR,CONST 30 CONTINUE STOP 9000 FORMAT (5X,'n',4X,'m',8X,'integral',9X,'rel error',3X,'err const', + /) 9010 FORMAT (1X,'ierab=',I3) 9020 FORMAT (1X,'iergq=',I3) 9030 FORMAT (1X,2I5,D23.14,E12.4,D12.4) END DOUBLE PRECISION FUNCTION F(T,OM) C .. Scalar Arguments .. DOUBLE PRECISION OM,T C .. C .. Local Scalars .. DOUBLE PRECISION PI C .. C .. Intrinsic Functions .. INTRINSIC ATAN,SIN C .. PI = 4.D0*ATAN(1.D0) F = 1.D0 IF (T.NE.0.D0) F = (PI*T/OM)/SIN(PI*T/OM) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'driver2.f' then echo shar: will not over-write existing file "'driver2.f'" else cat << SHAR_EOF > 'driver2.f' c c PROGRAM DTEST2 c c This is a double-precision version of test2 c C .. Parameters .. INTEGER NMAX,NCAPM PARAMETER (NMAX=10,NCAPM=100) DOUBLE PRECISION OM PARAMETER (OM=1.001D0) C .. C .. Local Scalars .. DOUBLE PRECISION CONST,EPS,EXACT,HN,HP,P,SGN,SUM,SUMN,SUMP,X1 REAL ERR INTEGER IERAB,IERGC,IERGQ,IERR,IROUT,K,KOUNT,M,MU,N,NCAP,NU C .. C .. Local Arrays .. DOUBLE PRECISION A(NCAPM),A1(NMAX+1),ALPHA(NMAX+1),B(NCAPM), + B1(NMAX+1),BE(NMAX+1),BETA(NMAX+1),E(NCAPM), + P0(NCAPM),P1(NCAPM),P2(NCAPM),W(NCAPM),WG(NMAX), + X(NCAPM),XII(2*NMAX),XIR(2*NMAX),ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. DOUBLE PRECISION D1MACH,F EXTERNAL D1MACH,F C .. C .. External Subroutines .. EXTERNAL DABMOD,DGCHRS,DGQRAT,DRECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,LOG,REAL integer i1mach, nout nout = i1mach(2) C .. IROUT = 1 c irout=0 EPS = D1MACH(3)*1.D2 CALL DRECUR(NCAPM,1,0.D0,0.D0,A,B,IERR) WRITE(NOUT,FMT=9000) DO 50 N = 1,NMAX M = 2*N c m=n c m=4 c m=2 SGN = 1.D0 DO 10 NU = 1,M SGN = -SGN IF (NU.LE.M-2) XIR(NU) = SGN/DBLE((NU+3)/2) XII(NU) = 0.D0 IS(NU) = 1 10 CONTINUE CALL DABMOD(N+1,NCAPM,M-2,EPS,IROUT,A,B,XIR,XII,IS,A1,B1,NCAP, + KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9010) IERAB XIR(M-1) = 1.D0/OM XIR(M) = -1.D0/OM X1 = OM HP = LOG(ABS((X1+1.D0)/ (X1-1.D0))) HN = -HP IF (M.GT.2) THEN SUMP = 0.D0 SUMN = 0.D0 DO 30 NU = 1,M - 2 P = 1.D0 DO 20 MU = 1,M - 2 IF (MU.EQ.NU) GO TO 20 P = (1.D0-XIR(MU)/XIR(NU))*P 20 CONTINUE SUMP = SUMP + LOG(ABS((X1+1.D0)* (XIR(NU)+1.D0)/ ((X1- + 1.D0)* (XIR(NU)-1.D0))))/ ((1.D0+XIR(NU)*X1)*P) SUMN = SUMN + LOG(ABS((X1-1.D0)* (XIR(NU)+1.D0)/ ((X1+ + 1.D0)* (XIR(NU)-1.D0))))/ ((1.D0-XIR(NU)*X1)*P) 30 CONTINUE HP = SUMP HN = SUMN END IF CALL DGCHRS(N+1,2,A1,B1,X1,HP,HN,ALPHA,BETA,IERGC) BETA(1) = - (OM**2)*BETA(1) CALL DGQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9020) IERGQ SUM = 0.D0 DO 40 K = 1,N SUM = SUM + WG(K)*F(ZG(K),OM) 40 CONTINUE EXACT = 12.929256850003D0 ERR = ABS(REAL((SUM-EXACT)/EXACT)) WRITE(NOUT,FMT=9030) N,M,SUM,ERR,CONST 50 CONTINUE STOP 9000 FORMAT (5X,'n',4X,'m',8X,'integral',10X,'rel error',4X, + 'err const',/) 9010 FORMAT (1X,'ierab=',I3) 9020 FORMAT (1X,'iergq=',I3) 9030 FORMAT (1X,2I5,D23.14,E13.4,D13.4) END DOUBLE PRECISION FUNCTION F(T,OM) C .. Scalar Arguments .. DOUBLE PRECISION OM,T C .. C .. Local Scalars .. DOUBLE PRECISION PI C .. C .. Intrinsic Functions .. INTRINSIC ATAN,SIN C .. PI = 4.D0*ATAN(1.D0) F = 1.D0 IF (T.NE.0.D0) F = (PI*T/OM)/SIN(PI*T/OM) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'driver3.f' then echo shar: will not over-write existing file "'driver3.f'" else cat << SHAR_EOF > 'driver3.f' c c PROGRAM DTEST3 c C .. Parameters .. INTEGER NMAX,NCAPM PARAMETER (NMAX=10,NCAPM=100) C .. C .. Local Scalars .. DOUBLE PRECISION CONST,EPS,OM,SUM REAL ERR INTEGER IERAB,IERGQ,IERR,IOM,IROUT,K,KOUNT,M,MU,N,NCAP C .. C .. Local Arrays .. DOUBLE PRECISION A(NCAPM),ALPHA(NMAX+1),B(NCAPM),BE(NMAX+1), + BETA(NMAX+1),E(NCAPM),EXACT(9),OOM(9),P0(NCAPM), + P1(NCAPM),P2(NCAPM),W(NCAPM),WG(NMAX),X(NCAPM), + XII(2*NMAX),XIR(2*NMAX),ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. DOUBLE PRECISION D1MACH,F EXTERNAL D1MACH,F C .. C .. External Subroutines .. EXTERNAL DABMOD,DGQRAT,DRECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,REAL,SQRT integer i1mach, nout C .. C .. Data statements .. DATA OOM/.1D0,.5D0,1.D0,1.999D0,2.D0,3.D0,3.001D0,10.5D0,50.D0/ DATA EXACT/7.6724625863706D0,2.5531371574419D0,1.4784672063106D0, + .81771417926000D0,.81735216597811D0,.56726444593607D0, + .56709144671961D0,.17313056604184D0,.037222082318054D0/ C .. nout = i1mach(2) WRITE(NOUT,FMT=9000) IROUT = 1 c irout=0 EPS = D1MACH(3)*1.D2 CALL DRECUR(NCAPM,6,0.D0,-.5D0,A,B,IERR) DO 10 K = 1,NCAPM A(K) = .5D0* (1.D0+A(K)) B(K) = .25D0*B(K) 10 CONTINUE B(1) = SQRT(8.D0)*B(1) WRITE(NOUT,FMT=9010) DO 50 IOM = 1,9 OM = OOM(IOM) DO 40 N = 1,NMAX M = 2*N c m=n c m=2 c m=0 IF (M.GT.0) THEN DO 20 MU = 1,M IF (MU.EQ.1) THEN XIR(MU) = 1.D0/OM ELSE XIR(MU) = 1.D0/DBLE(MU-1) END IF XII(MU) = 0.D0 IS(MU) = 1 20 CONTINUE END IF CALL DABMOD(N+1,NCAPM,M,EPS,IROUT,A,B,XIR,XII,IS,ALPHA, + BETA,NCAP,KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9020) IERAB CALL DGQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9030) IERGQ SUM = 0.D0 DO 30 K = 1,N SUM = SUM + WG(K)*F(ZG(K),OM) 30 CONTINUE ERR = ABS(REAL((SUM-EXACT(IOM))/EXACT(IOM))) IF (N.EQ.1) THEN WRITE(NOUT,FMT=9040) N,M,SUM,ERR,CONST,REAL(OM) ELSE WRITE(NOUT,FMT=9050) N,M,SUM,ERR,CONST END IF 40 CONTINUE WRITE(NOUT,FMT=9000) 50 CONTINUE STOP 9000 FORMAT (/) 9010 FORMAT (5X,'n',4X,'m',8X,'integral',10X,'rel error',4X, + 'err const',/) 9020 FORMAT (1X,'ierab=',I3) 9030 FORMAT (1X,'iergq=',I3) 9040 FORMAT (1X,2I5,D23.14,E13.4,D13.4,' om=',F6.3) 9050 FORMAT (1X,2I5,D23.14,E13.4,D13.4) END DOUBLE PRECISION FUNCTION F(T,OM) C .. Scalar Arguments .. DOUBLE PRECISION OM,T C .. C .. Local Scalars .. INTEGER IERR C .. C .. External Functions .. DOUBLE PRECISION DGAMMA EXTERNAL DGAMMA C .. F = DGAMMA(1.D0+T,IERR)/ (T+OM) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'driver5.f' then echo shar: will not over-write existing file "'driver5.f'" else cat << SHAR_EOF > 'driver5.f' c c PROGRAM DTEST5 c c This is a double-precision version of test5 c c parameter(eta=1.d0) c parameter(ak=1.5d0) c parameter(ak=2.5d0) C .. Parameters .. INTEGER NMAX,NCAPM PARAMETER (NMAX=10,NCAPM=200) DOUBLE PRECISION THETA PARAMETER (THETA=1.D-4) DOUBLE PRECISION ETA PARAMETER (ETA=-1.D0) DOUBLE PRECISION AK PARAMETER (AK=.5D0) C .. C .. Local Scalars .. DOUBLE PRECISION CONST,DEN,EPS,EXACT,PI,SUM REAL ERR INTEGER IERAB,IERGQ,IERR,IROUT,K,KOUNT,M,MU,N,NCAP C .. C .. Local Arrays .. DOUBLE PRECISION A(NCAPM),ALPHA(NMAX+1),B(NCAPM),BE(NMAX+1), + BETA(NMAX+1),E(NCAPM),P0(NCAPM),P1(NCAPM), + P2(NCAPM),W(NCAPM),WG(NMAX),X(NCAPM),XII(2*NMAX), + XIR(2*NMAX),ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. DOUBLE PRECISION D1MACH,F EXTERNAL D1MACH,F C .. C .. External Subroutines .. EXTERNAL DABMOD,DGQRAT,DRECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,DBLE,REAL integer i1mach, nout nout = i1mach(2) C .. IROUT = 1 c irout=0 EPS = D1MACH(3)*1.D2 CALL DRECUR(NCAPM,7,AK,0.D0,A,B,IERR) PI = 4.D0*ATAN(1.D0) WRITE(NOUT,FMT=9000) DO 30 N = 1,NMAX M = 2*N c m=2*((n+1)/2) c m=2 c m=0 IF (M.GT.0) THEN DO 10 MU = 1,M - 1,2 DEN = ETA**2 + DBLE(MU**2)* (PI**2) XIR(MU) = -ETA/DEN XIR(MU+1) = XIR(MU) XII(MU) = DBLE(MU)*PI/DEN XII(MU+1) = -XII(MU) IS(MU) = 1 IS(MU+1) = 1 10 CONTINUE END IF CALL DABMOD(N+1,NCAPM,M,EPS,IROUT,A,B,XIR,XII,IS,ALPHA,BETA, + NCAP,KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9010) IERAB CALL DGQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9020) IERGQ SUM = 0.D0 DO 20 K = 1,N SUM = SUM + WG(K)*F(ZG(K),ETA,THETA) 20 CONTINUE EXACT = .29051241701949D0 c exact=.46087845417799d0 c exact=1.18607350107576d0 c exact=1.39644182034912d0 c exact=2.66187327910715d0 c exact=7.62725609565345d0 ERR = ABS(REAL((SUM-EXACT)/EXACT)) WRITE(NOUT,FMT=9030) N,M,SUM,ERR,CONST 30 CONTINUE STOP 9000 FORMAT (5X,'n',4X,'m',8X,'integral',10X,'rel error',4X, + 'err const',/) 9010 FORMAT (1X,'ierab=',I3) 9020 FORMAT (1X,'iergq=',I3) 9030 FORMAT (1X,2I5,D23.14,E13.4,D13.4) END DOUBLE PRECISION FUNCTION F(T,ETA,THETA) C .. Scalar Arguments .. DOUBLE PRECISION ETA,T,THETA C .. C .. Intrinsic Functions .. INTRINSIC EXP,SQRT C .. F = SQRT(1.D0+.5D0*THETA*T)/ (EXP(-ETA)+EXP(-T)) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'driver6.f' then echo shar: will not over-write existing file "'driver6.f'" else cat << SHAR_EOF > 'driver6.f' c c PROGRAM DTEST6 c c This program is a double-precision version of test6 c c parameter(ak=1.5d0) c parameter(ak=2.5d0) C .. Parameters .. INTEGER NMAX,NCAPM PARAMETER (NMAX=10,NCAPM=200) DOUBLE PRECISION THETA PARAMETER (THETA=1.D-4) DOUBLE PRECISION ETA PARAMETER (ETA=-1.D0) DOUBLE PRECISION AK PARAMETER (AK=.5D0) C .. C .. Local Scalars .. DOUBLE PRECISION CONST,DEN,EPS,EXACT,PI,SUM REAL ERR INTEGER IERAB,IERGQ,IERR,IROUT,K,KOUNT,M,MU,N,NCAP C .. C .. Local Arrays .. DOUBLE PRECISION A(NCAPM),ALPHA(NMAX+1),B(NCAPM),BE(NMAX+1), + BETA(NMAX+1),E(NCAPM),P0(NCAPM),P1(NCAPM), + P2(NCAPM),W(NCAPM),WG(NMAX),X(NCAPM),XII(2*NMAX), + XIR(2*NMAX),ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. DOUBLE PRECISION D1MACH,F EXTERNAL D1MACH,F C .. C .. External Subroutines .. EXTERNAL DABMOD,DGQRAT,DRECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,DBLE,REAL integer i1mach, nout nout = i1mach(2) C .. IROUT = 1 c irout=0 EPS = D1MACH(3)*1.D2 CALL DRECUR(NCAPM,7,AK-1.D0,0.D0,A,B,IERR) PI = 4.D0*ATAN(1.D0) WRITE(NOUT,FMT=9000) DO 30 N = 1,NMAX M = 2*N - 1 c m=2*((n+1)/2)-1 c m=1 c m=0 IF (M.GE.3) THEN DO 10 MU = 1,M - 2,2 DEN = ETA**2 + DBLE((MU+1)**2)* (PI**2) XIR(MU) = -ETA/DEN XIR(MU+1) = XIR(MU) XII(MU) = DBLE(MU+1)*PI/DEN XII(MU+1) = -XII(MU) IS(MU) = 1 IS(MU+1) = 1 10 CONTINUE END IF IF (M.GT.0) THEN XIR(M) = -1.D0/ETA XII(M) = 0.D0 IS(M) = 1 END IF CALL DABMOD(N+1,NCAPM,M,EPS,IROUT,A,B,XIR,XII,IS,ALPHA,BETA, + NCAP,KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9010) IERAB CALL DGQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9020) IERGQ SUM = 0.D0 DO 20 K = 1,N SUM = SUM + WG(K)*F(ZG(K),ETA,THETA) 20 CONTINUE EXACT = .379708865998074D0 c exact=.52608888707965d0 c exact=1.266569126543118d0 ERR = ABS(REAL((SUM-EXACT)/EXACT)) WRITE(NOUT,FMT=9030) N,M,SUM,ERR,CONST 30 CONTINUE STOP 9000 FORMAT (5X,'n',4X,'m',8X,'integral',10X,'rel error',4X, + 'err const',/) 9010 FORMAT (1X,'ierab=',I3) 9020 FORMAT (1X,'iergq=',I3) 9030 FORMAT (1X,2I5,D23.14,E13.4,D13.4) END DOUBLE PRECISION FUNCTION F(T,ETA,THETA) C .. Scalar Arguments .. DOUBLE PRECISION ETA,T,THETA C .. C .. Local Scalars .. DOUBLE PRECISION S,S1,TERM,X INTEGER L C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,EXP,SQRT integer i1mach, nout C .. X = T - ETA IF (ABS(X).LE.1.D0) THEN L = 0 TERM = X S1 = TERM 10 S = S1 L = L + 1 IF (L.GT.200) GO TO 20 TERM = X*TERM/DBLE(L+1) S1 = S + TERM IF (S1.NE.S) GO TO 10 F = T*SQRT(1.D0+.5D0*THETA*T)/ (EXP(-T)*S) RETURN ELSE F = T*SQRT(1.D0+.5D0*THETA*T)/ (EXP(-ETA)-EXP(-T)) RETURN END IF 20 nout = i1mach(2) WRITE(NOUT,FMT=9000) RETURN 9000 FORMAT (1X,'exp series does not converge') END SHAR_EOF fi # end of overwriting check if test -f 'driver7.f' then echo shar: will not over-write existing file "'driver7.f'" else cat << SHAR_EOF > 'driver7.f' c c PROGRAM DTEST7 c c This is a double-precision version of test7 c c c Be sure that ncapm is greater or equalt to numax c c parameter(ak=1.5d0) c parameter(ak=2.5d0) C .. Parameters .. INTEGER NMAX,NCAPM,NUMAX PARAMETER (NMAX=10,NCAPM=100,NUMAX=50) DOUBLE PRECISION THETA PARAMETER (THETA=1.D-4) DOUBLE PRECISION AK PARAMETER (AK=.5D0) C .. C .. Local Scalars .. DOUBLE PRECISION C1,C2,CONST,DX,DY,EPS,ETA,EXACT,FMU,FMU2,FNU, + FNU2,HN,HP,P,PI,SQPI,SUM,SUM0,T,X1,XA REAL ERR INTEGER IERAB,IERGC,IERGQ,IERK,IERR,IROUT,J,K,KOUNT,M,MU,N,NCAP, + NU,NU0,NU1 C .. C .. Local Arrays .. DOUBLE PRECISION A(NCAPM),A1(NMAX+1),ALPHA(NMAX+1),B(NCAPM), + B1(NMAX+1),BE(NMAX+1),BETA(NMAX+1),E(NCAPM), + P0(NCAPM),P1(NCAPM),P2(NCAPM),RHOI(1),RHOR(1), + ROLDI(1),ROLDR(1),W(NCAPM),WG(NMAX),X(NCAPM), + XII(2*NMAX),XIR(2*NMAX),ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. DOUBLE PRECISION D1MACH,F EXTERNAL D1MACH,F C .. C .. External Subroutines .. EXTERNAL DABMOD,DGCHRS,DGQRAT,DKNUM,DRECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,DBLE,EXP,REAL,SQRT integer i1mach, nout nout = i1mach(2) C .. PI = 4.D0*ATAN(1.D0) SQPI = SQRT(PI) C1 = PI c c1=-pi c c1=pi C2 = SQPI c c2=.5*sqpi c c2=.75*sqpi IROUT = 1 c irout=0 EPS = D1MACH(3)*1.D3 CALL DRECUR(NCAPM,7,AK-1.D0,0.D0,A,B,IERR) WRITE(NOUT,FMT=9000) DO 60 N = 1,NMAX M = 2*N - 1 c m=2*((n+1)/2)-1 c m=1 IF (M.GE.3) THEN DO 10 MU = 1,M - 2,2 XIR(MU) = 0.D0 XIR(MU+1) = 0.D0 XII(MU) = 1.D0/ (DBLE(MU+1)*PI) XII(MU+1) = -XII(MU) IS(MU) = 1 IS(MU+1) = 1 10 CONTINUE END IF CALL DABMOD(N+1,NCAPM,M-1,EPS,IROUT,A,B,XIR,XII,IS,A1,B1,NCAP, + KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9010) IERAB ETA = -.001D0 XIR(M) = -1.D0/ETA XII(M) = 0.D0 IS(M) = 1 X1 = ETA XA = ABS(ETA) T = 1.D0 SUM0 = 0.D0 SUM = -1.D0/ (AK-1.D0) J = 0 20 J = J + 1 IF (J.GT.100) THEN WRITE(NOUT,FMT=9020) GO TO 60 END IF T = -XA*T/DBLE(J) SUM = SUM + T/ (DBLE(J+1)-AK) IF (SUM.NE.SUM0) THEN SUM0 = SUM GO TO 20 END IF HP = -EXP(XA)* (C1* (XA** (AK-1.D0))-C2*SUM) HN = 0.D0 IF (M.GE.3) THEN SUM = 0.D0 DO 40 NU = 1, (M-1)/2 FNU = DBLE(NU) FNU2 = FNU**2 P = 1.D0 DO 30 MU = 1, (M-1)/2 IF (MU.EQ.NU) GO TO 30 FMU = DBLE(MU) FMU2 = FMU**2 P = (FMU2/ (FMU2-FNU2))*P 30 CONTINUE NU0 = 10 DX = 0.D0 DY = 2.D0*FNU*PI CALL DKNUM(0,NU0,NUMAX,DX,DY,EPS,A,B,RHOR,RHOI,NU1, + IERK,ROLDR,ROLDI) IF (IERK.NE.0) WRITE(NOUT,FMT=9030) IERK SUM = SUM + 2.D0*FNU*PI*P* + (2.D0*FNU*PI*HP- (2.D0*FNU*PI*RHOR(1)+ + X1*RHOI(1)))/ (X1**2+4.D0*FNU2* (PI**2)) 40 CONTINUE HP = SUM END IF CALL DGCHRS(N+1,1,A1,B1,X1,HP,HN,ALPHA,BETA,IERGC) IF (IERGC.NE.0) WRITE(NOUT,FMT=9040) IERGC BETA(1) = -ETA*BETA(1) CALL DGQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9050) IERGQ SUM = 0.D0 DO 50 K = 1,N SUM = SUM + WG(K)*F(ZG(K),ETA,THETA) 50 CONTINUE EXACT = 2.2171501009112D0 c exact=1.7800123286291d0 c exact=3.7403844583705d0 ERR = ABS(REAL((SUM-EXACT)/EXACT)) WRITE(NOUT,FMT=9060) N,M,SUM,ERR,CONST 60 CONTINUE STOP 9000 FORMAT (5X,'n',4X,'m',8X,'integral',10X,'rel error',4X, + 'err const',/) 9010 FORMAT (1X,'ierab=',I3) 9020 FORMAT (1X,'power series for gamma does not converge') 9030 FORMAT (1X,'ierk for complex z =',I3) 9040 FORMAT (1X,'iergc=',I3) 9050 FORMAT (1X,'iergq=',I3) 9060 FORMAT (1X,2I5,D23.14,E13.4,D13.4) END DOUBLE PRECISION FUNCTION F(T,ETA,THETA) C .. Scalar Arguments .. DOUBLE PRECISION ETA,T,THETA C .. C .. Local Scalars .. DOUBLE PRECISION S,S1,TERM,X INTEGER L C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,EXP,SQRT integer i1mach, nout C .. X = T - ETA IF (ABS(X).LE.1.D0) THEN L = 0 TERM = X S1 = TERM 10 S = S1 L = L + 1 IF (L.GT.200) GO TO 20 TERM = X*TERM/DBLE(L+1) S1 = S + TERM IF (S1.NE.S) GO TO 10 F = T*SQRT(1.D0+.5D0*THETA*T)/ (EXP(-T)*S) RETURN ELSE F = T*SQRT(1.D0+.5D0*THETA*T)/ (EXP(-ETA)-EXP(-T)) RETURN END IF 20 nout = i1mach(2) WRITE(NOUT,FMT=9000) RETURN 9000 FORMAT (1X,'exp series does not converge') END SHAR_EOF fi # end of overwriting check if test -f 'res1' then echo shar: will not over-write existing file "'res1'" else cat << SHAR_EOF > 'res1' n m integral rel error err const 1 2 0.33489746814958D+01 0.2504E+00 0.8161D+00 2 4 0.44369266068371D+01 0.6904E-02 0.1730D-01 3 6 0.44674130628545D+01 0.8071E-04 0.1524D-03 4 8 0.44677713487831D+01 0.5143E-06 0.7027D-06 5 10 0.44677736372189D+01 0.2052E-08 0.1995D-08 6 12 0.44677736463627D+01 0.5600E-11 0.3838D-11 7 14 0.44677736463877D+01 0.1511E-13 0.5334D-14 8 16 0.44677736463878D+01 0.5964E-15 0.5607D-17 9 18 0.44677736463878D+01 0.1789E-14 0.4614D-20 10 20 0.44677736463878D+01 0.3976E-15 0.3054D-23 SHAR_EOF fi # end of overwriting check if test -f 'res2' then echo shar: will not over-write existing file "'res2'" else cat << SHAR_EOF > 'res2' n m integral rel error err const 1 2 0.76090037369184D+01 0.4115E+00 0.2810D+01 2 4 0.12820104334736D+02 0.8442E-02 0.3754D-01 3 6 0.12927897586232D+02 0.1051E-03 0.3170D-03 4 8 0.12929249074785D+02 0.6014E-06 0.1432D-05 5 10 0.12929256953344D+02 0.7993E-08 0.4018D-08 6 12 0.12929256859428D+02 0.7289E-09 0.7674D-11 7 14 0.12929256850633D+02 0.4870E-10 0.1061D-13 8 16 0.12929256850045D+02 0.3243E-11 0.1112D-16 9 18 0.12929256850005D+02 0.1870E-12 0.9128D-20 10 20 0.12929256850003D+02 0.2446E-13 0.6029D-23 SHAR_EOF fi # end of overwriting check if test -f 'res3' then echo shar: will not over-write existing file "'res3'" else cat << SHAR_EOF > 'res3' n m integral rel error err const 1 2 0.75464569543048D+01 0.1642E-01 0.1204D-01 om= 0.100 2 4 0.76676611672816D+01 0.6258E-03 0.4843D-04 3 6 0.76723634334878D+01 0.1292E-04 0.8441D-07 4 8 0.76724613152951D+01 0.1657E-06 0.8250D-10 5 10 0.76724625749979D+01 0.1482E-08 0.5157D-13 6 12 0.76724625862945D+01 0.9915E-11 0.2238D-16 7 14 0.76724625863702D+01 0.5800E-13 0.7136D-20 8 16 0.76724625863705D+01 0.6714E-14 0.1742D-23 9 18 0.76724625863706D+01 0.3820E-14 0.3360D-27 10 20 0.76724625863705D+01 0.7525E-14 0.5249D-31 1 2 0.24850912542716D+01 0.2665E-01 0.3101D-01 om= 0.500 2 4 0.25505961352628D+01 0.9953E-03 0.1217D-03 3 6 0.25530853848234D+01 0.2028E-04 0.2113D-06 4 8 0.25531364985082D+01 0.2581E-06 0.2062D-09 5 10 0.25531371515746D+01 0.2298E-08 0.1289D-12 6 12 0.25531371574028D+01 0.1530E-10 0.5591D-16 7 14 0.25531371574417D+01 0.7775E-13 0.1782D-19 8 16 0.25531371574419D+01 0.1218E-14 0.4350D-23 9 18 0.25531371574419D+01 0.1044E-14 0.8390D-27 10 20 0.25531371574419D+01 0.1913E-14 0.1311D-30 1 2 0.14335383524051D+01 0.3039E-01 0.4022D-01 om= 1.000 2 4 0.14768045526316D+01 0.1125E-02 0.1563D-03 3 6 0.14784335049727D+01 0.2279E-04 0.2709D-06 4 8 0.14784667785431D+01 0.2893E-06 0.2643D-09 5 10 0.14784672025085D+01 0.2572E-08 0.1651D-12 6 12 0.14784672062853D+01 0.1712E-10 0.7162D-16 7 14 0.14784672063104D+01 0.1023E-12 0.2283D-19 8 16 0.14784672063106D+01 0.1036E-13 0.5572D-23 9 18 0.14784672063106D+01 0.1141E-13 0.1075D-26 10 20 0.14784672063106D+01 0.1111E-13 0.1679D-30 1 2 0.79062183660300D+00 0.3313E-01 0.4783D-01 om= 1.999 2 4 0.81671889752784D+00 0.1217E-02 0.1844D-03 3 6 0.81769408432290D+00 0.2457E-04 0.3193D-06 4 8 0.81771392471853D+00 0.3113E-06 0.3114D-09 5 10 0.81771417700064D+00 0.2763E-08 0.1944D-12 6 12 0.81771417924499D+00 0.1835E-10 0.8434D-16 7 14 0.81771417925993D+00 0.9015E-13 0.2688D-19 8 16 0.81771417926000D+00 0.5974E-14 0.6561D-23 9 18 0.81771417926000D+00 0.4752E-14 0.1265D-26 10 20 0.81771417926000D+00 0.4752E-14 0.1977D-30 1 2 0.79027050508127D+00 0.3313E-01 0.4783D-01 om= 2.000 2 4 0.81635728109842D+00 0.1217E-02 0.1844D-03 3 6 0.81733207910043D+00 0.2458E-04 0.3193D-06 4 8 0.81735191153903D+00 0.3113E-06 0.3114D-09 5 10 0.81735216371965D+00 0.2763E-08 0.1945D-12 6 12 0.81735216596310D+00 0.1837E-10 0.8435D-16 7 14 0.81735216597803D+00 0.1009E-12 0.2689D-19 8 16 0.81735216597811D+00 0.4618E-14 0.6562D-23 9 18 0.81735216597811D+00 0.4347E-14 0.1265D-26 10 20 0.81735216597810D+00 0.6520E-14 0.1977D-30 1 2 0.54782395490603D+00 0.3427E-01 0.5121D-01 om= 3.000 2 4 0.56655256475675D+00 0.1255E-02 0.1968D-03 3 6 0.56725009686416D+00 0.2530E-04 0.3405D-06 4 8 0.56726426433430D+00 0.3201E-06 0.3320D-09 5 10 0.56726444432504D+00 0.2840E-08 0.2073D-12 6 12 0.56726444592537D+00 0.1887E-10 0.8993D-16 7 14 0.56726444593601D+00 0.1057E-12 0.2867D-19 8 16 0.56726444593607D+00 0.6067E-14 0.6996D-23 9 18 0.56726444593607D+00 0.6263E-14 0.1349D-26 10 20 0.56726444593607D+00 0.6850E-14 0.2108D-30 1 2 0.54765642732475D+00 0.3427E-01 0.5121D-01 om= 3.001 2 4 0.56637976755518D+00 0.1255E-02 0.1968D-03 3 6 0.56707710173669D+00 0.2530E-04 0.3405D-06 4 8 0.56709126516970D+00 0.3201E-06 0.3321D-09 5 10 0.56709144510904D+00 0.2840E-08 0.2073D-12 6 12 0.56709144670891D+00 0.1887E-10 0.8994D-16 7 14 0.56709144671955D+00 0.1010E-12 0.2867D-19 8 16 0.56709144671961D+00 0.2545E-14 0.6996D-23 9 18 0.56709144671961D+00 0.1175E-14 0.1349D-26 10 20 0.56709144671961D+00 0.2545E-14 0.2108D-30 1 2 0.16686705970311D+00 0.3618E-01 0.5719D-01 om=10.500 2 4 0.17290250057023D+00 0.1317E-02 0.2186D-03 3 6 0.17312598211460D+00 0.2648E-04 0.3778D-06 4 8 0.17313050811304D+00 0.3346E-06 0.3682D-09 5 10 0.17313056552844D+00 0.2965E-08 0.2299D-12 6 12 0.17313056603844D+00 0.1966E-10 0.9971D-16 7 14 0.17313056604183D+00 0.8385E-13 0.3178D-19 8 16 0.17313056604184D+00 0.1812E-13 0.7756D-23 9 18 0.17313056604184D+00 0.1924E-13 0.1496D-26 10 20 0.17313056604184D+00 0.1683E-13 0.2336D-30 1 2 0.35849757449253D-01 0.3687E-01 0.5945D-01 om=50.000 2 4 0.37172220145055D-01 0.1340E-02 0.2267D-03 3 6 0.37221081183120D-01 0.2690E-04 0.3917D-06 4 8 0.37222069673182D-01 0.3397E-06 0.3818D-09 5 10 0.37222082206027D-01 0.3010E-08 0.2384D-12 6 12 0.37222082317311D-01 0.1997E-10 0.1034D-15 7 14 0.37222082318050D-01 0.1016E-12 0.3295D-19 8 16 0.37222082318054D-01 0.1864E-14 0.8041D-23 9 18 0.37222082318054D-01 0.9321E-15 0.1551D-26 10 20 0.37222082318054D-01 0.7457E-15 0.2422D-30 SHAR_EOF fi # end of overwriting check if test -f 'res5' then echo shar: will not over-write existing file "'res5'" else cat << SHAR_EOF > 'res5' n m integral rel error err const 1 2 0.27540902998306D+00 0.5199E-01 0.2782D+00 2 4 0.29012218422625D+00 0.1343E-02 0.5457D-01 3 6 0.29050451350811D+00 0.2721E-04 0.9089D-02 4 8 0.29051227564405D+00 0.4866E-06 0.1394D-02 5 10 0.29051241467479D+00 0.8071E-08 0.2033D-03 6 12 0.29051241698254D+00 0.1272E-09 0.2868D-04 7 14 0.29051241701893D+00 0.1925E-11 0.3950D-05 8 16 0.29051241701948D+00 0.2217E-13 0.5343D-06 9 18 0.29051241701949D+00 0.8216E-14 0.7128D-07 10 20 0.29051241701949D+00 0.7070E-14 0.9403D-08 SHAR_EOF fi # end of overwriting check if test -f 'res6' then echo shar: will not over-write existing file "'res6'" else cat << SHAR_EOF > 'res6' n m integral rel error err const 1 1 0.28430903156845D+00 0.2512E+00 0.1600D+00 2 3 0.37673641714434D+00 0.7828E-02 0.1738D-01 3 5 0.37964411940764D+00 0.1705E-03 0.2003D-02 4 7 0.37970766344352D+00 0.3167E-05 0.2344D-03 5 9 0.37970884559901D+00 0.5372E-07 0.2762D-04 6 11 0.37970886567178D+00 0.8593E-09 0.3265D-05 7 13 0.37970886599306D+00 0.1319E-10 0.3869D-06 8 15 0.37970886599800D+00 0.1977E-12 0.4592D-07 9 17 0.37970886599807D+00 0.3070E-14 0.5454D-08 10 19 0.37970886599807D+00 0.3070E-14 0.6482D-09 SHAR_EOF fi # end of overwriting check if test -f 'res7' then echo shar: will not over-write existing file "'res7'" else cat << SHAR_EOF > 'res7' n m integral rel error err const 1 1 0.16904066964116D+01 0.2376E+00 0.4276D-03 2 3 0.22011951368759D+01 0.7196E-02 0.3995D-04 3 5 0.22168095642127D+01 0.1536E-03 0.4279D-05 4 7 0.22171438690047D+01 0.2811E-05 0.4785D-06 5 9 0.22171499969530D+01 0.4689E-07 0.5462D-07 6 11 0.22171500993576D+01 0.7007E-09 0.6305D-08 7 13 0.22171501008772D+01 0.1533E-10 0.7332D-09 8 15 0.22171501009074D+01 0.1703E-11 0.8568D-10 9 17 0.22171501009114D+01 0.7892E-13 0.1005D-10 10 19 0.22171501009114D+01 0.1022E-12 0.1181D-11 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'driver1.f' then echo shar: will not over-write existing file "'driver1.f'" else cat << SHAR_EOF > 'driver1.f' c c PROGRAM TEST1 c c This program computes the n-point rational Gauss quadrature c approximations, n=1,2,...,nmax, to the integral in Example 4.1 c using rational functions that share with the integrand the m c poles closest to the interval [-1,1]. c C .. Parameters .. INTEGER NMAX,NCAPM REAL OM PARAMETER (NMAX=10,NCAPM=50,OM=1.1) C .. C .. Local Scalars .. REAL CONST,EPS,ERR,EXACT,SGN,SUM INTEGER IERAB,IERGQ,IERR,IROUT,K,KOUNT,M,MU,N,NCAP C .. C .. Local Arrays .. REAL A(NCAPM),ALPHA(NMAX+1),B(NCAPM),BE(NMAX+1),BETA(NMAX+1), + E(NCAPM),P0(NCAPM),P1(NCAPM),P2(NCAPM),W(NCAPM),WG(NMAX), + X(NCAPM),XII(2*NMAX),XIR(2*NMAX),ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. REAL F,R1MACH EXTERNAL F,R1MACH C .. C .. External Subroutines .. EXTERNAL ABMOD,GQRAT,RECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,REAL integer i1mach, nout nout = i1mach(2) C .. IROUT = 1 c irout=0 EPS = R1MACH(3)*1.e2 CALL RECUR(NCAPM,1,0.,0.,A,B,IERR) WRITE(NOUT,FMT=9000) DO 30 N = 1,NMAX M = 2*N c m=2*((n+1)/2) c m=2 c m=0 SGN = 1. DO 10 MU = 1,M SGN = -SGN XIR(MU) = SGN/ (OM*REAL((MU+1)/2)) XII(MU) = 0. IS(MU) = 1 10 CONTINUE CALL ABMOD(N+1,NCAPM,M,EPS,IROUT,A,B,XIR,XII,IS,ALPHA,BETA, + NCAP,KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9010) IERAB CALL GQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9020) IERGQ SUM = 0. DO 20 K = 1,N SUM = SUM + WG(K)*F(ZG(K),OM) 20 CONTINUE EXACT = 4.4677736 ERR = ABS((SUM-EXACT)/EXACT) WRITE(NOUT,FMT=9030) N,M,SUM,ERR,CONST 30 CONTINUE STOP 9000 FORMAT (5X,'n',4X,'m',6X,'integral',7X,'rel error',5X,'err const', + /) 9010 FORMAT (1X,'ierab=',I3) 9020 FORMAT (1X,'iergq=',I3) 9030 FORMAT (1X,2I5,E17.7,2E14.4) END REAL FUNCTION F(T,OM) C .. Scalar Arguments .. REAL OM,T C .. C .. Local Scalars .. REAL PI C .. C .. Intrinsic Functions .. INTRINSIC ATAN,SIN C .. PI = 4.*ATAN(1.) F = 1. IF (T.NE.0.) F = (PI*T/OM)/SIN(PI*T/OM) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'driver2.f' then echo shar: will not over-write existing file "'driver2.f'" else cat << SHAR_EOF > 'driver2.f' c c PROGRAM TEST2 c c This program does the same as test1 but treats the two c (symmetric) poles closest to the interval [-1,1] separately. c C .. Parameters .. INTEGER NMAX,NCAPM REAL OM PARAMETER (NMAX=10,NCAPM=50,OM=1.001) C .. C .. Local Scalars .. REAL CONST,EPS,ERR,EXACT,HN,HP,P,SGN,SUM,SUMN,SUMP,X1 INTEGER IERAB,IERGC,IERGQ,IERR,IROUT,K,KOUNT,M,MU,N,NCAP,NU C .. C .. Local Arrays .. REAL A(NCAPM),A1(NMAX+1),ALPHA(NMAX+1),B(NCAPM),B1(NMAX+1), + BE(NMAX+1),BETA(NMAX+1),E(NCAPM),P0(NCAPM),P1(NCAPM), + P2(NCAPM),W(NCAPM),WG(NMAX),X(NCAPM),XII(2*NMAX),XIR(2*NMAX), + ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. REAL F,R1MACH EXTERNAL F,R1MACH C .. C .. External Subroutines .. EXTERNAL ABMOD,GCHRIS,GQRAT,RECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,LOG,REAL integer i1mach, nout nout = i1mach(2) C .. IROUT = 1 c irout=0 EPS = R1MACH(3)*1.e2 CALL RECUR(NCAPM,1,0.,0.,A,B,IERR) WRITE(NOUT,FMT=9000) DO 50 N = 1,NMAX M = 2*N c m=n c m=4 c m=2 SGN = 1. DO 10 NU = 1,M SGN = -SGN IF (NU.LE.M-2) XIR(NU) = SGN/REAL((NU+3)/2) XII(NU) = 0. IS(NU) = 1 10 CONTINUE CALL ABMOD(N+1,NCAPM,M-2,EPS,IROUT,A,B,XIR,XII,IS,A1,B1,NCAP, + KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9010) IERAB XIR(M-1) = 1./OM XIR(M) = -1./OM X1 = OM HP = LOG(ABS((X1+1.)/ (X1-1.))) HN = -HP IF (M.GT.2) THEN SUMP = 0. SUMN = 0. DO 30 NU = 1,M - 2 P = 1. DO 20 MU = 1,M - 2 IF (MU.EQ.NU) GO TO 20 P = (1.-XIR(MU)/XIR(NU))*P 20 CONTINUE SUMP = SUMP + LOG(ABS((X1+1.)* (XIR(NU)+1.)/ ((X1- + 1.)* (XIR(NU)-1.))))/ ((1.+XIR(NU)*X1)*P) SUMN = SUMN + LOG(ABS((X1-1.)* (XIR(NU)+1.)/ ((X1+ + 1.)* (XIR(NU)-1.))))/ ((1.-XIR(NU)*X1)*P) 30 CONTINUE HP = SUMP HN = SUMN END IF CALL GCHRIS(N+1,2,A1,B1,X1,HP,HN,ALPHA,BETA,IERGC) BETA(1) = - (OM**2)*BETA(1) CALL GQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9020) IERGQ SUM = 0. DO 40 K = 1,N SUM = SUM + WG(K)*F(ZG(K),OM) 40 CONTINUE EXACT = 12.929257 ERR = ABS((SUM-EXACT)/EXACT) WRITE(NOUT,FMT=9030) N,M,SUM,ERR,CONST 50 CONTINUE STOP 9000 FORMAT (5X,'n',4X,'m',6X,'integral',7X,'rel error',5X,'err const', + /) 9010 FORMAT (1X,'ierab=',I3) 9020 FORMAT (1X,'iergq=',I3) 9030 FORMAT (1X,2I5,E17.7,2E14.4) END REAL FUNCTION F(T,OM) C .. Scalar Arguments .. REAL OM,T C .. C .. Local Scalars .. REAL PI C .. C .. Intrinsic Functions .. INTRINSIC ATAN,SIN C .. PI = 4.*ATAN(1.) F = 1. IF (T.NE.0.) F = (PI*T/OM)/SIN(PI*T/OM) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'driver3.f' then echo shar: will not over-write existing file "'driver3.f'" else cat << SHAR_EOF > 'driver3.f' c c PROGRAM TEST3 c c This program computes the n-point rational Gauss quadrature c approximations, n=1,2,...,nmax, to the integral in Example 4.2 c for om=.1, .5, 1., 1.999, 2., 3., 3.001, 10.5, and 50., using c rational functions that share with the integrand the m poles c closest to, and to the left of, the interval [0,1]. c C .. Parameters .. INTEGER NMAX,NCAPM PARAMETER (NMAX=10,NCAPM=50) C .. C .. Local Scalars .. REAL CONST,EPS,ERR,OM,SUM INTEGER IERAB,IERGQ,IERR,IOM,IROUT,K,KOUNT,M,MU,N,NCAP C .. C .. Local Arrays .. REAL A(NCAPM),ALPHA(NMAX+1),B(NCAPM),BE(NMAX+1),BETA(NMAX+1), + E(NCAPM),EXACT(9),OOM(9),P0(NCAPM),P1(NCAPM),P2(NCAPM), + W(NCAPM),WG(NMAX),X(NCAPM),XII(2*NMAX),XIR(2*NMAX),ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. REAL F,R1MACH EXTERNAL F,R1MACH C .. C .. External Subroutines .. EXTERNAL ABMOD,GQRAT,RECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,REAL,SQRT integer i1mach, nout C .. C .. Data statements .. DATA OOM/.1,.5,1.,1.999,2.,3.,3.001,10.5,50./ DATA EXACT/7.6724626,2.5531372,1.4784672,.81771418,.81735217, + .56726445,.56709145,.17313057,.037222082/ C .. nout = i1mach(2) WRITE(NOUT,FMT=9000) IROUT = 1 c irout=0 EPS = R1MACH(3)*1.e2 CALL RECUR(NCAPM,6,0.,-.5,A,B,IERR) DO 10 K = 1,NCAPM A(K) = .5* (1.+A(K)) B(K) = .25*B(K) 10 CONTINUE B(1) = SQRT(8.)*B(1) WRITE(NOUT,FMT=9010) DO 50 IOM = 1,9 OM = OOM(IOM) DO 40 N = 1,NMAX M = 2*N c m=n c m=2 c m=0 IF (M.GT.0) THEN DO 20 MU = 1,M IF (MU.EQ.1) THEN XIR(MU) = 1./OM ELSE XIR(MU) = 1./REAL(MU-1) END IF XII(MU) = 0. IS(MU) = 1 20 CONTINUE END IF CALL ABMOD(N+1,NCAPM,M,EPS,IROUT,A,B,XIR,XII,IS,ALPHA, + BETA,NCAP,KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9020) IERAB CALL GQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9030) IERGQ SUM = 0. DO 30 K = 1,N SUM = SUM + WG(K)*F(ZG(K),OM) 30 CONTINUE ERR = ABS((SUM-EXACT(IOM))/EXACT(IOM)) IF (N.EQ.1) THEN WRITE(NOUT,FMT=9040) N,M,SUM,ERR,CONST,OM ELSE WRITE(NOUT,FMT=9050) N,M,SUM,ERR,CONST END IF 40 CONTINUE WRITE(NOUT,FMT=9000) 50 CONTINUE STOP 9000 FORMAT (/) 9010 FORMAT (5X,'n',4X,'m',6X,'integral',7X,'rel error',5X,'err const', + /) 9020 FORMAT (1X,'ierab=',I3) 9030 FORMAT (1X,'iergq=',I3) 9040 FORMAT (1X,2I5,E17.7,2E14.4,' om=',F6.3) 9050 FORMAT (1X,2I5,E17.7,2E14.4) END REAL FUNCTION F(T,OM) C .. Scalar Arguments .. REAL OM,T C .. C .. Local Scalars .. INTEGER IERR C .. C .. External Functions .. REAL GAMMA EXTERNAL GAMMA C .. F = GAMMA(1.+T,IERR)/ (T+OM) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'driver4.f' then echo shar: will not over-write existing file "'driver4.f'" else cat << SHAR_EOF > 'driver4.f' c c PROGRAM TEST4 c c This program does the same as test3 but treats the pole closest c to the interval [0,1] separately. c C .. Parameters .. INTEGER NMAX,NCAPM PARAMETER (NMAX=10,NCAPM=50) C .. C .. Local Scalars .. DOUBLE PRECISION DHP,DOM,DP,DSUM,DX1 REAL CONST,EPS,ERR,EXACT,HN,HP,OM,SUM,X1 INTEGER IERAB,IERGC,IERGQ,IERR,IROUT,K,KOUNT,M,MU,N,NCAP,NU C .. C .. Local Arrays .. DOUBLE PRECISION DXIR(2*NMAX) REAL A(NCAPM),A1(NMAX+1),ALPHA(NMAX+1),B(NCAPM),B1(NMAX+1), + BE(NMAX+1),BETA(NMAX+1),E(NCAPM),P0(NCAPM),P1(NCAPM), + P2(NCAPM),W(NCAPM),WG(NMAX),X(NCAPM),XII(2*NMAX),XIR(2*NMAX), + ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. REAL F,R1MACH EXTERNAL F,R1MACH C .. C .. External Subroutines .. EXTERNAL ABMOD,GCHRIS,GQRAT,RECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,DBLE,REAL,SQRT integer i1mach, nout nout = i1mach(2) C .. IROUT = 1 c irout=0 EPS = R1MACH(3)*1.e2 CALL RECUR(NCAPM,6,0.,-.5,A,B,IERR) DO 10 K = 1,NCAPM A(K) = .5* (1.+A(K)) B(K) = .25*B(K) 10 CONTINUE B(1) = SQRT(8.)*B(1) WRITE(NOUT,FMT=9000) DO 60 N = 1,NMAX M = 2*N c m=n c m=1 IF (M.GT.0) THEN DO 20 NU = 1,M DXIR(NU) = 1.D0/DBLE(NU) XIR(NU) = REAL(DXIR(NU)) XII(NU) = 0. IS(NU) = 1 20 CONTINUE END IF CALL ABMOD(N+1,NCAPM,M-1,EPS,IROUT,A,B,XIR,XII,IS,A1,B1,NCAP, + KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9010) IERAB,N DOM = .001D0 OM = REAL(DOM) DXIR(M) = 1.D0/DOM XIR(M) = REAL(DXIR(M)) DX1 = -DOM X1 = REAL(DX1) DHP = -2.D0*ATAN(1.D0/SQRT(ABS(DX1)))/SQRT(ABS(DX1)) HP = REAL(DHP) HN = 0. IF (M.GT.1) THEN DSUM = 0.D0 DO 40 NU = 1,M - 1 DP = 1.D0 DO 30 MU = 1,M - 1 IF (MU.EQ.NU) GO TO 30 DP = (1.D0-DXIR(MU)/DXIR(NU))*DP 30 CONTINUE DSUM = DSUM + (DHP+2.D0*SQRT(DXIR(NU))* + ATAN(SQRT(DXIR(NU))))/ ((1.D0+DXIR(NU)*DX1)*DP) 40 CONTINUE DHP = DSUM HP = REAL(DHP) END IF CALL GCHRIS(N+1,1,A1,B1,X1,HP,HN,ALPHA,BETA,IERGC) IF (IERGC.NE.0) WRITE(NOUT,FMT=9020) IERGC,N BETA(1) = OM*BETA(1) CALL GQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9030) IERGQ,N SUM = 0. DO 50 K = 1,N SUM = SUM + WG(K)*F(ZG(K),OM) 50 CONTINUE EXACT = 96.703688 ERR = ABS((SUM-EXACT)/EXACT) WRITE(NOUT,FMT=9040) N,M,SUM,ERR,CONST 60 CONTINUE STOP 9000 FORMAT (5X,'n',4X,'m',6X,'integral',7X,'rel error',5X,'err const', + /) 9010 FORMAT (1X,'ierab=',I3,' for n=',I2) 9020 FORMAT (1X,'iergc=',I3,' for n=',I2) 9030 FORMAT (1X,'iergq=',I3,' for n=',I2) 9040 FORMAT (1X,2I5,E17.7,2E14.4) END REAL FUNCTION F(T,OM) C .. Scalar Arguments .. REAL OM,T C .. C .. Local Scalars .. INTEGER IERR C .. C .. External Functions .. REAL GAMMA EXTERNAL GAMMA C .. F = GAMMA(1.+T,IERR)/ (T+OM) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'driver5.f' then echo shar: will not over-write existing file "'driver5.f'" else cat << SHAR_EOF > 'driver5.f' c c PROGRAM TEST5 c c This program computes the n-point rational Gauss quadrature c approximations, n=1,2,...,nmax, to the Fermi-Dirac integral c using rational functions that share with the integrand the m c conjugate complex poles closest to the interval [0,oo]. c c parameter(eta=1.) c parameter(ak=1.5) c parameter(ak=2.5) C .. Parameters .. INTEGER NMAX,NCAPM REAL THETA PARAMETER (NMAX=10,NCAPM=100,THETA=1.e-4) REAL ETA PARAMETER (ETA=-1.) REAL AK PARAMETER (AK=.5) C .. C .. Local Scalars .. REAL CONST,DEN,EPS,ERR,EXACT,PI,SUM INTEGER IERAB,IERGQ,IERR,IROUT,K,KOUNT,M,MU,N,NCAP C .. C .. Local Arrays .. REAL A(NCAPM),ALPHA(NMAX+1),B(NCAPM),BE(NMAX+1),BETA(NMAX+1), + E(NCAPM),P0(NCAPM),P1(NCAPM),P2(NCAPM),W(NCAPM),WG(NMAX), + X(NCAPM),XII(2*NMAX),XIR(2*NMAX),ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. REAL F,R1MACH EXTERNAL F,R1MACH C .. C .. External Subroutines .. EXTERNAL ABMOD,GQRAT,RECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,REAL integer i1mach, nout nout = i1mach(2) C .. IROUT = 1 c irout=0 EPS = R1MACH(3)*1.e2 CALL RECUR(NCAPM,7,AK,0.,A,B,IERR) PI = 4.*ATAN(1.) WRITE(NOUT,FMT=9000) DO 30 N = 1,NMAX M = 2*N c m=2*((n+1)/2) c m=2 c m=0 IF (M.GT.0) THEN DO 10 MU = 1,M - 1,2 DEN = ETA**2 + REAL(MU**2)* (PI**2) XIR(MU) = -ETA/DEN XIR(MU+1) = XIR(MU) XII(MU) = REAL(MU)*PI/DEN XII(MU+1) = -XII(MU) IS(MU) = 1 IS(MU+1) = 1 10 CONTINUE END IF CALL ABMOD(N+1,NCAPM,M,EPS,IROUT,A,B,XIR,XII,IS,ALPHA,BETA, + NCAP,KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9010) IERAB CALL GQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9020) IERGQ SUM = 0. DO 20 K = 1,N SUM = SUM + WG(K)*F(ZG(K),ETA,THETA) 20 CONTINUE EXACT = .29051242 c exact=.46087845 c exact=1.18607350 c exact=1.39644182 c exact=2.66187328 c exact=7.62725610 ERR = ABS((SUM-EXACT)/EXACT) WRITE(NOUT,FMT=9030) N,M,SUM,ERR,CONST 30 CONTINUE STOP 9000 FORMAT (5X,'n',4X,'m',6X,'integral',7X,'rel error',5X,'err const', + /) 9010 FORMAT (1X,'ierab=',I3) 9020 FORMAT (1X,'iergq=',I3) 9030 FORMAT (1X,2I5,E17.7,2E14.4) END REAL FUNCTION F(T,ETA,THETA) C .. Scalar Arguments .. REAL ETA,T,THETA C .. C .. Intrinsic Functions .. INTRINSIC EXP,SQRT C .. F = SQRT(1.+.5*THETA*T)/ (EXP(-ETA)+EXP(-T)) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'driver6.f' then echo shar: will not over-write existing file "'driver6.f'" else cat << SHAR_EOF > 'driver6.f' c c PROGRAM TEST6 c c This program computes the n-point rational Gauss quadrature c approximations, n=1,2,...,nmax, to the Bose-Einstein integral c using rational functions that share with the integrand m-1 c conjugate complex poles closest to the interval [0,oo] and the c real pole on the negative real axis. c c parameter(ak=1.5) c parameter(ak=2.5) C .. Parameters .. INTEGER NMAX,NCAPM REAL THETA PARAMETER (NMAX=10,NCAPM=100,THETA=1.e-4) REAL ETA PARAMETER (ETA=-1.) REAL AK PARAMETER (AK=.5) C .. C .. Local Scalars .. REAL CONST,DEN,EPS,ERR,EXACT,PI,SUM INTEGER IERAB,IERGQ,IERR,IROUT,K,KOUNT,M,MU,N,NCAP C .. C .. Local Arrays .. REAL A(NCAPM),ALPHA(NMAX+1),B(NCAPM),BE(NMAX+1),BETA(NMAX+1), + E(NCAPM),P0(NCAPM),P1(NCAPM),P2(NCAPM),W(NCAPM),WG(NMAX), + X(NCAPM),XII(2*NMAX),XIR(2*NMAX),ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. REAL F,R1MACH EXTERNAL F,R1MACH C .. C .. External Subroutines .. EXTERNAL ABMOD,GQRAT,RECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,REAL integer i1mach, nout nout = i1mach(2) C .. IROUT = 1 c irout=0 EPS = R1MACH(3)*1.e2 CALL RECUR(NCAPM,7,AK-1.,0.,A,B,IERR) PI = 4.*ATAN(1.) WRITE(NOUT,FMT=9000) DO 30 N = 1,NMAX M = 2*N - 1 c m=2*((n+1)/2)-1 c m=1 c m=0 IF (M.GE.3) THEN DO 10 MU = 1,M - 2,2 DEN = ETA**2 + REAL((MU+1)**2)* (PI**2) XIR(MU) = -ETA/DEN XIR(MU+1) = XIR(MU) XII(MU) = REAL(MU+1)*PI/DEN XII(MU+1) = -XII(MU) IS(MU) = 1 IS(MU+1) = 1 10 CONTINUE END IF IF (M.GT.0) THEN XIR(M) = -1./ETA XII(M) = 0. IS(M) = 1 END IF CALL ABMOD(N+1,NCAPM,M,EPS,IROUT,A,B,XIR,XII,IS,ALPHA,BETA, + NCAP,KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9010) IERAB CALL GQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9020) IERGQ SUM = 0. DO 20 K = 1,N SUM = SUM + WG(K)*F(ZG(K),ETA,THETA) 20 CONTINUE EXACT = .37970887 c exact=.52608889 c exact=1.26656913 ERR = ABS((SUM-EXACT)/EXACT) WRITE(NOUT,FMT=9030) N,M,SUM,ERR,CONST 30 CONTINUE STOP 9000 FORMAT (5X,'n',4X,'m',6X,'integral',7X,'rel error',5X,'err const', + /) 9010 FORMAT (1X,'ierab=',I3) 9020 FORMAT (1X,'iergq=',I3) 9030 FORMAT (1X,2I5,E17.7,2E14.4) END REAL FUNCTION F(T,ETA,THETA) C .. Scalar Arguments .. REAL ETA,T,THETA C .. C .. Local Scalars .. REAL S,S1,TERM,X INTEGER L C .. C .. Intrinsic Functions .. INTRINSIC ABS,EXP,REAL,SQRT integer i1mach, nout C .. X = T - ETA IF (ABS(X).LE.1.) THEN L = 0 TERM = X S1 = TERM 10 S = S1 L = L + 1 IF (L.GT.200) GO TO 20 TERM = X*TERM/REAL(L+1) S1 = S + TERM IF (S1.NE.S) GO TO 10 F = T*SQRT(1.+.5*THETA*T)/ (EXP(-T)*S) RETURN ELSE F = T*SQRT(1.+.5*THETA*T)/ (EXP(-ETA)-EXP(-T)) RETURN END IF 20 nout = i1mach(2) WRITE(NOUT,FMT=9000) RETURN 9000 FORMAT (1X,'exp series does not converge') END SHAR_EOF fi # end of overwriting check if test -f 'driver7.f' then echo shar: will not over-write existing file "'driver7.f'" else cat << SHAR_EOF > 'driver7.f' c c PROGRAM TEST7 c c This program does the same as test6 but treats the real pole c close to the origin separately. c c c Be sure that ncapm is greater or equal to kapmax c c parameter(ak=1.5) c parameter(ak=2.5) C .. Parameters .. INTEGER NMAX,NCAPM,KAPMAX REAL THETA PARAMETER (NMAX=10,NCAPM=50,KAPMAX=20,THETA=1.e-4) REAL AK PARAMETER (AK=.5) C .. C .. Local Scalars .. COMPLEX Z REAL C1,C2,CONST,EPS,ERR,ETA,EXACT,FMU,FMU2,FNU,FNU2,HN,HP,P,PI, + SQPI,SUM,SUM0,T,X1,XA INTEGER IERAB,IERGC,IERGQ,IERK,IERR,IROUT,J,K,KOUNT,M,MU,N,NCAP, + NU,NU0,NU1 C .. C .. Local Arrays .. COMPLEX RHO(1),ROLD(1) REAL A(NCAPM),A1(NMAX+1),ALPHA(NMAX+1),B(NCAPM),B1(NMAX+1), + BE(NMAX+1),BETA(NMAX+1),E(NCAPM),P0(NCAPM),P1(NCAPM), + P2(NCAPM),W(NCAPM),WG(NMAX),X(NCAPM),XII(2*NMAX),XIR(2*NMAX), + ZG(NMAX) INTEGER IS(2*NMAX) C .. C .. External Functions .. REAL F,R1MACH EXTERNAL F,R1MACH C .. C .. External Subroutines .. EXTERNAL ABMOD,GCHRIS,GQRAT,KNUM,RECUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,CMPLX,EXP,REAL,SQRT integer i1mach, nout nout = i1mach(2) C .. PI = 4.*ATAN(1.) SQPI = SQRT(PI) C1 = PI c c1=-pi c c1=pi C2 = SQPI c c2=.5*sqpi c c2=.75*sqpi IROUT = 1 c irout=0 EPS = R1MACH(3)*1.e2 CALL RECUR(NCAPM,7,AK-1.,0.,A,B,IERR) WRITE(NOUT,FMT=9000) DO 60 N = 1,NMAX M = 2*N - 1 c m=2*((n+1)/2)-1 c m=1 IF (M.GE.3) THEN DO 10 MU = 1,M - 2,2 XIR(MU) = 0. XIR(MU+1) = 0. XII(MU) = 1./ (REAL(MU+1)*PI) XII(MU+1) = -XII(MU) IS(MU) = 1 IS(MU+1) = 1 10 CONTINUE END IF CALL ABMOD(N+1,NCAPM,M-1,EPS,IROUT,A,B,XIR,XII,IS,A1,B1,NCAP, + KOUNT,IERAB,BE,X,W,E,P0,P1,P2) IF (IERAB.NE.0) WRITE(NOUT,FMT=9010) IERAB ETA = -.001 XIR(M) = -1./ETA XII(M) = 0. IS(M) = 1 X1 = ETA XA = ABS(ETA) T = 1. SUM0 = 0. SUM = -1./ (AK-1.) J = 0 20 J = J + 1 IF (J.GT.100) THEN WRITE(NOUT,FMT=9020) GO TO 60 END IF T = -XA*T/REAL(J) SUM = SUM + T/ (REAL(J+1)-AK) IF (SUM.NE.SUM0) THEN SUM0 = SUM GO TO 20 END IF HP = -EXP(XA)* (C1* (XA** (AK-1.))-C2*SUM) HN = 0. IF (M.GE.3) THEN SUM = 0. DO 40 NU = 1, (M-1)/2 FNU = REAL(NU) FNU2 = FNU**2 P = 1. DO 30 MU = 1, (M-1)/2 IF (MU.EQ.NU) GO TO 30 FMU = REAL(MU) FMU2 = FMU**2 P = (FMU2/ (FMU2-FNU2))*P 30 CONTINUE NU0 = 10 Z = CMPLX(0.,2.*FNU*PI) CALL KNUM(0,NU0,KAPMAX,Z,EPS,A,B,RHO,NU1,IERK,ROLD) IF (IERK.NE.0) WRITE(NOUT,FMT=9030) IERK SUM = SUM + 2.*FNU*PI*P* (2.*FNU*PI*HP- + REAL(CMPLX(2.*FNU*PI,-X1)*RHO(1)))/ + (X1**2+4.*FNU2* (PI**2)) 40 CONTINUE HP = SUM END IF CALL GCHRIS(N+1,1,A1,B1,X1,HP,HN,ALPHA,BETA,IERGC) IF (IERGC.NE.0) WRITE(NOUT,FMT=9040) IERGC BETA(1) = -ETA*BETA(1) CALL GQRAT(N,M,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERGQ,E) IF (IERGQ.NE.0) WRITE(NOUT,FMT=9050) IERGQ SUM = 0. DO 50 K = 1,N SUM = SUM + WG(K)*F(ZG(K),ETA,THETA) 50 CONTINUE EXACT = 2.2171501 c exact=1.7800123 c exact=3.7403845 ERR = ABS((SUM-EXACT)/EXACT) WRITE(NOUT,FMT=9060) N,M,SUM,ERR,CONST 60 CONTINUE STOP 9000 FORMAT (5X,'n',4X,'m',6X,'integral',7X,'rel error',5X,'err const', + /) 9010 FORMAT (1X,'ierab=',I3) 9020 FORMAT (1X,'power series for gamma does not converge') 9030 FORMAT (1X,'ierk for complex z =',I3) 9040 FORMAT (1X,'iergc=',I3) 9050 FORMAT (1X,'iergq=',I3) 9060 FORMAT (1X,2I5,E17.7,2E14.4) END REAL FUNCTION F(T,ETA,THETA) C .. Scalar Arguments .. REAL ETA,T,THETA C .. C .. Local Scalars .. REAL S,S1,TERM,X INTEGER L C .. C .. Intrinsic Functions .. INTRINSIC ABS,EXP,REAL,SQRT integer i1mach, nout C .. X = T - ETA IF (ABS(X).LE.1.) THEN L = 0 TERM = X S1 = TERM 10 S = S1 L = L + 1 IF (L.GT.200) GO TO 20 TERM = X*TERM/REAL(L+1) S1 = S + TERM IF (S1.NE.S) GO TO 10 F = T*SQRT(1.+.5*THETA*T)/ (EXP(-T)*S) RETURN ELSE F = T*SQRT(1.+.5*THETA*T)/ (EXP(-ETA)-EXP(-T)) RETURN END IF 20 nout = i1mach(2) WRITE(NOUT,FMT=9000) RETURN 9000 FORMAT (1X,'exp series does not converge') END SHAR_EOF fi # end of overwriting check if test -f 'res1' then echo shar: will not over-write existing file "'res1'" else cat << SHAR_EOF > 'res1' n m integral rel error err const 1 2 0.3348973E+01 0.2504E+00 0.8161E+00 2 4 0.4436926E+01 0.6904E-02 0.1730E-01 3 6 0.4467412E+01 0.8090E-04 0.1524E-03 4 8 0.4467774E+01 0.2135E-06 0.7027E-06 5 10 0.4467771E+01 0.5336E-06 0.1995E-08 6 12 0.4467773E+01 0.1067E-06 0.3838E-11 7 14 0.4467776E+01 0.5336E-06 0.5334E-14 8 16 0.4467773E+01 0.0000E+00 0.5607E-17 9 18 0.4467775E+01 0.4269E-06 0.4614E-20 10 20 0.4467773E+01 0.1067E-06 0.3054E-23 SHAR_EOF fi # end of overwriting check if test -f 'res2' then echo shar: will not over-write existing file "'res2'" else cat << SHAR_EOF > 'res2' n m integral rel error err const 1 2 0.7608957E+01 0.4115E+00 0.2810E+01 2 4 0.1282002E+02 0.8449E-02 0.3754E-01 3 6 0.1292781E+02 0.1122E-03 0.3170E-03 4 8 0.1292914E+02 0.9073E-05 0.1432E-05 5 10 0.1292918E+02 0.5975E-05 0.4018E-08 6 12 0.1292916E+02 0.7450E-05 0.7674E-11 7 14 0.1292916E+02 0.7229E-05 0.1061E-13 8 16 0.1292918E+02 0.5680E-05 0.1112E-16 9 18 0.1292918E+02 0.5753E-05 0.9128E-20 10 20 0.1292914E+02 0.9073E-05 0.6029E-23 SHAR_EOF fi # end of overwriting check if test -f 'res3' then echo shar: will not over-write existing file "'res3'" else cat << SHAR_EOF > 'res3' n m integral rel error err const 1 2 0.7546452E+01 0.1642E-01 0.1204E-01 om= 0.100 2 4 0.7667665E+01 0.6252E-03 0.4843E-04 3 6 0.7672359E+01 0.1355E-04 0.8441E-07 4 8 0.7672458E+01 0.5593E-06 0.8250E-10 5 10 0.7672466E+01 0.4350E-06 0.5157E-13 6 12 0.7672466E+01 0.4972E-06 0.2238E-16 7 14 0.7672463E+01 0.1243E-06 0.7136E-20 8 16 0.7672465E+01 0.3107E-06 0.1742E-23 9 18 0.7672470E+01 0.9944E-06 0.3360E-27 10 20 0.7672466E+01 0.4350E-06 0.5249E-31 1 2 0.2485091E+01 0.2665E-01 0.3101E-01 om= 0.500 2 4 0.2550596E+01 0.9954E-03 0.1217E-03 3 6 0.2553083E+01 0.2129E-04 0.2113E-06 4 8 0.2553135E+01 0.8404E-06 0.2062E-09 5 10 0.2553137E+01 0.0000E+00 0.1289E-12 6 12 0.2553138E+01 0.2801E-06 0.5591E-16 7 14 0.2553138E+01 0.1868E-06 0.1782E-19 8 16 0.2553138E+01 0.1868E-06 0.4350E-23 9 18 0.2553138E+01 0.1868E-06 0.8390E-27 10 20 0.2553139E+01 0.5603E-06 0.1311E-30 1 2 0.1433538E+01 0.3039E-01 0.4022E-01 om= 1.000 2 4 0.1476804E+01 0.1125E-02 0.1563E-03 3 6 0.1478432E+01 0.2371E-04 0.2709E-06 4 8 0.1478467E+01 0.4032E-06 0.2643E-09 5 10 0.1478468E+01 0.2419E-06 0.1651E-12 6 12 0.1478468E+01 0.4838E-06 0.7162E-16 7 14 0.1478468E+01 0.3225E-06 0.2283E-19 8 16 0.1478468E+01 0.2419E-06 0.5572E-23 9 18 0.1478467E+01 0.8063E-07 0.1075E-26 10 20 0.1478467E+01 0.1613E-06 0.1679E-30 1 2 0.7906218E+00 0.3313E-01 0.4783E-01 om= 1.999 2 4 0.8167189E+00 0.1217E-02 0.1844E-03 3 6 0.8176936E+00 0.2515E-04 0.3193E-06 4 8 0.8177140E+00 0.2187E-06 0.3114E-09 5 10 0.8177141E+00 0.7289E-07 0.1944E-12 6 12 0.8177143E+00 0.2187E-06 0.8434E-16 7 14 0.8177143E+00 0.1458E-06 0.2688E-19 8 16 0.8177145E+00 0.3645E-06 0.6561E-23 9 18 0.8177142E+00 0.0000E+00 0.1265E-26 10 20 0.8177140E+00 0.1458E-06 0.1977E-30 1 2 0.7902704E+00 0.3313E-01 0.4783E-01 om= 2.000 2 4 0.8163572E+00 0.1217E-02 0.1844E-03 3 6 0.8173314E+00 0.2538E-04 0.3193E-06 4 8 0.8173519E+00 0.2917E-06 0.3114E-09 5 10 0.8173526E+00 0.5105E-06 0.1945E-12 6 12 0.8173519E+00 0.2917E-06 0.8435E-16 7 14 0.8173525E+00 0.3646E-06 0.2689E-19 8 16 0.8173525E+00 0.4375E-06 0.6562E-23 9 18 0.8173521E+00 0.7292E-07 0.1265E-26 10 20 0.8173527E+00 0.5834E-06 0.1977E-30 1 2 0.5478238E+00 0.3427E-01 0.5121E-01 om= 3.000 2 4 0.5665526E+00 0.1255E-02 0.1968E-03 3 6 0.5672499E+00 0.2564E-04 0.3405E-06 4 8 0.5672644E+00 0.0000E+00 0.3320E-09 5 10 0.5672650E+00 0.1051E-05 0.2073E-12 6 12 0.5672646E+00 0.3152E-06 0.8993E-16 7 14 0.5672645E+00 0.1051E-06 0.2867E-19 8 16 0.5672647E+00 0.4203E-06 0.6996E-23 9 18 0.5672642E+00 0.4203E-06 0.1349E-26 10 20 0.5672643E+00 0.2101E-06 0.2108E-30 1 2 0.5476564E+00 0.3427E-01 0.5121E-01 om= 3.001 2 4 0.5663798E+00 0.1255E-02 0.1968E-03 3 6 0.5670768E+00 0.2586E-04 0.3405E-06 4 8 0.5670917E+00 0.4204E-06 0.3321E-09 5 10 0.5670917E+00 0.4204E-06 0.2073E-12 6 12 0.5670913E+00 0.3153E-06 0.8994E-16 7 14 0.5670916E+00 0.3153E-06 0.2867E-19 8 16 0.5670917E+00 0.4204E-06 0.6996E-23 9 18 0.5670916E+00 0.3153E-06 0.1349E-26 10 20 0.5670913E+00 0.2102E-06 0.2108E-30 1 2 0.1668670E+00 0.3618E-01 0.5719E-01 om=10.500 2 4 0.1729025E+00 0.1317E-02 0.2186E-03 3 6 0.1731258E+00 0.2728E-04 0.3778E-06 4 8 0.1731305E+00 0.4303E-06 0.3682E-09 5 10 0.1731307E+00 0.5164E-06 0.2299E-12 6 12 0.1731306E+00 0.3443E-06 0.9971E-16 7 14 0.1731307E+00 0.7746E-06 0.3178E-19 8 16 0.1731306E+00 0.1721E-06 0.7756E-23 9 18 0.1731305E+00 0.1721E-06 0.1496E-26 10 20 0.1731305E+00 0.3443E-06 0.2336E-30 1 2 0.3584975E-01 0.3687E-01 0.5945E-01 om=50.000 2 4 0.3717221E-01 0.1340E-02 0.2267E-03 3 6 0.3722105E-01 0.2772E-04 0.3917E-06 4 8 0.3722207E-01 0.4003E-06 0.3818E-09 5 10 0.3722209E-01 0.2002E-06 0.2384E-12 6 12 0.3722208E-01 0.2002E-06 0.1034E-15 7 14 0.3722208E-01 0.0000E+00 0.3295E-19 8 16 0.3722209E-01 0.3002E-06 0.8041E-23 9 18 0.3722208E-01 0.0000E+00 0.1551E-26 10 20 0.3722207E-01 0.5004E-06 0.2422E-30 SHAR_EOF fi # end of overwriting check if test -f 'res4' then echo shar: will not over-write existing file "'res4'" else cat << SHAR_EOF > 'res4' n m integral rel error err const 1 2 0.9650458E+02 0.2059E-02 0.2025E-03 2 4 0.9669596E+02 0.7992E-04 0.8408E-06 3 6 0.9670349E+02 0.2051E-05 0.1472E-08 4 8 0.9670368E+02 0.7889E-07 0.1441E-11 5 10 0.9670367E+02 0.2367E-06 0.9012E-15 6 12 0.9670368E+02 0.7889E-07 0.3913E-18 7 14 0.9670372E+02 0.3156E-06 0.1248E-21 8 16 0.9670367E+02 0.1578E-06 0.3046E-25 9 18 0.9670367E+02 0.2367E-06 0.5875E-29 10 20 0.9670371E+02 0.2367E-06 0.9180E-33 SHAR_EOF fi # end of overwriting check if test -f 'res5' then echo shar: will not over-write existing file "'res5'" else cat << SHAR_EOF > 'res5' n m integral rel error err const 1 2 0.2754090E+00 0.5199E-01 0.2782E+00 2 4 0.2901222E+00 0.1343E-02 0.5457E-01 3 6 0.2905045E+00 0.2739E-04 0.9089E-02 4 8 0.2905124E+00 0.1026E-06 0.1394E-02 5 10 0.2905125E+00 0.3078E-06 0.2033E-03 6 12 0.2905124E+00 0.0000E+00 0.2868E-04 7 14 0.2905126E+00 0.6155E-06 0.3950E-05 8 16 0.2905124E+00 0.0000E+00 0.5343E-06 9 18 0.2905126E+00 0.7181E-06 0.7128E-07 10 20 0.2905125E+00 0.2052E-06 0.9403E-08 SHAR_EOF fi # end of overwriting check if test -f 'res6' then echo shar: will not over-write existing file "'res6'" else cat << SHAR_EOF > 'res6' n m integral rel error err const 1 1 0.2843093E+00 0.2512E+00 0.1600E+00 2 3 0.3767365E+00 0.7828E-02 0.1738E-01 3 5 0.3796442E+00 0.1702E-03 0.2003E-02 4 7 0.3797079E+00 0.2512E-05 0.2344E-03 5 9 0.3797088E+00 0.1570E-06 0.2762E-04 6 11 0.3797089E+00 0.0000E+00 0.3265E-05 7 13 0.3797089E+00 0.2355E-06 0.3869E-06 8 15 0.3797089E+00 0.2355E-06 0.4592E-07 9 17 0.3797090E+00 0.3924E-06 0.5454E-08 10 19 0.3797091E+00 0.5494E-06 0.6482E-09 SHAR_EOF fi # end of overwriting check if test -f 'res7' then echo shar: will not over-write existing file "'res7'" else cat << SHAR_EOF > 'res7' n m integral rel error err const 1 1 0.1690406E+01 0.2376E+00 0.4276E-03 2 3 0.2201195E+01 0.7196E-02 0.3995E-04 3 5 0.2216809E+01 0.1540E-03 0.4279E-05 4 7 0.2217144E+01 0.2581E-05 0.4785E-06 5 9 0.2217151E+01 0.2151E-06 0.5462E-07 6 11 0.2217149E+01 0.3226E-06 0.6306E-08 7 13 0.2217151E+01 0.3226E-06 0.7332E-09 8 15 0.2217151E+01 0.4301E-06 0.8568E-10 9 17 0.2217151E+01 0.3226E-06 0.1005E-10 10 19 0.2217150E+01 0.0000E+00 0.1181E-11 SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << SHAR_EOF > 'src.f' c c SUBROUTINE DABMOD(N,NCAPMD,MCAP,DEPS,IROUT,DA,DB,DXIR,DXII,IS, + DALPHA,DBETA,NCAPD,KOUNTD,IERRD,DBE,DX,DW,DE, + DP0,DP1,DP2) c double precision theta c common/par/theta c .. Scalar Arguments .. DOUBLE PRECISION DEPS INTEGER IERRD,IROUT,KOUNTD,MCAP,N,NCAPD,NCAPMD c .. c .. Array Arguments .. The upper bounds DXII, DXIR and IS are all c MCAP, but this doesn't work in Fortran 77 c if MCAP = 0. DOUBLE PRECISION DA(NCAPMD),DALPHA(N),DB(NCAPMD),DBE(N),DBETA(N), + DE(NCAPMD),DP0(NCAPMD),DP1(NCAPMD),DP2(NCAPMD), + DW(NCAPMD),DX(NCAPMD),DXII(*),DXIR(*) INTEGER IS(*) c .. c .. Local Scalars .. DOUBLE PRECISION DEPSM,DP INTEGER IED,IERRGD,IMU,INCAP,K,MU c .. c .. External Functions .. DOUBLE PRECISION D1MACH EXTERNAL D1MACH c .. c .. External Subroutines .. EXTERNAL DGAUSS,DLANCZ,DSTI c .. c .. Intrinsic Functions .. INTRINSIC ABS c .. integer i1mach, nout nout = i1mach(2) IERRD = 0 DEPSM = D1MACH(3) INCAP = 1 DO 10 K = 1,N DALPHA(K) = DA(K) DBETA(K) = DB(K) 10 CONTINUE IF (MCAP.GT.0) THEN KOUNTD = -1 DO 20 K = 1,N DBETA(K) = 0.D0 20 CONTINUE NCAPD = (2*N-1)/2 30 DO 40 K = 1,N DBE(K) = DBETA(K) 40 CONTINUE KOUNTD = KOUNTD + 1 IF (KOUNTD.GT.1) INCAP = 2** (KOUNTD/5)*N NCAPD = NCAPD + INCAP IF (NCAPD.GT.NCAPMD) THEN IERRD = NCAPMD RETURN END IF CALL DGAUSS(NCAPD,DA,DB,DEPSM,DX,DW,IERRGD,DE) IF (IERRGD.NE.0) THEN WRITE (NOUT,FMT=9000) IERRGD IERRD = 1 RETURN END IF DO 60 K = 1,NCAPD DP = 1.D0 IMU = 0 DO 50 MU = 1,MCAP IF (IMU.EQ.0) THEN IF (DXII(MU).EQ.0.D0) THEN DP = ((1.D0+DXIR(MU)*DX(K))**IS(MU))*DP ELSE DP = (((1.D0+DXIR(MU)*DX(K))**2+ + (DXII(MU)*DX(K))**2)**IS(MU))*DP IMU = 1 END IF ELSE IMU = 0 END IF 50 CONTINUE DW(K) = DW(K)/DP c dw(k)=dsqrt(1.d0+.5d0*theta*dx(k))*dw(k)/dp 60 CONTINUE IF (IROUT.EQ.1) THEN CALL DSTI(N,NCAPD,DX,DW,DALPHA,DBETA,IED,DP0,DP1,DP2) IF (IED.NE.0) THEN IERRD = 2 RETURN END IF ELSE CALL DLANCZ(N,NCAPD,DX,DW,DALPHA,DBETA,IED,DP0,DP1) IF (IED.NE.0) THEN IERRD = 2 RETURN END IF END IF DO 70 K = 1,N IF (ABS(DBETA(K)-DBE(K)).GT.DEPS*DBETA(K)) GO TO 30 70 CONTINUE END IF RETURN 9000 FORMAT (1X,'ierrgd in dgauss=',I5) END c c SUBROUTINE DGAUSS(N,DALPHA,DBETA,DEPS,DZERO,DWEIGH,IERR,DE) c c This is a double-precision version of the routine gauss. c C .. Scalar Arguments .. DOUBLE PRECISION DEPS INTEGER IERR,N C .. C .. Array Arguments .. DOUBLE PRECISION DALPHA(N),DBETA(N),DE(N),DWEIGH(N),DZERO(N) C .. C .. Local Scalars .. DOUBLE PRECISION DB,DC,DF,DG,DP,DR,DS INTEGER I,II,J,K,L,M,MML C .. C .. Intrinsic Functions .. INTRINSIC ABS,SIGN,SQRT C .. IF (N.LT.1) THEN IERR = -1 RETURN END IF IERR = 0 DZERO(1) = DALPHA(1) IF (DBETA(1).LT.0.D0) THEN IERR = -2 RETURN END IF DWEIGH(1) = DBETA(1) IF (N.EQ.1) RETURN DWEIGH(1) = 1.D0 DE(N) = 0.D0 DO 10 K = 2,N DZERO(K) = DALPHA(K) IF (DBETA(K).LT.0.D0) THEN IERR = -2 RETURN END IF DE(K-1) = SQRT(DBETA(K)) DWEIGH(K) = 0.D0 10 CONTINUE DO 80 L = 1,N J = 0 20 DO 30 M = L,N IF (M.EQ.N) GO TO 40 IF (ABS(DE(M)).LE.DEPS* (ABS(DZERO(M))+ABS(DZERO(M+ + 1)))) GO TO 40 30 CONTINUE 40 DP = DZERO(L) IF (M.EQ.L) GO TO 80 IF (J.EQ.30) GO TO 120 J = J + 1 DG = (DZERO(L+1)-DP)/ (2.D0*DE(L)) DR = SQRT(DG*DG+1.D0) DG = DZERO(M) - DP + DE(L)/ (DG+SIGN(DR,DG)) DS = 1.D0 DC = 1.D0 DP = 0.D0 MML = M - L DO 70 II = 1,MML I = M - II DF = DS*DE(I) DB = DC*DE(I) IF (ABS(DF).LT.ABS(DG)) GO TO 50 DC = DG/DF DR = SQRT(DC*DC+1.D0) DE(I+1) = DF*DR DS = 1.D0/DR DC = DC*DS GO TO 60 50 DS = DF/DG DR = SQRT(DS*DS+1.D0) DE(I+1) = DG*DR DC = 1.D0/DR DS = DS*DC 60 DG = DZERO(I+1) - DP DR = (DZERO(I)-DG)*DS + 2.D0*DC*DB DP = DS*DR DZERO(I+1) = DG + DP DG = DC*DR - DB DF = DWEIGH(I+1) DWEIGH(I+1) = DS*DWEIGH(I) + DC*DF DWEIGH(I) = DC*DWEIGH(I) - DS*DF 70 CONTINUE DZERO(L) = DZERO(L) - DP DE(L) = DG DE(M) = 0.D0 GO TO 20 80 CONTINUE DO 100 II = 2,N I = II - 1 K = I DP = DZERO(I) DO 90 J = II,N IF (DZERO(J).GE.DP) GO TO 90 K = J DP = DZERO(J) 90 CONTINUE IF (K.EQ.I) GO TO 100 DZERO(K) = DZERO(I) DZERO(I) = DP DP = DWEIGH(I) DWEIGH(I) = DWEIGH(K) DWEIGH(K) = DP 100 CONTINUE DO 110 K = 1,N DWEIGH(K) = DBETA(1)*DWEIGH(K)*DWEIGH(K) 110 CONTINUE RETURN 120 IERR = L RETURN END c c SUBROUTINE DGCHRS(N,IOPT,DA,DB,DX,DHP,DHN,DALPHA,DBETA,IERRD) C .. Scalar Arguments .. DOUBLE PRECISION DHN,DHP,DX INTEGER IERRD,IOPT,N C .. C .. Array Arguments .. DOUBLE PRECISION DA(N),DALPHA(N),DB(N),DBETA(N) C .. C .. Local Scalars .. DOUBLE PRECISION DE,DEH,DQ,DQH INTEGER K C .. IERRD = 0 IF (IOPT.EQ.1) THEN DALPHA(1) = DX - DB(1)/DHP DBETA(1) = -DHP DQ = -DB(1)/DHP DO 10 K = 2,N DE = DA(K-1) - DX - DQ DBETA(K) = DQ*DE DQ = DB(K)/DE DALPHA(K) = DQ + DE + DX 10 CONTINUE RETURN ELSE IF (IOPT.EQ.2) THEN DALPHA(1) = DX* (DHP+DHN)/ (DHP-DHN) DBETA(1) = - (DHP-DHN)/ (2.D0*DX) DQ = -DB(1)/DHP DQH = -DHP/DBETA(1) DE = 0.D0 DO 20 K = 2,N DEH = DQ + DE + 2.D0*DX - DQH DBETA(K) = DQH*DEH DE = DA(K-1) - DX - DQ DQH = DQ*DE/DEH DALPHA(K) = DQH + DEH - DX DQ = DB(K)/DE 20 CONTINUE RETURN ELSE IERRD = 1 RETURN END IF END c c SUBROUTINE DGQRAT(N,MCAP,DALPHA,DBETA,DXIR,DXII,IS,DZG,DWG,DCONST, + IERRD,DE) c .. Scalar Arguments .. DOUBLE PRECISION DCONST INTEGER IERRD,MCAP,N c .. c .. Array Arguments .. The upper bounds DXII, DXIR and IS are all c MCAP, but this doesn't work in Fortran 77 c if MCAP = 0. DOUBLE PRECISION DALPHA(N+1),DBETA(N+1),DE(N),DWG(N),DXII(*), + DXIR(*),DZG(N) INTEGER IS(*) c .. c .. Local Scalars .. DOUBLE PRECISION DEPSM,DP INTEGER IERRGD,IMU,K,MU c .. c .. External Functions .. DOUBLE PRECISION D1MACH EXTERNAL D1MACH c .. c .. External Subroutines .. EXTERNAL DGAUSS c .. c .. Intrinsic Functions .. INTRINSIC DBLE c .. integer i1mach, nout nout = i1mach(2) IERRD = 0 DEPSM = D1MACH(3) CALL DGAUSS(N,DALPHA,DBETA,DEPSM,DZG,DWG,IERRGD,DE) IF (IERRGD.NE.0) THEN WRITE (nout,FMT=9000) IERRGD IERRD = 1 RETURN END IF DCONST = DBETA(1) DO 20 K = 1,N DCONST = DBETA(K+1)*DCONST/DBLE(2*K* (2*K-1)) IF (MCAP.GT.0) THEN DP = 1.D0 IMU = 0 DO 10 MU = 1,MCAP IF (IMU.EQ.0) THEN IF (DXII(MU).EQ.0.D0) THEN DP = ((1.D0+DXIR(MU)*DZG(K))**IS(MU))*DP ELSE DP = (((1.D0+DXIR(MU)*DZG(K))**2+ + (DXII(MU)*DZG(K))**2)**IS(MU))*DP IMU = 1 END IF ELSE IMU = 0 END IF 10 CONTINUE DWG(K) = DP*DWG(K) END IF 20 CONTINUE RETURN 9000 FORMAT (1X,'ierrgd in dgauss=',I5) END c c SUBROUTINE DKNUM(N,NU0,NUMAX,DX,DY,DEPS,DA,DB,DRHOR,DRHOI,NU,IERR, + DROLDR,DROLDI) c c This is a double-precision version of the routine knum. c c c The arrays drhor,drhoi,droldr,droldi are assumed to have c dimension n+1. c C .. Scalar Arguments .. DOUBLE PRECISION DEPS,DX,DY INTEGER IERR,N,NU,NU0,NUMAX C .. C .. Array Arguments .. DOUBLE PRECISION DA(NUMAX),DB(NUMAX),DRHOI(*),DRHOR(*),DROLDI(*), + DROLDR(*) C .. C .. Local Scalars .. DOUBLE PRECISION DDEN,DRI,DRR,DT INTEGER J,J1,K,NP1 C .. IERR = 0 NP1 = N + 1 IF (NU0.GT.NUMAX) THEN IERR = NU0 RETURN END IF IF (NU0.LT.NP1) NU0 = NP1 NU = NU0 - 5 DO 10 K = 1,NP1 DRHOR(K) = 0.D0 DRHOI(K) = 0.D0 10 CONTINUE 20 NU = NU + 5 IF (NU.GT.NUMAX) THEN IERR = NUMAX GO TO 60 END IF DO 30 K = 1,NP1 DROLDR(K) = DRHOR(K) DROLDI(K) = DRHOI(K) 30 CONTINUE DRR = 0.D0 DRI = 0.D0 DO 40 J = 1,NU J1 = NU - J + 1 DDEN = (DX-DA(J1)-DRR)**2 + (DY-DRI)**2 DRR = DB(J1)* (DX-DA(J1)-DRR)/DDEN DRI = -DB(J1)* (DY-DRI)/DDEN IF (J1.LE.NP1) THEN DRHOR(J1) = DRR DRHOI(J1) = DRI END IF 40 CONTINUE DO 50 K = 1,NP1 IF ((DRHOR(K)-DROLDR(K))**2+ (DRHOI(K)-DROLDI(K))**2.GT. + DEPS* (DRHOR(K)**2+DRHOI(K)**2)) GO TO 20 50 CONTINUE 60 IF (N.EQ.0) RETURN DO 70 K = 2,NP1 DT = DRHOR(K)*DRHOR(K-1) - DRHOI(K)*DRHOI(K-1) DRHOI(K) = DRHOR(K)*DRHOI(K-1) + DRHOI(K)*DRHOR(K-1) DRHOR(K) = DT 70 CONTINUE RETURN END c c SUBROUTINE DLANCZ(N,NCAP,DX,DW,DALPHA,DBETA,IERR,DP0,DP1) c c This is a double-precision version of the routine lancz. c C .. Scalar Arguments .. INTEGER IERR,N,NCAP C .. C .. Array Arguments .. DOUBLE PRECISION DALPHA(N),DBETA(N),DP0(NCAP),DP1(NCAP),DW(NCAP), + DX(NCAP) C .. C .. Local Scalars .. DOUBLE PRECISION DGAM,DPI,DRHO,DSIG,DT,DTK,DTMP,DTSIG,DXLAM INTEGER I,K C .. IF (N.LE.0 .OR. N.GT.NCAP) THEN IERR = 1 RETURN ELSE IERR = 0 END IF DO 10 I = 1,NCAP DP0(I) = DX(I) DP1(I) = 0.D0 10 CONTINUE DP1(1) = DW(1) DO 30 I = 1,NCAP - 1 DPI = DW(I+1) DGAM = 1.D0 DSIG = 0.D0 DT = 0.D0 DXLAM = DX(I+1) DO 20 K = 1,I + 1 DRHO = DP1(K) + DPI DTMP = DGAM*DRHO DTSIG = DSIG IF (DRHO.LE.0.D0) THEN DGAM = 1.D0 DSIG = 0.D0 ELSE DGAM = DP1(K)/DRHO DSIG = DPI/DRHO END IF DTK = DSIG* (DP0(K)-DXLAM) - DGAM*DT DP0(K) = DP0(K) - (DTK-DT) DT = DTK IF (DSIG.LE.0.D0) THEN DPI = DTSIG*DP1(K) ELSE DPI = (DT**2)/DSIG END IF DTSIG = DSIG DP1(K) = DTMP 20 CONTINUE 30 CONTINUE DO 40 K = 1,N DALPHA(K) = DP0(K) DBETA(K) = DP1(K) 40 CONTINUE RETURN END c c SUBROUTINE DRECUR(N,IPOLY,DAL,DBE,DA,DB,IDERR) c c This is a double-precision version of the routine recur. c C .. Scalar Arguments .. DOUBLE PRECISION DAL,DBE INTEGER IDERR,IPOLY,N C .. C .. Array Arguments .. DOUBLE PRECISION DA(N),DB(N) C .. C .. Local Scalars .. DOUBLE PRECISION DAL2,DALPBE,DBE2,DKM1,DLMACH,DT INTEGER K C .. C .. External Functions .. DOUBLE PRECISION D1MACH,DGAMMA,DLGA EXTERNAL D1MACH,DGAMMA,DLGA C .. C .. Intrinsic Functions .. INTRINSIC ATAN,DBLE,EXP,LOG,SQRT C .. IF (N.LT.1) THEN IDERR = 3 RETURN END IF DLMACH = LOG(D1MACH(2)) IDERR = 0 DO 10 K = 1,N DA(K) = 0.D0 10 CONTINUE IF (IPOLY.EQ.1) THEN DB(1) = 2.D0 IF (N.EQ.1) RETURN DO 20 K = 2,N DKM1 = DBLE(K-1) DB(K) = 1.D0/ (4.D0-1.D0/ (DKM1*DKM1)) 20 CONTINUE RETURN ELSE IF (IPOLY.EQ.2) THEN DA(1) = .5D0 DB(1) = 1.D0 IF (N.EQ.1) RETURN DO 30 K = 2,N DA(K) = .5D0 DKM1 = DBLE(K-1) DB(K) = .25D0/ (4.D0-1.D0/ (DKM1*DKM1)) 30 CONTINUE RETURN ELSE IF (IPOLY.EQ.3) THEN DB(1) = 4.D0*ATAN(1.D0) IF (N.EQ.1) RETURN DB(2) = .5D0 IF (N.EQ.2) RETURN DO 40 K = 3,N DB(K) = .25D0 40 CONTINUE RETURN ELSE IF (IPOLY.EQ.4) THEN DB(1) = 2.D0*ATAN(1.D0) IF (N.EQ.1) RETURN DO 50 K = 2,N DB(K) = .25D0 50 CONTINUE RETURN ELSE IF (IPOLY.EQ.5) THEN DB(1) = 4.D0*ATAN(1.D0) DA(1) = .5D0 IF (N.EQ.1) RETURN DO 60 K = 2,N DB(K) = .25D0 60 CONTINUE RETURN ELSE IF (IPOLY.EQ.6) THEN IF (DAL.LE.-1.D0 .OR. DBE.LE.-1.D0) THEN IDERR = 1 RETURN ELSE DALPBE = DAL + DBE DA(1) = (DBE-DAL)/ (DALPBE+2.D0) DT = (DALPBE+1.D0)*LOG(2.D0) + DLGA(DAL+1.D0) + + DLGA(DBE+1.D0) - DLGA(DALPBE+2.D0) IF (DT.GT.DLMACH) THEN IDERR = 2 DB(1) = D1MACH(2) ELSE DB(1) = EXP(DT) END IF IF (N.EQ.1) RETURN DAL2 = DAL*DAL DBE2 = DBE*DBE DA(2) = (DBE2-DAL2)/ ((DALPBE+2.D0)* (DALPBE+4.D0)) DB(2) = 4.D0* (DAL+1.D0)* (DBE+1.D0)/ + ((DALPBE+3.D0)* (DALPBE+2.D0)**2) IF (N.EQ.2) RETURN DO 70 K = 3,N DKM1 = DBLE(K-1) DA(K) = .25D0* (DBE2-DAL2)/ + (DKM1*DKM1* (1.D0+.5D0*DALPBE/DKM1)* + (1.D0+.5D0* (DALPBE+2.D0)/DKM1)) DB(K) = .25D0* (1.D0+DAL/DKM1)* (1.D0+DBE/DKM1)* + (1.D0+DALPBE/DKM1)/ ((1.D0+.5D0* (DALPBE+ + 1.D0)/DKM1)* (1.D0+.5D0* (DALPBE-1.D0)/DKM1)* + (1.D0+.5D0*DALPBE/DKM1)**2) 70 CONTINUE RETURN END IF ELSE IF (IPOLY.EQ.7) THEN IF (DAL.LE.-1.D0) THEN IDERR = 1 RETURN ELSE DA(1) = DAL + 1.D0 DB(1) = DGAMMA(DAL+1.D0,IDERR) IF (IDERR.EQ.2) DB(1) = D1MACH(2) IF (N.EQ.1) RETURN DO 80 K = 2,N DKM1 = DBLE(K-1) DA(K) = 2.D0*DKM1 + DAL + 1.D0 DB(K) = DKM1* (DKM1+DAL) 80 CONTINUE RETURN END IF ELSE IF (IPOLY.EQ.8) THEN DB(1) = SQRT(4.D0*ATAN(1.D0)) IF (N.EQ.1) RETURN DO 90 K = 2,N DB(K) = .5D0*DBLE(K-1) 90 CONTINUE RETURN ELSE IDERR = 4 END IF END DOUBLE PRECISION FUNCTION DLGA(DX) C .. Scalar Arguments .. DOUBLE PRECISION DX C .. C .. Local Scalars .. DOUBLE PRECISION DC,DP,DS,DT,DY REAL DPREC,Y,Y0 INTEGER I C .. C .. Local Arrays .. DOUBLE PRECISION DBDEN(8),DBNUM(8) C .. C .. External Functions .. DOUBLE PRECISION D1MACH EXTERNAL D1MACH C .. C .. Intrinsic Functions .. INTRINSIC ATAN,EXP,LOG,LOG10,REAL C .. C .. Data statements .. c c This routine evaluates the logarithm of the gamma function by a c combination of recurrence and asymptotic approximation. c c The entries in the next data statement are the numerators and c denominators, respectively, of the quantities B[16]/(16*15), c B[14]/(14*13),..., B[2]/(2*1), where B[2n] are the Bernoulli c numbers. c DATA DBNUM/-3.617D3,1.D0,-6.91D2,1.D0,-1.D0,1.D0,-1.D0,1.D0/, + DBDEN/1.224D5,1.56D2,3.6036D5,1.188D3,1.68D3,1.26D3,3.6D2, + 1.2D1/ C .. c c The quantity dprec in the next statement is the number of decimal c digits carried in double-precision floating-point arithmetic. c DPREC = -LOG10(REAL(D1MACH(3))) DC = .5D0*LOG(8.D0*ATAN(1.D0)) DP = 1.D0 DY = DX Y = REAL(DY) c c The quantity y0 below is the threshold value beyond which asymptotic c evaluation gives sufficient accuracy; see Eq. 6.1.42 in M. Abramowitz c and I.A. Stegun,Handbook of Mathematical Functions''. The constants c are .12118868... = ln(10)/19 and .05390522... = ln(|B[20]|/190)/19. c Y0 = EXP(.121189*DPREC+.053905) 10 IF (Y.GT.Y0) GO TO 20 DP = DY*DP DY = DY + 1.D0 Y = REAL(DY) GO TO 10 20 DT = 1.D0/ (DY*DY) c c The right-hand side of the next assignment statement is B[18]/(18*17). c DS = 4.3867D4/2.44188D5 DO 30 I = 1,8 DS = DT*DS + DBNUM(I)/DBDEN(I) 30 CONTINUE DLGA = (DY-.5D0)*LOG(DY) - DY + DC + DS/DY - LOG(DP) RETURN END DOUBLE PRECISION FUNCTION DGAMMA(DX,IDERR) c c This evaluates the gamma function for real positive dx, using the c function subroutine dlga. c C .. Scalar Arguments .. DOUBLE PRECISION DX INTEGER IDERR C .. C .. Local Scalars .. DOUBLE PRECISION DLMACH,DT C .. C .. External Functions .. DOUBLE PRECISION D1MACH,DLGA EXTERNAL D1MACH,DLGA C .. C .. Intrinsic Functions .. INTRINSIC EXP,LOG C .. DLMACH = LOG(D1MACH(2)) IDERR = 0 DT = DLGA(DX) IF (DT.GE.DLMACH) THEN IDERR = 2 DGAMMA = D1MACH(2) RETURN ELSE DGAMMA = EXP(DT) RETURN END IF END c c SUBROUTINE DSTI(N,NCAP,DX,DW,DALPHA,DBETA,IERR,DP0,DP1,DP2) c c This is a double-precision version of the routine sti. c C .. Scalar Arguments .. INTEGER IERR,N,NCAP C .. C .. Array Arguments .. DOUBLE PRECISION DALPHA(N),DBETA(N),DP0(NCAP),DP1(NCAP),DP2(NCAP), + DW(NCAP),DX(NCAP) C .. C .. Local Scalars .. DOUBLE PRECISION DHUGE,DSUM0,DSUM1,DSUM2,DT,DTINY INTEGER K,M,NM1 C .. C .. External Functions .. DOUBLE PRECISION D1MACH EXTERNAL D1MACH C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. DTINY = 10.D0*D1MACH(1) DHUGE = .1D0*D1MACH(2) IERR = 0 IF (N.LE.0 .OR. N.GT.NCAP) THEN IERR = 1 RETURN END IF NM1 = N - 1 DSUM0 = 0.D0 DSUM1 = 0.D0 DO 10 M = 1,NCAP DSUM0 = DSUM0 + DW(M) DSUM1 = DSUM1 + DW(M)*DX(M) 10 CONTINUE DALPHA(1) = DSUM1/DSUM0 DBETA(1) = DSUM0 IF (N.EQ.1) RETURN DO 20 M = 1,NCAP DP1(M) = 0.D0 DP2(M) = 1.D0 20 CONTINUE DO 40 K = 1,NM1 DSUM1 = 0.D0 DSUM2 = 0.D0 DO 30 M = 1,NCAP IF (DW(M).EQ.0.D0) GO TO 30 DP0(M) = DP1(M) DP1(M) = DP2(M) DP2(M) = (DX(M)-DALPHA(K))*DP1(M) - DBETA(K)*DP0(M) IF (ABS(DP2(M)).GT.DHUGE .OR. ABS(DSUM2).GT.DHUGE) THEN IERR = K RETURN END IF DT = DW(M)*DP2(M)*DP2(M) DSUM1 = DSUM1 + DT DSUM2 = DSUM2 + DT*DX(M) 30 CONTINUE IF (ABS(DSUM1).LT.DTINY) THEN IERR = -K RETURN END IF DALPHA(K+1) = DSUM2/DSUM1 DBETA(K+1) = DSUM1/DSUM0 DSUM0 = DSUM1 40 CONTINUE RETURN END c c INTEGER FUNCTION NU0HER(N,Z,EPS) c c This is an auxiliary function routine providing a starting backward c recurrence index for the Hermite measure that can be used in place c of nu0 in the routines knum and dknum. c C .. Scalar Arguments .. COMPLEX Z REAL EPS INTEGER N C .. C .. Intrinsic Functions .. INTRINSIC ABS,AIMAG,LOG,REAL,SQRT C .. NU0HER = 2.* (SQRT(.5*REAL(N+1))+.25*LOG(1./EPS)/ABS(AIMAG(Z)))**2 RETURN END c c INTEGER FUNCTION NU0JAC(N,Z,EPS) c c This is an auxiliary function routine providing a starting backward c recurrence index for the Jacobi measure that can be used in place c of nu0 in the routines knum and dknum. c C .. Scalar Arguments .. COMPLEX Z REAL EPS INTEGER N C .. C .. Local Scalars .. REAL ANGLE,PI,R,X,X2,Y,Y2 C .. C .. Intrinsic Functions .. INTRINSIC ABS,AIMAG,ATAN,COS,LOG,REAL,SIN,SQRT C .. PI = 4.*ATAN(1.) X = REAL(Z) Y = ABS(AIMAG(Z)) IF (X.LT.1.) THEN IF (X.LT.-1.) ANGLE = .5* (2.*PI+ATAN(Y/ (X-1.))+ + ATAN(Y/ (X+1.))) IF (X.EQ.-1.) ANGLE = .5* (1.5*PI-ATAN(.5*Y)) IF (X.GT.-1.) ANGLE = .5* (PI+ATAN(Y/ (X-1.))+ATAN(Y/ (X+1.))) ELSE IF (X.EQ.1.) ANGLE = .5* (.5*PI+ATAN(.5*Y)) IF (X.GT.1.) ANGLE = .5* (ATAN(Y/ (X-1.))+ATAN(Y/ (X+1.))) END IF X2 = X*X Y2 = Y*Y R = ((X2-Y2-1.)**2+4.*X2*Y2)**.25 R = SQRT((X+R*COS(ANGLE))**2+ (Y+R*SIN(ANGLE))**2) NU0JAC = REAL(N+1) + .5*LOG(1./EPS)/LOG(R) RETURN END c c INTEGER FUNCTION NU0LAG(N,Z,AL,EPS) c c This is an auxiliary function routine providing a starting backward c recurrence index for the Laguerre measure that can be used in place c of nu0 in the routines knum and dknum. c C .. Scalar Arguments .. COMPLEX Z REAL AL,EPS INTEGER N C .. C .. Local Scalars .. REAL PHI,PI,X,Y C .. C .. Intrinsic Functions .. INTRINSIC AIMAG,ATAN,COS,LOG,REAL,SQRT C .. PI = 4.*ATAN(1.) X = REAL(Z) Y = AIMAG(Z) PHI = .5*PI IF (Y.LT.0.) PHI = 1.5*PI IF (X.EQ.0.) GO TO 10 PHI = ATAN(Y/X) IF (Y.GT.0. .AND. X.GT.0.) GO TO 10 PHI = PHI + PI IF (X.LT.0.) GO TO 10 PHI = PHI + PI 10 NU0LAG = (SQRT(REAL(N+1)+.5* (AL+1.))+ + LOG(1./EPS)/ (4.* (X*X+Y*Y)**.25*COS(.5* (PHI-PI))))**2 - + .5* (AL+1.) RETURN END SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << SHAR_EOF > 'src.f' c c SUBROUTINE ABMOD(N,NCAPM,MCAP,EPS,IROUT,A,B,XIR,XII,IS,ALPHA,BETA, + NCAP,KOUNT,IERR,BE,X,W,E,P0,P1,P2) c c This routine generates recursion coefficients for the modified c measure c c dlambda(t)/p(t), c c where c c p(t) = product[(1+zeta(mu)*t)**is(mu)], c c the product being extended from mu=1 to mu=mcap. (If mcap=0, c the routine simply returns the recursion coefficients of the measure c dlambda.) This is done by a discretization method based on Gauss c quadrature rules relative to the measure dlambda. It is assumed c that all zeta's are either real or occurring in conjugate complex c pairs which involve consecutive values of mu. c c Input: n - - - the number of recursion coefficients desired; c type integer c ncapm - - an upper bound on the number of terms in the c discrete inner product used in the discretization c method; normally the choice ncapm=400 should be c satisfactory; type integer c mcap - - the number of poles (not counting multiplicities); c type integer c eps - - - a (relative) error tolerance; type real c irout - - an integer selecting the routine for generating c the recursion coefficients of the discretized c inner product; if irout=1, the ORTHPOL routine c sti is selected, otherwise the routine lancz c a,b - - - real arrays of dimension ncapm containing the c recursion coefficients of the measure dlambda c xir,xii - real arrays of dimension mcap containing the c real and imaginary parts, respectively, of c zeta(mu); if zeta(mu) is real, put xii(mu)=0. c is - - - an integer array of dimension mcap containing c the respective multiplicities of the poles c c Output: alpha,beta - arrays of dimension n containing the c recursion coefficients of the modified measure c ncap - - the number of terms in the discrete inner product c that yields convergence; type integer c kount - - the number of iterations required in the c discretization method; type integer c ierr - - an error flag, where c ierr=0 on normal return c ierr=1 indicates an error condition in the c ORTHPOL routine gauss c ierr=2 indicates an error condition in the c ORTHPOL routine sti or lancz, whichever c is used c ierr=ncapm if the method does not converge with c ncap less than or equal to ncapm c c All remaining arrays in the calling sequence are internal working c arrays. c c c The common declaration below is to be used only if the parameter c theta in Fermi-Dirac and Bose-Einstein integrals is incorporated c into a weight function as described at the end of Section 5 of the c companion paper c c common/par/theta C .. Scalar Arguments .. REAL EPS INTEGER IERR,IROUT,KOUNT,MCAP,N,NCAP,NCAPM C .. C .. Array Arguments .. The upper bounds XII, XIR and IS are all c MCAP, but this doesn't work in Fortran 77 c if MCAP = 0. REAL A(NCAPM),ALPHA(N),B(NCAPM),BE(N),BETA(N),E(NCAPM),P0(NCAPM), + P1(NCAPM),P2(NCAPM),W(NCAPM),X(NCAPM),XII(*),XIR(*) INTEGER IS(*) C .. C .. Local Scalars .. REAL EPSM,P INTEGER IE,IERRG,IMU,INCAP,K,MU C .. C .. External Functions .. REAL R1MACH EXTERNAL R1MACH C .. C .. External Subroutines .. EXTERNAL GAUSS,LANCZ,STI C .. C .. Intrinsic Functions .. INTRINSIC ABS integer i1mach, nout nout = i1mach(2) C .. IERR = 0 EPSM = R1MACH(3) INCAP = 1 DO 10 K = 1,N ALPHA(K) = A(K) BETA(K) = B(K) 10 CONTINUE IF (MCAP.GT.0) THEN KOUNT = -1 DO 20 K = 1,N BETA(K) = 0. 20 CONTINUE NCAP = (2*N-1)/2 30 DO 40 K = 1,N BE(K) = BETA(K) 40 CONTINUE KOUNT = KOUNT + 1 IF (KOUNT.GT.1) INCAP = 2** (KOUNT/5)*N NCAP = NCAP + INCAP IF (NCAP.GT.NCAPM) THEN IERR = NCAPM RETURN END IF c c Discretization of the inner product c CALL GAUSS(NCAP,A,B,EPSM,X,W,IERRG,E) IF (IERRG.NE.0) THEN IERR = 1 RETURN END IF DO 60 K = 1,NCAP P = 1. IMU = 0 DO 50 MU = 1,MCAP IF (IMU.EQ.0) THEN IF (XII(MU).EQ.0.) THEN P = ((1.+XIR(MU)*X(K))**IS(MU))*P ELSE P = (((1.+XIR(MU)*X(K))**2+ + (XII(MU)*X(K))**2)**IS(MU))*P IMU = 1 END IF ELSE IMU = 0 END IF 50 CONTINUE W(K) = W(K)/P c c The following statement incorporates the parameter theta c into a weight function as described at the end of Section 5 c of the companion paper. It replaces the statement immediately c preceding it. c c w(k)=sqrt(1.+.5*theta*x(k))*w(k)/p c 60 CONTINUE c c Computation of the desired recursion coefficients c IF (IROUT.EQ.1) THEN CALL STI(N,NCAP,X,W,ALPHA,BETA,IE,P0,P1,P2) IF (IE.NE.0) THEN WRITE (nout,FMT=9000) IE IERR = 2 RETURN END IF ELSE CALL LANCZ(N,NCAP,X,W,ALPHA,BETA,IE,P0,P1) IF (IE.NE.0) THEN IERR = 2 RETURN END IF END IF DO 70 K = 1,N IF (ABS(BETA(K)-BE(K)).GT.EPS*BETA(K)) GO TO 30 70 CONTINUE END IF RETURN 9000 FORMAT (1X,'ie in sti=',I3) END c c SUBROUTINE GAUSS(N,ALPHA,BETA,EPS,ZERO,WEIGHT,IERR,E) c c Given n and a measure dlambda, this routine generates the n-point c Gaussian quadrature formula c c integral over supp(dlambda) of f(x)dlambda(x) c c = sum from k=1 to k=n of w(k)f(x(k)) + R(n;f). c c The nodes are returned as zero(k)=x(k) and the weights as c weight(k)=w(k), k=1,2,...,n. The user has to supply the recursion c coefficients alpha(k), beta(k), k=0,1,2,...,n-1, for the measure c dlambda. The routine computes the nodes as eigenvalues, and the c weights in term of the first component of the respective normalized c eigenvectors of the n-th order Jacobi matrix associated with dlambda. c It uses a translation and adaptation of the algol procedure imtql2, c Numer. Math. 12, 1968, 377-383, by Martin and Wilkinson, as modified c by Dubrulle, Numer. Math. 15, 1970, 450. See also Handbook for c Autom. Comput., vol. 2 - Linear Algebra, pp.241-248, and the eispack c routine imtql2. c c Input: n - - the number of points in the Gaussian quadrature c formula; type integer c alpha,beta - - arrays of dimension n to be filled c with the values of alpha(k-1), beta(k-1), k=1,2, c ...,n c eps - the relative accuracy desired in the nodes c and weights c c Output: zero- array of dimension n containing the Gaussian c nodes (in increasing order) zero(k)=x(k), k=1,2, c ...,n c weight - array of dimension n containing the c Gaussian weights weight(k)=w(k), k=1,2,...,n c ierr- an error flag equal to 0 on normal return, c equal to i if the QR algorithm does not c converge within 30 iterations on evaluating the c i-th eigenvalue, equal to -1 if n is not in c range, and equal to -2 if one of the beta's is c negative. c c The array e is needed for working space. c C .. Scalar Arguments .. REAL EPS INTEGER IERR,N C .. C .. Array Arguments .. REAL ALPHA(N),BETA(N),E(N),WEIGHT(N),ZERO(N) C .. C .. Local Scalars .. REAL B,C,F,G,P,R,S INTEGER I,II,J,K,L,M,MML C .. C .. Intrinsic Functions .. INTRINSIC ABS,SIGN,SQRT C .. IF (N.LT.1) THEN IERR = -1 RETURN END IF IERR = 0 ZERO(1) = ALPHA(1) IF (BETA(1).LT.0.) THEN IERR = -2 RETURN END IF WEIGHT(1) = BETA(1) IF (N.EQ.1) RETURN WEIGHT(1) = 1. E(N) = 0. DO 10 K = 2,N ZERO(K) = ALPHA(K) IF (BETA(K).LT.0.) THEN IERR = -2 RETURN END IF E(K-1) = SQRT(BETA(K)) WEIGHT(K) = 0. 10 CONTINUE DO 80 L = 1,N J = 0 c c Look for a small subdiagonal element. c 20 DO 30 M = L,N IF (M.EQ.N) GO TO 40 IF (ABS(E(M)).LE.EPS* (ABS(ZERO(M))+ABS(ZERO(M+ + 1)))) GO TO 40 30 CONTINUE 40 P = ZERO(L) IF (M.EQ.L) GO TO 80 IF (J.EQ.30) GO TO 120 J = J + 1 c c Form shift. c G = (ZERO(L+1)-P)/ (2.*E(L)) R = SQRT(G*G+1.) G = ZERO(M) - P + E(L)/ (G+SIGN(R,G)) S = 1. C = 1. P = 0. MML = M - L c c For i=m-1 step -1 until l do ... c DO 70 II = 1,MML I = M - II F = S*E(I) B = C*E(I) IF (ABS(F).LT.ABS(G)) GO TO 50 C = G/F R = SQRT(C*C+1.) E(I+1) = F*R S = 1./R C = C*S GO TO 60 50 S = F/G R = SQRT(S*S+1.) E(I+1) = G*R C = 1./R S = S*C 60 G = ZERO(I+1) - P R = (ZERO(I)-G)*S + 2.*C*B P = S*R ZERO(I+1) = G + P G = C*R - B c c Form first component of vector. c F = WEIGHT(I+1) WEIGHT(I+1) = S*WEIGHT(I) + C*F WEIGHT(I) = C*WEIGHT(I) - S*F 70 CONTINUE ZERO(L) = ZERO(L) - P E(L) = G E(M) = 0. GO TO 20 80 CONTINUE c c Order eigenvalues and eigenvectors. c DO 100 II = 2,N I = II - 1 K = I P = ZERO(I) DO 90 J = II,N IF (ZERO(J).GE.P) GO TO 90 K = J P = ZERO(J) 90 CONTINUE IF (K.EQ.I) GO TO 100 ZERO(K) = ZERO(I) ZERO(I) = P P = WEIGHT(I) WEIGHT(I) = WEIGHT(K) WEIGHT(K) = P 100 CONTINUE DO 110 K = 1,N WEIGHT(K) = BETA(1)*WEIGHT(K)*WEIGHT(K) 110 CONTINUE RETURN c c Set error - no convergence to an eigenvalue after 30 iterations. c 120 IERR = L RETURN END c c SUBROUTINE GCHRIS(N,IOPT,A,B,X,HP,HN,ALPHA,BETA,IERR) c c This routine implements the generalized Christoffel theorem in two c special cases. Given the first n recursion coefficients for some c measure dlambda(t), it generates the first n recursion coefficients c for the measure c c dlambda(t)/(t-x) if iopt=1 c dlambda(t)/(t**2-x**2) if iopt=2 c c where x is a real number outside the support interval of dlambda. c The case iopt=1 is identical with option 4 of the ORTHPOL routine c chri. c c On entry, c c n - - is the number of recursion coefficients desired; c type integer c iopt - is an integer selecting the desired modification c of the measure dlambda c a,b - are real arrays of dimension n containing the c first n recursion coefficients for the (monic) c orthogonal polynomials relative to the measure c dlambda c x - - is a real parameter defining the linear resp. c quadratic modification of the measure dlambda c hp - - is the value of the integral of dlambda(t)/(x-t); c type real c hn - - is the value of the integral of dlmabda(t)/(-x-t); c type real c c On retrun, c c alpha,beta - are real arrays of dimension n containing c the first n recursion coefficients of the modified c measure c ierr - is an error flag, type integer, where c ierr=0 on normal return c ierr-1 if iopt is not in range c C .. Scalar Arguments .. REAL HN,HP,X INTEGER IERR,IOPT,N C .. C .. Array Arguments .. REAL A(N),ALPHA(N),B(N),BETA(N) C .. C .. Local Scalars .. REAL E,EH,Q,QH INTEGER K C .. IERR = 0 IF (IOPT.EQ.1) THEN ALPHA(1) = X - B(1)/HP BETA(1) = -HP Q = -B(1)/HP DO 10 K = 2,N E = A(K-1) - X - Q BETA(K) = Q*E Q = B(K)/E ALPHA(K) = Q + E + X 10 CONTINUE RETURN ELSE IF (IOPT.EQ.2) THEN ALPHA(1) = X* (HP+HN)/ (HP-HN) BETA(1) = - (HP-HN)/ (2.*X) Q = -B(1)/HP QH = -HP/BETA(1) E = 0. DO 20 K = 2,N EH = Q + E + 2.*X - QH BETA(K) = QH*EH E = A(K-1) - X - Q QH = Q*E/EH ALPHA(K) = QH + EH - X Q = B(K)/E 20 CONTINUE RETURN ELSE IERR = 1 RETURN END IF END c c SUBROUTINE GQRAT(N,MCAP,ALPHA,BETA,XIR,XII,IS,ZG,WG,CONST,IERR,E) c c Given the first n+1 recursion coefficients alpha,beta of the c modified measure dlambda(t)/p(t), this routine generates the n-point c Gauss-type rational quadrature rule for the measure dlambda, c including its error constant c c .. Scalar Arguments .. REAL CONST INTEGER IERR,MCAP,N c .. c .. Array Arguments .. The upper bounds XII, XIR and IS are all c MCAP, but this doesn't work in Fortran 77 c if MCAP = 0. REAL ALPHA(N+1),BETA(N+1),E(N),WG(N),XII(*),XIR(*),ZG(N) INTEGER IS(*) c .. c .. Local Scalars .. REAL EPSM,P INTEGER IERRG,IMU,K,MU c .. c .. External Functions .. REAL R1MACH EXTERNAL R1MACH c .. c .. External Subroutines .. EXTERNAL GAUSS c .. c .. Intrinsic Functions .. INTRINSIC REAL integer i1mach, nout nout = i1mach(2) c .. IERR = 0 EPSM = R1MACH(3) c c compute the n-point Gaussian quadrature rule for the modified measure c CALL GAUSS(N,ALPHA,BETA,EPSM,ZG,WG,IERRG,E) IF (IERRG.NE.0) THEN WRITE (nout,FMT=9000) IERRG IERR = 1 RETURN END IF c c compute the error constant and the weights of the desired rational c Gauss-type quadrature rule c CONST = BETA(1) DO 20 K = 1,N CONST = BETA(K+1)*CONST/REAL(2*K* (2*K-1)) IF (MCAP.GT.0) THEN P = 1. IMU = 0 DO 10 MU = 1,MCAP IF (IMU.EQ.0) THEN IF (XII(MU).EQ.0.) THEN P = ((1.+XIR(MU)*ZG(K))**IS(MU))*P ELSE P = (((1.+XIR(MU)*ZG(K))**2+ + (XII(MU)*ZG(K))**2)**IS(MU))*P IMU = 1 END IF ELSE IMU = 0 END IF 10 CONTINUE WG(K) = P*WG(K) END IF 20 CONTINUE RETURN 9000 FORMAT (1X,'ierrg in gauss=',I5) END c c SUBROUTINE KNUM(N,NU0,NUMAX,Z,EPS,A,B,RHO,NU,IERR,ROLD) c c This routine generates c c rho(k)(z)=integral pi(k)(t)dlambda(t)/(z-t), k=0,1,2,...,n, c c where pi(k)(t) is the (monic) k-th degree orthogonal polynomial c with respect to the measure dlambda(t), and the integral is extended c over the support of dlambda. It is assumed that z is a complex c number outside the smallest interval containing the support of c dlambda. The quantities rho(k)(z) are computed as the first n+1 c members of the minimal solution of the basic three-term recurrence c relation c c y(k+1)(z)=(z-a(k))y(k)(z)-b(k)y(k-1)(z), k=0,1,2,..., c c satisfied by the orthogonal polynomials pi(k)(z). c c Input: n - - the largest integer k for which rho(k) is c desired c nu0 - an estimate of the starting backward recurrence c index; if no better estimate is known, set c nu0 = 3*n/2; for Jacobi, Laguerre and Hermite c weight functions, estimates of nu0 are generated c respectively by the routines nu0jac,nu0lag and c nu0her c numax - an integer larger than n cutting off backward c recursion in case of nonconvergence; if nu0 c exceeds numax, then the routine aborts with the c error flag ierr set equal to nu0 c z - - - the variable in rho(k)(z); type complex c eps - - the relative accuracy to which the rho(k) are c desired c a,b - - arrays of dimension numax to be supplied with the c recurrence coefficients a(k-1), b(k-1), k=1,2,..., c numax. c c Output: rho - - an array of dimension n+1 containing the results c rho(k)=rho(k-1)(z), k=1,2,...,n+1; type complex c nu - - the starting backward recurrence index that yields c convergence c ierr - an error flag equal to zero on normal return, equal c to nu0 if nu0 > numax, and equal to numax in c case of nonconvergence. c c The complex array rold of dimension n+1 is used for working space. c c c The arrays rho,rold are assumed to have dimension n+1. c C .. Scalar Arguments .. COMPLEX Z REAL EPS INTEGER IERR,N,NU,NU0,NUMAX C .. C .. Array Arguments .. COMPLEX RHO(*),ROLD(*) REAL A(NUMAX),B(NUMAX) C .. C .. Local Scalars .. COMPLEX R INTEGER J,J1,K,NP1 C .. C .. Intrinsic Functions .. INTRINSIC ABS,CMPLX C .. IERR = 0 NP1 = N + 1 IF (NU0.GT.NUMAX) THEN IERR = NU0 RETURN END IF IF (NU0.LT.NP1) NU0 = NP1 NU = NU0 - 5 DO 10 K = 1,NP1 RHO(K) = (0.,0.) 10 CONTINUE 20 NU = NU + 5 IF (NU.GT.NUMAX) THEN IERR = NUMAX GO TO 60 END IF DO 30 K = 1,NP1 ROLD(K) = RHO(K) 30 CONTINUE R = (0.,0.) DO 40 J = 1,NU J1 = NU - J + 1 R = CMPLX(B(J1),0.)/ (Z-CMPLX(A(J1),0.)-R) IF (J1.LE.NP1) RHO(J1) = R 40 CONTINUE DO 50 K = 1,NP1 IF (ABS(RHO(K)-ROLD(K)).GT.EPS*ABS(RHO(K))) GO TO 20 50 CONTINUE 60 IF (N.EQ.0) RETURN DO 70 K = 2,NP1 RHO(K) = RHO(K)*RHO(K-1) 70 CONTINUE RETURN END c c SUBROUTINE LANCZ(N,NCAP,X,W,ALPHA,BETA,IERR,P0,P1) c c This routine carries out the same task as the routine sti, but c uses the more stable Lanczos method. The meaning of the input c and output parameters is the same as in the routine sti. (This c routine is adapted from the routine RKPW in W.B. Gragg and c W.J. Harrod,The numerically stable reconstruction of Jacobi c matrices from spectral data'', Numer. Math. 44, 1984, 317-335.) c C .. Scalar Arguments .. INTEGER IERR,N,NCAP C .. C .. Array Arguments .. REAL ALPHA(N),BETA(N),P0(NCAP),P1(NCAP),W(NCAP),X(NCAP) C .. C .. Local Scalars .. REAL GAM,PI,RHO,SIG,T,TK,TMP,TSIG,XLAM INTEGER I,K C .. IF (N.LE.0 .OR. N.GT.NCAP) THEN IERR = 1 RETURN ELSE IERR = 0 END IF DO 10 I = 1,NCAP P0(I) = X(I) P1(I) = 0. 10 CONTINUE P1(1) = W(1) DO 30 I = 1,NCAP - 1 PI = W(I+1) GAM = 1. SIG = 0. T = 0. XLAM = X(I+1) DO 20 K = 1,I + 1 RHO = P1(K) + PI TMP = GAM*RHO TSIG = SIG IF (RHO.LE.0.) THEN GAM = 1. SIG = 0. ELSE GAM = P1(K)/RHO SIG = PI/RHO END IF TK = SIG* (P0(K)-XLAM) - GAM*T P0(K) = P0(K) - (TK-T) T = TK IF (SIG.LE.0.) THEN PI = TSIG*P1(K) ELSE PI = (T**2)/SIG END IF TSIG = SIG P1(K) = TMP 20 CONTINUE 30 CONTINUE DO 40 K = 1,N ALPHA(K) = P0(K) BETA(K) = P1(K) 40 CONTINUE RETURN END c c SUBROUTINE RECUR(N,IPOLY,AL,BE,A,B,IERR) c c This subroutine generates the coefficients a(k),b(k), k=0,1,...,n-1, c in the recurrence relation c c p(k+1)(x)=(x-a(k))*p(k)(x)-b(k)*p(k-1)(x), c k=0,1,...,n-1, c c p(-1)(x)=0, p(0)(x)=1, c c for some classical (monic) orthogonal polynomials, and sets b(0) c equal to the total mass of the weight distribution. The results are c stored in the arrays a,b, which hold, respectively, the coefficients c a(k-1),b(k-1), k=1,2,...,n. c c Input: n - - the number of recursion coefficients desired c ipoly-integer identifying the polynomial as follows: c 1=Legendre polynomial on (-1,1) c 2=Legendre polynomial on (0,1) c 3=Chebyshev polynomial of the first kind c 4=Chebyshev polynomial of the second kind c 5=Jacobi polynomial with parameters al=-.5,be=.5 c 6=Jacobi polynomial with parameters al,be c 7=generalized Laguerre polynomial with c parameter al c 8=Hermite polynomial c al,be-input parameters for Jacobi and generalized c Laguerre polynomials c c Output: a,b - arrays containing, respectively, the recursion c coefficients a(k-1),b(k-1), k=1,2,...,n. c ierr -an error flag, equal to 0 on normal return, c equal to 1 if al or be are out of range c when ipoly=6 or ipoly=7, equal to 2 if b(0) c overflows when ipoly=6 or ipoly=7, equal to 3 c if n is out of range, and equal to 4 if ipoly c is not an admissible integer. In the case ierr=2, c the coefficient b(0) is set equal to the largest c machine-representable number. c c The subroutine calls for the function subroutines r1mach,gamma and c alga. The routines gamma and alga , which are included in this c file, evaluate respectively the gamma function and its logarithm for c positive arguments. They are used only in the cases ipoly=6 and c ipoly=7. c C .. Scalar Arguments .. REAL AL,BE INTEGER IERR,IPOLY,N C .. C .. Array Arguments .. REAL A(N),B(N) C .. C .. Local Scalars .. REAL AL2,ALMACH,ALPBE,BE2,FKM1,T INTEGER K C .. C .. External Functions .. REAL ALGA,GAMMA,R1MACH EXTERNAL ALGA,GAMMA,R1MACH C .. C .. Intrinsic Functions .. INTRINSIC ATAN,EXP,LOG,REAL,SQRT C .. IF (N.LT.1) THEN IERR = 3 RETURN END IF ALMACH = LOG(R1MACH(2)) IERR = 0 DO 10 K = 1,N A(K) = 0. 10 CONTINUE IF (IPOLY.EQ.1) THEN B(1) = 2. IF (N.EQ.1) RETURN DO 20 K = 2,N FKM1 = REAL(K-1) B(K) = 1./ (4.-1./ (FKM1*FKM1)) 20 CONTINUE RETURN ELSE IF (IPOLY.EQ.2) THEN A(1) = .5 B(1) = 1. IF (N.EQ.1) RETURN DO 30 K = 2,N A(K) = .5 FKM1 = REAL(K-1) B(K) = .25/ (4.-1./ (FKM1*FKM1)) 30 CONTINUE RETURN ELSE IF (IPOLY.EQ.3) THEN B(1) = 4.*ATAN(1.) IF (N.EQ.1) RETURN B(2) = .5 IF (N.EQ.2) RETURN DO 40 K = 3,N B(K) = .25 40 CONTINUE RETURN ELSE IF (IPOLY.EQ.4) THEN B(1) = 2.*ATAN(1.) IF (N.EQ.1) RETURN DO 50 K = 2,N B(K) = .25 50 CONTINUE RETURN ELSE IF (IPOLY.EQ.5) THEN B(1) = 4.*ATAN(1.) A(1) = .5 IF (N.EQ.1) RETURN DO 60 K = 2,N B(K) = .25 60 CONTINUE RETURN ELSE IF (IPOLY.EQ.6) THEN IF (AL.LE.-1. .OR. BE.LE.-1.) THEN IERR = 1 RETURN ELSE ALPBE = AL + BE A(1) = (BE-AL)/ (ALPBE+2.) T = (ALPBE+1.)*LOG(2.) + ALGA(AL+1.) + ALGA(BE+1.) - + ALGA(ALPBE+2.) IF (T.GT.ALMACH) THEN IERR = 2 B(1) = R1MACH(2) ELSE B(1) = EXP(T) END IF IF (N.EQ.1) RETURN AL2 = AL*AL BE2 = BE*BE A(2) = (BE2-AL2)/ ((ALPBE+2.)* (ALPBE+4.)) B(2) = 4.* (AL+1.)* (BE+1.)/ ((ALPBE+3.)* (ALPBE+2.)**2) IF (N.EQ.2) RETURN DO 70 K = 3,N FKM1 = REAL(K-1) A(K) = .25* (BE2-AL2)/ (FKM1*FKM1* (1.+.5*ALPBE/FKM1)* + (1.+.5* (ALPBE+2.)/FKM1)) B(K) = .25* (1.+AL/FKM1)* (1.+BE/FKM1)* + (1.+ALPBE/FKM1)/ ((1.+.5* (ALPBE+1.)/FKM1)* + (1.+.5* (ALPBE-1.)/FKM1)* + (1.+.5*ALPBE/FKM1)**2) 70 CONTINUE RETURN END IF ELSE IF (IPOLY.EQ.7) THEN IF (AL.LE.-1.) THEN IERR = 1 RETURN ELSE A(1) = AL + 1. B(1) = GAMMA(AL+1.,IERR) IF (IERR.EQ.2) B(1) = R1MACH(2) IF (N.EQ.1) RETURN DO 80 K = 2,N FKM1 = REAL(K-1) A(K) = 2.*FKM1 + AL + 1. B(K) = FKM1* (FKM1+AL) 80 CONTINUE RETURN END IF ELSE IF (IPOLY.EQ.8) THEN B(1) = SQRT(4.*ATAN(1.)) IF (N.EQ.1) RETURN DO 90 K = 2,N B(K) = .5*REAL(K-1) 90 CONTINUE RETURN ELSE IERR = 4 END IF END REAL FUNCTION ALGA(X) c c This is an auxiliary function subroutine (not optimized in any c sense) evaluating the logarithm of the gamma function for positive c arguments x. It is called by the subroutine gamma. The integer m0 c in the first executable statement is the smallest integer m such c that 1*3*5* ... *(2*m+1)/(2**m) is greater than or equal to the c largest machine-representable number. The routine is based on a c rational approximation valid on [.5,1.5] due to W.J. Cody and c K.E. Hillstrom; see Math. Comp. 21, 1967, 198-203, in particular the c case n=7 in Table II. For the computation of m0 it calls upon the c function subroutines t and r1mach. The former, appended below, c evaluates the inverse function t = t(y) of y = t ln t. c C .. Scalar Arguments .. REAL X C .. C .. Local Scalars .. REAL P,SDEN,SNUM,XE,XI INTEGER K,M,M0,MM1 C .. C .. Local Arrays .. REAL CDEN(8),CNUM(8) C .. C .. External Functions .. REAL R1MACH,T EXTERNAL R1MACH,T C .. C .. Intrinsic Functions .. INTRINSIC AINT,INT,LOG,REAL C .. C .. Data statements .. DATA CNUM/4.120843185847770,85.68982062831317,243.175243524421, + -261.7218583856145,-922.2613728801522,-517.6383498023218, + -77.41064071332953,-2.208843997216182/,CDEN/1., + 45.64677187585908,377.8372484823942,951.323597679706, + 846.0755362020782,262.3083470269460,24.43519662506312, + .4097792921092615/ C .. c c The constants in the statement below are exp(1.) and .5*alog(8.). c M0 = 2.71828*T((LOG(R1MACH(2))-1.03972)/2.71828) XI = AINT(X) IF (X-XI.GT..5) XI = XI + 1. M = INT(XI) - 1 c c Computation of log gamma on the standard interval (1/2,3/2] c XE = X - REAL(M) SNUM = CNUM(1) SDEN = CDEN(1) DO 10 K = 2,8 SNUM = XE*SNUM + CNUM(K) SDEN = XE*SDEN + CDEN(K) 10 CONTINUE ALGA = (XE-1.)*SNUM/SDEN c c Computation of log gamma on (0,1/2] c IF (M.EQ.-1) THEN ALGA = ALGA - LOG(X) RETURN ELSE IF (M.EQ.0) THEN RETURN ELSE c c Computation of log gamma on (3/2,5/2] c P = XE IF (M.EQ.1) THEN ALGA = ALGA + LOG(P) RETURN ELSE c c Computation of log gamma for arguments larger than 5/2 c MM1 = M - 1 c c The else-clause in the next statement is designed to avoid possible c overflow in the computation of p in the if-clause, at the expense c of computing many logarithms. c IF (M.LT.M0) THEN DO 20 K = 1,MM1 P = (XE+REAL(K))*P 20 CONTINUE ALGA = ALGA + LOG(P) RETURN ELSE ALGA = ALGA + LOG(XE) DO 30 K = 1,MM1 ALGA = ALGA + LOG(XE+REAL(K)) 30 CONTINUE RETURN END IF END IF END IF END REAL FUNCTION GAMMA(X,IERR) c c This evaluates the gamma function for real positive x, using the c function subroutines alga and r1mach. In case of overflow, the c routine returns the largest machine-representable number and the c error flag ierr=2. c C .. Scalar Arguments .. REAL X INTEGER IERR C .. C .. Local Scalars .. REAL ALMACH,T C .. C .. External Functions .. REAL ALGA,R1MACH EXTERNAL ALGA,R1MACH C .. C .. Intrinsic Functions .. INTRINSIC EXP,LOG C .. ALMACH = LOG(R1MACH(2)) IERR = 0 T = ALGA(X) IF (T.GE.ALMACH) THEN IERR = 2 GAMMA = R1MACH(2) RETURN ELSE GAMMA = EXP(T) RETURN END IF END REAL FUNCTION T(Y) c c This evaluates the inverse function t = t(y) of y = t ln t for c nonnegative y to an accuracy of about one percent. For the c approximation used, see pp. 51-52 in W. Gautschi,Computational c aspects of three-term recurrence relations'', SIAM Rev. 9, 1967, c 24-82. c C .. Scalar Arguments .. REAL Y C .. C .. Local Scalars .. REAL P,Z C .. C .. Intrinsic Functions .. INTRINSIC LOG C .. IF (Y.LE.10.) THEN P = .000057941*Y - .00176148 P = Y*P + .0208645 P = Y*P - .129013 P = Y*P + .85777 T = Y*P + 1.0125 ELSE Z = LOG(Y) - .775 P = (.775-LOG(Z))/ (1.+Z) P = 1./ (1.+P) T = Y*P/Z END IF RETURN END c c SUBROUTINE STI(N,NCAP,X,W,ALPHA,BETA,IERR,P0,P1,P2) c c This routine applies Stieltjes's procedure'' (cf. Section 2.1 of c W. Gautschi,On generating orthogonal polynomials'', SIAM J. Sci. c Statist. Comput. 3, 1982, 289-317) to generate the recursion c coefficients alpha(k), beta(k) , k=0,1,...,n-1, for the discrete c (monic) orthogonal polynomials associated with the inner product c c (f,g)=sum over k from 1 to ncap of w(k)*f(x(k))*g(x(k)). c c The integer n must be between 1 and ncap, inclusive; otherwise, c there is an error exit with ierr=1. The results are stored in the c arrays alpha, beta; the arrays p0, p1, p2 are working arrays. c c If there is a threat of underflow or overflow in the calculation c of the coefficients alpha(k) and beta(k), the routine exits with c the error flag ierr set equal to -k (in the case of underflow) c or +k (in the case of overflow), where k is the recursion index c for which the problem occurs. The former [latter] can often be avoided c by multiplying all weights w(k) by a sufficiently large [small] c scaling factor prior to entering the routine, and, upon exit, divide c the coefficient beta(0) by the same factor. c c This routine should be used with caution if n is relatively close c to ncap, since there is a distinct possibility of numerical c instability developing. (See W. Gautschi,Is the recurrence relation c for orthogonal polynomials always stable?'', BIT, 1993, to appear.) c In that case, the routine lancz should be used. c c The routine uses the function subroutine r1mach. c C .. Scalar Arguments .. INTEGER IERR,N,NCAP C .. C .. Array Arguments .. REAL ALPHA(N),BETA(N),P0(NCAP),P1(NCAP),P2(NCAP),W(NCAP),X(NCAP) C .. C .. Local Scalars .. REAL HUGE,SUM0,SUM1,SUM2,T,TINY INTEGER K,M,NM1 C .. C .. External Functions .. REAL R1MACH EXTERNAL R1MACH C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. TINY = 10.*R1MACH(1) HUGE = .1*R1MACH(2) IERR = 0 IF (N.LE.0 .OR. N.GT.NCAP) THEN IERR = 1 RETURN END IF NM1 = N - 1 c c Compute the first alpha- and beta-coefficient. c SUM0 = 0. SUM1 = 0. DO 10 M = 1,NCAP SUM0 = SUM0 + W(M) SUM1 = SUM1 + W(M)*X(M) 10 CONTINUE ALPHA(1) = SUM1/SUM0 BETA(1) = SUM0 IF (N.EQ.1) RETURN c c Compute the remaining alpha- and beta-coefficients. c DO 20 M = 1,NCAP P1(M) = 0. P2(M) = 1. 20 CONTINUE DO 40 K = 1,NM1 SUM1 = 0. SUM2 = 0. DO 30 M = 1,NCAP c c The following statement is designed to avoid an overflow condition c in the computation of p2(m) when the weights w(m) go to zero c faster (and underflow) than the p2(m) grow. c IF (W(M).EQ.0.) GO TO 30 P0(M) = P1(M) P1(M) = P2(M) P2(M) = (X(M)-ALPHA(K))*P1(M) - BETA(K)*P0(M) c c Check for impending overflow. c IF (ABS(P2(M)).GT.HUGE .OR. ABS(SUM2).GT.HUGE) THEN IERR = K RETURN END IF T = W(M)*P2(M)*P2(M) SUM1 = SUM1 + T SUM2 = SUM2 + T*X(M) 30 CONTINUE c c Check for impending underflow. c IF (ABS(SUM1).LT.TINY) THEN IERR = -K RETURN END IF ALPHA(K+1) = SUM2/SUM1 BETA(K+1) = SUM1/SUM0 SUM0 = SUM1 40 CONTINUE RETURN END c c INTEGER FUNCTION NU0HER(N,Z,EPS) c c This is an auxiliary function routine providing a starting backward c recurrence index for the Hermite measure that can be used in place c of nu0 in the routines knum and dknum. c C .. Scalar Arguments .. COMPLEX Z REAL EPS INTEGER N C .. C .. Intrinsic Functions .. INTRINSIC ABS,AIMAG,LOG,REAL,SQRT C .. NU0HER = 2.* (SQRT(.5*REAL(N+1))+.25*LOG(1./EPS)/ABS(AIMAG(Z)))**2 RETURN END c c INTEGER FUNCTION NU0JAC(N,Z,EPS) c c This is an auxiliary function routine providing a starting backward c recurrence index for the Jacobi measure that can be used in place c of nu0 in the routines knum and dknum. c C .. Scalar Arguments .. COMPLEX Z REAL EPS INTEGER N C .. C .. Local Scalars .. REAL ANGLE,PI,R,X,X2,Y,Y2 C .. C .. Intrinsic Functions .. INTRINSIC ABS,AIMAG,ATAN,COS,LOG,REAL,SIN,SQRT C .. PI = 4.*ATAN(1.) X = REAL(Z) Y = ABS(AIMAG(Z)) IF (X.LT.1.) THEN IF (X.LT.-1.) ANGLE = .5* (2.*PI+ATAN(Y/ (X-1.))+ + ATAN(Y/ (X+1.))) IF (X.EQ.-1.) ANGLE = .5* (1.5*PI-ATAN(.5*Y)) IF (X.GT.-1.) ANGLE = .5* (PI+ATAN(Y/ (X-1.))+ATAN(Y/ (X+1.))) ELSE IF (X.EQ.1.) ANGLE = .5* (.5*PI+ATAN(.5*Y)) IF (X.GT.1.) ANGLE = .5* (ATAN(Y/ (X-1.))+ATAN(Y/ (X+1.))) END IF X2 = X*X Y2 = Y*Y R = ((X2-Y2-1.)**2+4.*X2*Y2)**.25 R = SQRT((X+R*COS(ANGLE))**2+ (Y+R*SIN(ANGLE))**2) NU0JAC = REAL(N+1) + .5*LOG(1./EPS)/LOG(R) RETURN END c c INTEGER FUNCTION NU0LAG(N,Z,AL,EPS) c c This is an auxiliary function routine providing a starting backward c recurrence index for the Laguerre measure that can be used in place c of nu0 in the routines knum and dknum. c C .. Scalar Arguments .. COMPLEX Z REAL AL,EPS INTEGER N C .. C .. Local Scalars .. REAL PHI,PI,X,Y C .. C .. Intrinsic Functions .. INTRINSIC AIMAG,ATAN,COS,LOG,REAL,SQRT C .. PI = 4.*ATAN(1.) X = REAL(Z) Y = AIMAG(Z) PHI = .5*PI IF (Y.LT.0.) PHI = 1.5*PI IF (X.EQ.0.) GO TO 10 PHI = ATAN(Y/X) IF (Y.GT.0. .AND. X.GT.0.) GO TO 10 PHI = PHI + PI IF (X.LT.0.) GO TO 10 PHI = PHI + PI 10 NU0LAG = (SQRT(REAL(N+1)+.5* (AL+1.))+ + LOG(1./EPS)/ (4.* (X*X+Y*Y)**.25*COS(.5* (PHI-PI))))**2 - + .5* (AL+1.) RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. # End of shell archive exit 0