C ALGORITHM 460, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM C VOL. 16, NO. 10, October, 1973, PP.633--635. #! /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: # Fortran/ # Fortran/Dp/ # Fortran/Dp/Drivers/ # Fortran/Dp/Drivers/Makefile # Fortran/Dp/Drivers/driver.f # Fortran/Dp/Drivers/res # Fortran/Dp/Src/ # Fortran/Dp/Src/src.f # This archive created: Thu Dec 15 13:28:12 2005 export PATH; PATH=/bin:$PATH if test ! -d 'Fortran' then mkdir 'Fortran' fi cd 'Fortran' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << "SHAR_EOF" > 'Makefile' all: Res src.o: src.f $(F77) $(F77OPTS) -c src.f driver.o: driver.f $(F77) $(F77OPTS) -c driver.f DRIVERS= driver RESULTS= Res Objs1= driver.o src.o driver: $(Objs1) $(F77) $(F77OPTS) -o driver $(Objs1) $(SRCLIBS) Res: driver ./driver >Res diffres:Res res echo "Differences in results from driver" $(DIFF) Res res clean: rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS) SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << "SHAR_EOF" > 'driver.f' program main c*********************************************************************** c cc TOMS460_PRB tests ADIP. c implicit none integer n parameter ( n = 10 ) double precision a double precision b double precision c double precision d integer i integer ier integer iopt integer itns double precision dmu double precision omeh(n) double precision omev(n) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOMS460_PRB' write ( *, '(a)' ) ' Test TOMS algorithm 460, optimal' write ( *, '(a)' ) ' ADI parameters.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Given ITNS, compute DMU:' write ( *, '(a)' ) ' ' a = 0.5D+00 b = 2.0D+00 c = 0.25D+00 d = 0.75D+00 iopt = 1 do itns = 1, 5 call adip ( a, b, c, d, iopt, n, itns, dmu, omeh, & omev, ier ) write ( *, '(2x,i6,2x,g14.6,2x,i6)' ) itns, dmu, ier end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Given DMU, compute ITNS:' write ( *, '(a)' ) ' ' a = 0.5D+00 b = 2.0D+00 c = 0.25D+00 d = 0.75D+00 iopt = 2 dmu = 1.0D+00 do i = 1, 8 dmu = dmu / 10.0D+00 call adip ( a, b, c, d, iopt, n, itns, dmu, omeh, & omev, ier ) write ( *, '(2x,g14.6,2x,i6,2x,i6)' ) dmu, itns, ier end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOMS460_PRB' write ( *, '(a)' ) ' Normal end of execution.' stop end SHAR_EOF fi # end of overwriting check if test -f 'res' then echo shar: will not over-write existing file "'res'" else cat << "SHAR_EOF" > 'res' TOMS460_PRB Test TOMS algorithm 460, optimal ADI parameters. Given ITNS, compute DMU: 1 0.775958E-01 0 2 0.150983E-02 0 3 0.293334E-04 0 4 0.569896E-06 0 5 0.110721E-07 0 Given DMU, compute ITNS: 0.100000 1 0 0.100000E-01 2 0 0.100000E-02 3 0 0.100000E-03 3 0 0.100000E-04 4 0 0.100000E-05 4 0 0.100000E-06 5 0 0.100000E-07 6 0 TOMS460_PRB Normal end of execution. SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << "SHAR_EOF" > 'src.f' SUBROUTINE ADIP ( A, B, C, D, IOPT, N, ITNS, DMU, OMEH, & OMEV, IER ) INTEGER N, IER, IOPT, ITNS, J DOUBLE PRECISION A, ALFA, B, BETA, BMD, BPD, C, CMA, CPA, & D, DEL, DKPR, DM, DMU, DRJ, OJ, & OMEH(N), OMEV(N), PISQ, TEMP, TEMPA, TEMPB, TEMPC DATA PISQ / 9.869604401089359D0 / C GIVEN A MATRIX EQUATION GZ=S, WHERE G IS A REAL POSITIVE C DEFINITE MATRIX, S IS A KNOWN, AND Z THE UNKNOWN, VECTOR, C LET H AND V BE SYMMETRIC COMMUTING MATRICES SUCH THAT C G=H+V. BEGINNING WITH AN INITIAL APPROXIMATION Z(0), LET C Z(K+1/2)=(H+OMEH(K)*I)**(-1)*(S-(V-OMEH(K)*I)*Z(K), C Z(K+1) =(V+OMEGV(K)*I)**(-1)*(S-(H-OMEV(K)*I)*Z(K+1/2), C WHERE I IS THE IDENTITY MATRIX. FINALLY, LET ONE OF ITNS C AND DMU BE GIVEN. THEN THIS SUBROUTINE COMPUTES THE PAR- C AMETERS OMEH(K), OMEV(K) THAT MINIMIZE THE VALUE OF C DMU AND ITNS WHICH IS NOT GIVEN WHILE SATISFYING THE C INEQUALITY /Z-Z(ITNS)/.LE.DMU*/Z-Z(0)/, WHERE // DENOTES C THE EUCLIDEAN NORM. C THE SUBROUTINE ARGUMENTS HAVE THE FOLOWING MEANING. C A AND D ARE LOWER AND UPPER BOUNDS, RESPECTIVELY, ON THE C EIGENVALUES OF H. C AND D ARE LOWER AND UPPER BOUNDS, C RESPECTIVELY, ON THE EIGENVALUES OF V. THE VALUES OF A, B, C C AND D MUST SATISFY THE INEQUALITIES 0.LT.A.LE.B AND C 0.LT.C.LE.D. C IOPT DENOTES THE INPUT OPTION. IF IOPT=1, THEN THE VALUE OF C ITNS MUST BE SPECIFIED ON ENTRY AND DMU WILL BE COMPUTED. C IF IOPT = 2 THEN THE VALULE OF DMU MUST BE SPECIFIED ON ENTRY C AND ITNS WILL BE COMPUTED. C IF IOPT=1 THEN THE INEQULAITY 1.LE.ITNS.LE.N MUST BE C SATISFIED, WHILE IF IOPT=2 THEN THE INEQUALITIES N.GE.1 C AND DMU.GT.0 MUST BE SATISFIED. C N IS THE DIMENSION OF THE ARRAYS OMEV AND OMEH. C ITNS IS THE NUMBER OF ITERATIONS TO BE PERFORMED. C DMU IS A BOUND ON THE SPECTRAL NORM OF THE ITERATION C MATRIX TO THE ITNS POWER. IF IOPT=2 A VALUE FOR DMU MUST C BE SPECIFIED ON ENTRY, AND THIS VALUE MAY BE CHANGED BY C ADIP (SEE IER, BELOW). C THE VALUES OF THE REQUIRED PARAMETERS ARE CONTAINED IN THE C LOCATIONS OMEV(K), OMEH(K), K=1,...,ITNS ON EXIT FROM C ADIP. C IER IS A VARIABLE WHOSE VALUE ON EXIT FROM ADIP HAS THE C FOLLOWING MEANING C IER=0 SIGNIFIES COMPUTATION OF THE PARAMETERS HAS BEEN C PERFORMED WITH NO CHANGED OF THE VALUES SPECIFIED ON ENTRY. C IER=1 SIGNIFIES THAT SOME INPUT VALUE VIOLATES THE C CONSTRAINTS GIVEN ABOVE, AND HENCE THE PARAMETERS HAVE NOT C BEEN COMPUTED. C IER=2 (POSSIBLE ONLY IF IOPT=2) SIGNIFIES THAT FOR THE C GIVEN VALUE OF DMU, THE COMPUTED VALUE OF ITNS WOULD BE C GREATER THAN N, SO THAT ITNS HAS BEEN SET EQUAL TO N AND C DMU HAS BEEN RECOMPUTED AS FOR IOPT=1. C TEST INPUT VALUES FOR RANGE. IER = 1 IF ( .NOT. ( A.GT.0.D0 .AND. A.LE.B .AND. C.GT.0.D0 .AND. & C .LE. D ) ) GO TO 90 IF ( .NOT. ( IOPT .EQ. 1 .OR. IOPT .EQ. 2 ) ) GO TO 90 IF ( IOPT .EQ. 2 ) GO TO 10 IF ( .NOT. ( ITNS .GE. 1 .AND. ITNS .LE. N ) ) GO TO 90 GO TO 20 10 IF ( .NOT. (N.GE.1 .AND. ITNS.LE.N ) ) GO TO 90 C STATE 1 - PRELIMINARY COMPUTATIONS COMMON TO BOTH OPTIONS. 20 IER = 0 BPD = B + D BMD = B - D CPA = C + A CMA = C - A DM = 2.D0 * ( ( D - C ) * ( B - A ) ) / ( BPD * CPA ) DKPR = 1.D0 / ( 1.D0 + DM + DSQRT ( DM * ( DM + 2.D0 ) ) ) DEL = 0.D0 IF ( BMD .EQ. 0.D0 .AND. CMA .EQ. 0.D0 ) GO TO 30 TEMP = BPD * DKPR DEL = 2.D0 * ( TEMP - CPA ) / ( CPA * BMD + TEMP * CMA ) 30 ALFA = DKPR * ( CMA + 2.D0 * DEL * A * C ) / CPA BETA = ( 2.D0 + DEL * BMD ) / BPD TEMP = DKPR / 4.D0 C END OF STAGE 1 - COMPUTE ITNS FOR OPTION 2. IF ( IOPT .EQ. 1 ) GO TO 40 ITNS = ( DLOG ( DMU / 4.D0 ) * DLOG ( TEMP ) ) / PISQ + 1.D0 IF ( ITNS .LE. N ) GO TO 40 ITNS = N IER = 2 C STAGE 2 - COMPUTATION OF THE OPTIMAL PARAMETERS. 40 TEMPA = 2 * ITNS TEMPB = TEMP * TEMP DO 50 J = 1, ITNS DRJ = 2 * J - 1 DRJ = DRJ / TEMPA TEMPC = TEMP**DRJ OJ = 2.D0 + ( TEMPC+TEMPB/TEMPC ) / ( 1.D0 + TEMPC * TEMPC ) TEMPC = DEL + OJ OMEV(J) = ( OJ - ALFA ) / ( BETA - TEMPC ) OMEH(J) = ( OJ + ALFA ) / ( BETA + TEMPC ) 50 CONTINUE IF ( IOPT .EQ. 2 .AND. IER .EQ. 0 ) GO TO 90 C END OF STAGE 2 - COMPUTE DMU FOR OPTION 1. TEMPA = ITNS TEMP = PISQ * TEMPA / DLOG ( TEMPB * ( 1.D0 + 8.D0 * TEMPB ) ) C CHOOSE PROPER FORMULA TO AVOID UNDERFLOW OR OVERFLOW. IF ( TEMP .LE. -90.D0 .OR. TEMP .GE. 30.D0 ) GO TO 60 IF ( TEMP .LE. -10.D0 ) GO TO 70 IF ( TEMP .LT. 10.D0 ) GO TO 80 DMU = DEXP ( -6.D0 * TEMP ) GO TO 90 60 DMU = 0.D0 GO TO 90 70 DMU = 4.D0 * DEXP ( 2.0 * TEMP ) GO TO 90 80 TEMP = DEXP ( TEMP ) DMU = ( ( 2.D0 * TEMP ) / ( 1.D0 + 2.0D0 * TEMP**4 ) )**2 90 RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. # End of shell archive exit 0