C ALGORITHM 682, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 16, NO. 2, PP. 158-168. INTEGER FUNCTION I1MACH(I) C***BEGIN PROLOGUE I1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 860501 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE RETURN INTEGER MACHINE DEPENDENT CONSTANTS. C***DESCRIPTION C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C THESE MACHINE CONSTANT ROUTINES MUST BE ACTIVATED FOR C A PARTICULAR ENVIRONMENT. C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C I1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C K = I1MACH(I) C C WHERE I=1,...,16. THE (OUTPUT) VALUE OF K ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C I/O UNIT NUMBERS. C I1MACH( 1) = THE STANDARD INPUT UNIT. C I1MACH( 2) = THE STANDARD OUTPUT UNIT. C I1MACH( 3) = THE STANDARD PUNCH UNIT. C I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT. C C WORDS. C I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT. C I1MACH( 6) = THE NUMBER OF CHARACTERS PER INTEGER STORAGE UNIT. C C INTEGERS. C ASSUME INTEGERS ARE REPRESENTED IN THE S-DIGIT, BASE-A FORM C C SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) C C WHERE 0 .LE. X(I) .LT. A FOR I=0,...,S-1. C I1MACH( 7) = A, THE BASE. C I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS. C I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE. C C FLOATING-POINT NUMBERS. C ASSUME FLOATING-POINT NUMBERS ARE REPRESENTED IN THE T-DIGIT, C BASE-B FORM C SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T, C 0 .LT. X(1), AND EMIN .LE. E .LE. EMAX. C I1MACH(10) = B, THE BASE. C C SINGLE-PRECISION C I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS. C I1MACH(12) = EMIN, THE SMALLEST EXPONENT E. C I1MACH(13) = EMAX, THE LARGEST EXPONENT E. C C DOUBLE-PRECISION C I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS. C I1MACH(15) = EMIN, THE SMALLEST EXPONENT E. C I1MACH(16) = EMAX, THE LARGEST EXPONENT E. C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. ALSO, THE VALUES OF C I1MACH(1) - I1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY C WITH THE LOCAL OPERATING SYSTEM. C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***ROUTINES CALLED (NONE) C***END PROLOGUE I1MACH C INTEGER IMACH(16),OUTPUT EQUIVALENCE (IMACH(4),OUTPUT) C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA IMACH( 1) / 7 / C DATA IMACH( 2) / 2 / C DATA IMACH( 3) / 2 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 33 / C DATA IMACH( 9) / Z1FFFFFFFF / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -256 / C DATA IMACH(13) / 255 / C DATA IMACH(14) / 60 / C DATA IMACH(15) / -256 / C DATA IMACH(16) / 255 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -50 / C DATA IMACH(16) / 76 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -32754 / C DATA IMACH(16) / 32780 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C FOR FTN4 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6LOUTPUT / C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / 00007777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -929 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -929 / C DATA IMACH(16) / 1069 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C FOR FTN5 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / L"OUTPUT" / C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / O"00007777777777777777" / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -929 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -929 / C DATA IMACH(16) / 1069 / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA IMACH( 1) / 100 / C DATA IMACH( 2) / 101 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 101 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 777777777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -8189 / C DATA IMACH(13) / 8190 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -8099 / C DATA IMACH(16) / 8190 / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C DATA IMACH( 1) / 11 / C DATA IMACH( 2) / 12 / C DATA IMACH( 3) / 8 / C DATA IMACH( 4) / 10 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / 3 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 23 / C DATA IMACH( 9) / 8388607 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 38 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 43 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 63 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 4 / C DATA IMACH( 4) / 1 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 39 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 4 / C DATA IMACH( 4) / 1 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 55 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 9000 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 7 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 32 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1015 / C DATA IMACH(16) / 1017 / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 16 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z7FFFFFFF / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 54 / C DATA IMACH(15) / -101 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 62 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 1 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 60 / C DATA IMACH(15) / -1024 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE VAX 11/780 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE ELXSI 6400 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 32 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1022 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE Z80 MICROPROCESSOR C C DATA IMACH( 1) / 1 / C DATA IMACH( 2) / 1 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 1 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE IBM PC - MICROSOFT FORTRAN C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1022 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE IBM PC - PROFESSIONAL FORTRAN C AND LAHEY FORTRAN C C DATA IMACH( 1) / 4 / C DATA IMACH( 2) / 7 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1022 / C DATA IMACH(16) / 1023 / C C***FIRST EXECUTABLE STATEMENT I1MACH IF (I .LT. 1 .OR. I .GT. 16) GO TO 10 C I1MACH=IMACH(I) RETURN C 10 CONTINUE WRITE(OUTPUT,9000) 9000 FORMAT('1ERROR 1 IN I1MACH - I OUT OF BOUNDS ') C C CALL FDUMP C C STOP END DOUBLE PRECISION FUNCTION D1MACH(I) C***BEGIN PROLOGUE D1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 860501 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE RETURN DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS. C***DESCRIPTION C C D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C D = D1MACH(I) C C WHERE I=1,...,5. THE (OUTPUT) VALUE OF D ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C DOUBLE-PRECISION MACHINE CONSTANTS C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. C D1MACH( 5) = LOG10(B) C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***ROUTINES CALLED XERROR C***END PROLOGUE D1MACH C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) C DOUBLE PRECISION DMACH(5) C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C FOR FTN4 C C DATA SMALL(1) / 00564000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C C DATA LARGE(1) / 37757777777777777777B / C DATA LARGE(2) / 37157777777777777777B / C C DATA RIGHT(1) / 15624000000000000000B / C DATA RIGHT(2) / 00000000000000000000B / C C DATA DIVER(1) / 15634000000000000000B / C DATA DIVER(2) / 00000000000000000000B / C C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C FOR FTN5 C C DATA SMALL(1) / O"00564000000000000000" / C DATA SMALL(2) / O"00000000000000000000" / C C DATA LARGE(1) / O"37757777777777777777" / C DATA LARGE(2) / O"37157777777777777777" / C C DATA RIGHT(1) / O"15624000000000000000" / C DATA RIGHT(2) / O"00000000000000000000" / C C DATA DIVER(1) / O"15634000000000000000" / C DATA DIVER(2) / O"00000000000000000000" / C C DATA LOG10(1) / O"17164642023241175717" / C DATA LOG10(2) / O"16367571421742254654" / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777774B / C C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC DMACH(5) C C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ C DATA LOG10/40423K,42023K,50237K,74776K/ C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1), SMALL(2) / '20000000, '00000201 / C DATA LARGE(1), LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1), RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1), DIVER(2) / '20000000, '00000334 / C DATA LOG10(1), LOG10(2) / '23210115, '10237777 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA SMALL(1), SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1), LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1), RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1), DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1), LOG10(2) / O776464202324, O117571775714 / C C MACHINE CONSTANTS FOR THE HP 2100 C THREE WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B / C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B / C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C MACHINE CONSTANTS FOR THE HP 2100 C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177777B / C DATA LARGE(3), LARGE(4) / 177777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 225B / C DATA DIVER(1), DIVER(2) / 40000B, 0 / C DATA DIVER(3), DIVER(4) / 0, 227B / C DATA LOG10(1), LOG10(2) / 46420B, 46502B / C DATA LOG10(3), LOG10(4) / 76747B, 176377B / C C MACHINE CONSTANTS FOR THE HP 9000 C C D1MACH(1) = 2.8480954D-306 C D1MACH(2) = 1.40444776D+306 C D1MACH(3) = 2.22044605D-16 C D1MACH(4) = 4.44089210D-16 C D1MACH(5) = 3.01029996D-1 C C DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B / C DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B / C DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B / C DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B / C DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA SMALL(1), SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1), LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1), DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA SMALL(1), SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1), LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1), RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1), DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1), LOG10(2) / "177464202324, "144117571776 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA SMALL(1), SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1), LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1), RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1), DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1), LOG10(2) / "177464202324, "476747767461 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 8388608, 0 / C DATA LARGE(1), LARGE(2) / 2147483647, -1 / C DATA RIGHT(1), RIGHT(2) / 612368384, 0 / C DATA DIVER(1), DIVER(2) / 620756992, 0 / C DATA LOG10(1), LOG10(2) / 1067065498, -2063872008 / C C DATA SMALL(1), SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1), LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1), RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1), DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1), LOG10(2) / O07746420232, O20476747770 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA SMALL(3), SMALL(4) / 0, 0 / C C DATA LARGE(1), LARGE(2) / 32767, -1 / C DATA LARGE(3), LARGE(4) / -1, -1 / C C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 0 / C C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA DIVER(3), DIVER(4) / 0, 0 / C C DATA LOG10(1), LOG10(2) / 16282, 8346 / C DATA LOG10(3), LOG10(4) / -31493, -12296 / C C DATA SMALL(1), SMALL(2) / O000200, O000000 / C DATA SMALL(3), SMALL(4) / O000000, O000000 / C C DATA LARGE(1), LARGE(2) / O077777, O177777 / C DATA LARGE(3), LARGE(4) / O177777, O177777 / C C DATA RIGHT(1), RIGHT(2) / O022200, O000000 / C DATA RIGHT(3), RIGHT(4) / O000000, O000000 / C C DATA DIVER(1), DIVER(2) / O022400, O000000 / C DATA DIVER(3), DIVER(4) / O000000, O000000 / C C DATA LOG10(1), LOG10(2) / O037632, O020232 / C DATA LOG10(3), LOG10(4) / O102373, O147770 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1), LOG10(2) / O177746420232, O411757177572 / C C MACHINE CONSTANTS FOR VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA LOG10(1), LOG10(2) / 546979738, -805796613 / C C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB / C C MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING) C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1), SMALL(2) / 16, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 15552, 0 / C DATA DIVER(1), DIVER(2) / 15568, 0 / C DATA LOG10(1), LOG10(2) / 1142112243, 2046775455 / C C DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 / C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / C C MACHINE CONSTANTS FOR THE ELXSI 6400 C (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION) C C DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X / C DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X / C DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X / C DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X / C DATA LOG10(1), DIVER(2) / '3FD34413'X,'509F79FF'X / C C MACHINE CONSTANTS FOR THE IBM PC - MICROSOFT FORTRAN C C DATA SMALL(1), SMALL(2) / #00000000, #00100000 / C DATA LARGE(1), LARGE(2) / #FFFFFFFF, #7FEFFFFF / C DATA RIGHT(1), RIGHT(2) / #00000000, #3CA00000 / C DATA DIVER(1), DIVER(2) / #00000000, #3CB00000 / C DATA LOG10(1), LOG10(2) / #509F79FF, #3FD34413 / C C MACHINE CONSTANTS FOR THE IBM PC - PROFESSIONAL FORTRAN C AND LAHEY FORTRAN C C DATA SMALL(1), SMALL(2) / Z'00000000', Z'00100000' / C DATA LARGE(1), LARGE(2) / Z'FFFFFFFF', Z'7FEFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'00000000', Z'3CA00000' / C DATA DIVER(1), DIVER(2) / Z'00000000', Z'3CB00000' / C DATA LOG10(1), LOG10(2) / Z'509F79FF', Z'3FD34413' / C C***FIRST EXECUTABLE STATEMENT D1MACH IF (I .LT. 1 .OR. I .GT. 5) 1 CALL XERROR( 'D1MACH -- I OUT OF BOUNDS',25,1,2) C D1MACH = DMACH(I) RETURN C END REAL FUNCTION R1MACH(I) C***BEGIN PROLOGUE R1MACH C***DATE WRITTEN 790101 (YYMMDD) C***REVISION DATE 860501 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE RETURN SINGLE PRECISION MACHINE DEPENDENT CONSTANTS. C***DESCRIPTION C C R1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C A = R1MACH(I) C C WHERE I=1,...,5. THE (OUTPUT) VALUE OF A ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C SINGLE-PRECISION MACHINE CONSTANTS C R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING. C R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING. C R1MACH(5) = LOG10(B) C***REFERENCES FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR C A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE- C MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978, C PP. 177-188. C***ROUTINES CALLED XERROR C***END PROLOGUE R1MACH C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) C REAL RMACH(5) C EQUIVALENCE (RMACH(1),SMALL(1)) EQUIVALENCE (RMACH(2),LARGE(1)) EQUIVALENCE (RMACH(3),RIGHT(1)) EQUIVALENCE (RMACH(4),DIVER(1)) EQUIVALENCE (RMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA RMACH(1) / Z400800000 / C DATA RMACH(2) / Z5FFFFFFFF / C DATA RMACH(3) / Z4E9800000 / C DATA RMACH(4) / Z4EA800000 / C DATA RMACH(5) / Z500E730E8 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS. C C DATA RMACH(1) / O1771000000000000 / C DATA RMACH(2) / O0777777777777777 / C DATA RMACH(3) / O1311000000000000 / C DATA RMACH(4) / O1301000000000000 / C DATA RMACH(5) / O1157163034761675 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C FOR FTN4 C C DATA RMACH(1) / 00564000000000000000B / C DATA RMACH(2) / 37767777777777777776B / C DATA RMACH(3) / 16414000000000000000B / C DATA RMACH(4) / 16424000000000000000B / C DATA RMACH(5) / 17164642023241175720B / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C FOR FTN5 C C DATA RMACH(1) / O"00564000000000000000" / C DATA RMACH(2) / O"37767777777777777776" / C DATA RMACH(3) / O"16414000000000000000" / C DATA RMACH(4) / O"16424000000000000000" / C DATA RMACH(5) / O"17164642023241175720" / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA RMACH(1) / 200034000000000000000B / C DATA RMACH(2) / 577767777777777777776B / C DATA RMACH(3) / 377224000000000000000B / C DATA RMACH(4) / 377234000000000000000B / C DATA RMACH(5) / 377774642023241175720B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC RMACH(5) C C DATA SMALL/20K,0/,LARGE/77777K,177777K/ C DATA RIGHT/35420K,0/,DIVER/36020K,0/ C DATA LOG10/40423K,42023K/ C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1), SMALL(2) / '20000000, '00000201 / C DATA LARGE(1), LARGE(2) / '37777777, '00000177 / C DATA RIGHT(1), RIGHT(2) / '20000000, '00000352 / C DATA DIVER(1), DIVER(2) / '20000000, '00000353 / C DATA LOG10(1), LOG10(2) / '23210115, '00000377 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA RMACH(1) / O402400000000 / C DATA RMACH(2) / O376777777777 / C DATA RMACH(3) / O714400000000 / C DATA RMACH(4) / O716400000000 / C DATA RMACH(5) / O776464202324 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE91), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C MACHINE CONSTANTS FOR THE HP 9000 C C R1MACH(1) = 1.17549435E-38 C R1MACH(2) = 1.70141163E+38 C R1MACH(3) = 5.960464478E-8 C R1MACH(4) = 1.119209290E-7 C R1MACH(5) = 3.01030010E-1 C C DATA SMALL(1) / 00040000000B / C DATA LARGE(1) / 17677777777B / C DATA RIGHT(1) / 06340000000B / C DATA DIVER(1) / 06400000000B / C DATA LOG10(1) / 07646420233B / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86 AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA RMACH(1) / Z00100000 / C DATA RMACH(2) / Z7FFFFFFF / C DATA RMACH(3) / Z3B100000 / C DATA RMACH(4) / Z3C100000 / C DATA RMACH(5) / Z41134413 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR). C C DATA RMACH(1) / "000400000000 / C DATA RMACH(2) / "377777777777 / C DATA RMACH(3) / "146400000000 / C DATA RMACH(4) / "147400000000 / C DATA RMACH(5) / "177464202324 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1) / 8388608 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 880803840 / C DATA DIVER(1) / 889192448 / C DATA LOG10(1) / 1067065499 / C C DATA RMACH(1) / O00040000000 / C DATA RMACH(2) / O17777777777 / C DATA RMACH(3) / O06440000000 / C DATA RMACH(4) / O06500000000 / C DATA RMACH(5) / O07746420233 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA LARGE(1), LARGE(2) / 32767, -1 / C DATA RIGHT(1), RIGHT(2) / 13440, 0 / C DATA DIVER(1), DIVER(2) / 13568, 0 / C DATA LOG10(1), LOG10(2) / 16282, 8347 / C C DATA SMALL(1), SMALL(2) / O000200, O000000 / C DATA LARGE(1), LARGE(2) / O077777, O177777 / C DATA RIGHT(1), RIGHT(2) / O032200, O000000 / C DATA DIVER(1), DIVER(2) / O032400, O000000 / C DATA LOG10(1), LOG10(2) / O037632, O020233 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C DATA RMACH(1) / O000400000000 / C DATA RMACH(2) / O377777777777 / C DATA RMACH(3) / O146400000000 / C DATA RMACH(4) / O147400000000 / C DATA RMACH(5) / O177464202324 / C C MACHINE CONSTANTS FOR THE VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1) / 128 / C DATA LARGE(1) / -32769 / C DATA RIGHT(1) / 13440 / C DATA DIVER(1) / 13568 / C DATA LOG10(1) / 547045274 / C C DATA SMALL(1) / Z00000080 / C DATA LARGE(1) / ZFFFF7FFF / C DATA RIGHT(1) / Z00003480 / C DATA DIVER(1) / Z00003500 / C DATA LOG10(1) / Z209B3F9A / C C MACHINE CONSTANTS FOR THE ELXSI 6400 C ASSUMING REAL*4 IS THE DEFAULT REAL C C DATA SMALL(1) / '00800000'X / C DATA LARGE(1) / '7F7FFFFF'X / C DATA RIGHT(1) / '33800000'X / C DATA DIVER(1) / '34000000'X / C DATA LOG10(1) / '3E9A209B'X / C C MACHINE CONSTANTS FOR THE Z80 MICROPROCESSOR C C DATA SMALL(1), SMALL(2) / 0, 256 / C DATA LARGE(1), LARGE(2) / -1, -129 / C DATA RIGHT(1), RIGHT(2) / 0, 26880 / C DATA DIVER(1), DIVER(2) / 0, 27136 / C DATA LOG10(1), LOG10(2) / 8347, 32538 / C C MACHINE CONSTANTS FOR THE IBM PC - MICROSOFT FORTRAN C C DATA SMALL(1) / #00800000 / C DATA LARGE(1) / #7F7FFFFF / C DATA RIGHT(1) / #33800000 / C DATA DIVER(1) / #34000000 / C DATA LOG10(1) / #3E9A209A / C C MACHINE CONSTANTS FOR THE IBM PC - PROFESSIONAL FORTRAN C AND LAHEY FORTRAN C C DATA SMALL(1)/ Z'00800000' / C DATA LARGE(1)/ Z'7F7FFFFF' / C DATA RIGHT(1)/ Z'33800000' / C DATA DIVER(1)/ Z'34000000' / C DATA LOG10(1)/ Z'3E9A209A' / C C***FIRST EXECUTABLE STATEMENT R1MACH IF (I .LT. 1 .OR. I .GT. 5) 1 CALL XERROR ( 'R1MACH -- I OUT OF BOUNDS',25,1,2) C R1MACH = RMACH(I) RETURN C END SUBROUTINE FDUMP C***BEGIN PROLOGUE FDUMP C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 850801 (YYMMDD) C***CATEGORY NO. R3 C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE SYMBOLIC DUMP (SHOULD BE LOCALLY WRITTEN). C***DESCRIPTION C C ***NOTE*** MACHINE DEPENDENT ROUTINE C FDUMP IS INTENDED TO BE REPLACED BY A LOCALLY WRITTEN C VERSION WHICH PRODUCES A SYMBOLIC DUMP. FAILING THIS, C IT SHOULD BE REPLACED BY A VERSION WHICH PRINTS THE C SUBPROGRAM NESTING LIST. NOTE THAT THIS DUMP MUST BE C PRINTED ON EACH OF UP TO FIVE FILES, AS INDICATED BY THE C XGETUA ROUTINE. SEE XSETUA AND XGETUA FOR DETAILS. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C***REFERENCES (NONE) C***ROUTINES CALLED (NONE) C***END PROLOGUE FDUMP C***FIRST EXECUTABLE STATEMENT FDUMP RETURN END SUBROUTINE XERROR(MESS,NMESS,NERR,LEVEL) C***BEGIN PROLOGUE XERROR C***DATE WRITTEN 880401 (YYMMDD) C***REVISION DATE 880401 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR C***PURPOSE TO PRINT AN ERROR MESSAGE C***DESCRIPTION C C ABSTRACT C XERROR IS A DUMMY SLATEC LIBRARY ROUTINE TO PRINT AN ERROR C MESSAGE. THE CALL LIST IS CONSISTENT WITH THE CORRESPONDING C SLATEC LIBRARY SUBROUTINE. C C DESCRIPTION OF PARAMETERS C --INPUT-- C MESSG - THE HOLLERITH MESSAGE TO BE PROCESSED. C NMESSG- THE ACTUAL NUMBER OF CHARACTERS IN MESSG. C NERR - THE ERROR NUMBER ASSOCIATED WITH THIS MESSAGE. C (A DUMMY VARIABLE NOT USED, BUT NEEDED TO BE COMPATIBLE C WITH THE CORRESPONDING SLATEC ROUTINE) C LEVEL - ERROR CATEGORY. C (A DUMMY VARIABLE NOT USED, BUT NEEDED TO BE COMPATIBLE C WITH THE CORRESPONDING SLATEC ROUTINE) C C***REFERENCES JONES R.E., KAHANER D.K., 'XERROR, THE SLATEC ERROR- C HANDLING PACKAGE', SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED NONE C***END PROLOGUE XERROR CHARACTER*(*) MESS INTEGER NMESS,NERR,LEVEL,NN,NR,K,I,KMIN,MIN IF (NMESS.LE.0) THEN PRINT *,' IN XERROR, NMESS IS OUT OF RANGE' ELSE NN=NMESS/70 NR=NMESS-70*NN IF(NR.NE.0) NN=NN+1 K=1 PRINT 900 900 FORMAT(/) DO 10 I=1,NN KMIN=MIN(K+69,NMESS) PRINT *, MESS(K:KMIN) K=K+70 10 CONTINUE RETURN END IF C----------------------------------------------------------------------- C END XERROR C----------------------------------------------------------------------- END SUBROUTINE XERRWV(MESSG,NMESSG,NERR,LEVEL,NI,I1,I2,NR,R1,R2) C***BEGIN PROLOGUE XERRWV C***DATE WRITTEN 880401 (YYMMDD) C***REVISION DATE 880401 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR C***PURPOSE PROCESS AN ERROR MESSAGE ALLOWING 2 INTEGER AND 2 REAL C VALUES TO BE INCLUDED IN THE MESSAGE. C***DESCRIPTION C C ABSTRACT C XERRWV IS A DUMMY SLATEC LIBRARY ROUTINE TO PRINT AN ERROR C MESSAGE AND PRINT UP TO 2 INTEGER VARIABLES AND 2 REAL C VARIABLES. THE CALL LIST IS CONSISTENT WITH THE CORRESPONDING C SLATEC LIBRARY SUBROUTINE. C C DESCRIPTION OF PARAMETERS C --INPUT-- C MESSG - THE HOLLERITH MESSAGE TO BE PROCESSED. C NMESSG- THE ACTUAL NUMBER OF CHARACTERS IN MESSG. C NERR - THE ERROR NUMBER ASSOCIATED WITH THIS MESSAGE. C (A DUMMY VARIABLE NOT USED, BUT NEEDED TO BE COMPATIBLE C WITH THE CORRESPONDING SLATEC ROUTINE) C LEVEL - ERROR CATEGORY. C (A DUMMY VARIABLE NOT USED, BUT NEEDED TO BE COMPATIBLE C WITH THE CORRESPONDING SLATEC ROUTINE) C NI - NUMBER OF INTEGER VALUES TO BE PRINTED. (0 TO 2) C I1 - FIRST INTEGER VALUE. C I2 - SECOND INTEGER VALUE. C NR - NUMBER OF REAL VALUES TO BE PRINTED. (0 TO 2) C R1 - FIRST REAL VALUE. C R2 - SECOND REAL VALUE. C C***REFERENCES JONES R.E., KAHANER D.K., 'XERROR, THE SLATEC ERROR- C HANDLING PACKAGE', SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED XERROR C***END PROLOGUE XERRWV CHARACTER*(*) MESSG INTEGER I1,I2,NI,NR,NMESSG,NERR,LEVEL REAL R1,R2 CALL XERROR(MESSG,NMESSG,NERR,LEVEL) IF (NI.LT.0 .OR. NI.GT.2) THEN PRINT *,' IN XERRWV, NI IS OUT OF RANGE' ELSE IF (NI.NE.0) THEN IF (NI.EQ.1) THEN PRINT *,' I1 = ',I1 ELSE PRINT *,' I1 , I2 = ',I1,I2 END IF END IF END IF IF (NR.LT.0 .OR. NR.GT.2) THEN PRINT *,' IN XERRWV, NR IS OUT OF RANGE' ELSE IF (NR.NE.0) THEN IF (NR.EQ.1) THEN PRINT *,' R1 = ',R1 ELSE PRINT *,' R1 , R2 = ',R1,R2 END IF END IF END IF RETURN C----------------------------------------------------------------------- C END XERRWV C----------------------------------------------------------------------- END PROGRAM DRIVER C C C*********************************************************************** C * C * C DRIVER PROGRAM FOR TESTING SOFTWARE (TAPAR AND TSUM) IMPLEMENTING * C * C TALBOT'S METHOD FOR THE NUMERICAL INVERSION OF LAPLACE TRANSFORMS * C * C * C * C * C>>>>>>>>>>> VERSION 3.0 FEBRUARY 1, 1989 <<<<<<<<<<<* C * C * C * C * C AUTHORS: ALMERICO MURLI, MARIAROSARIA RIZZARDI * C * C * C DEPARTMENT OF MATHEMATIC AND APPLICATIONS * C UNIVERSITY OF NAPLES - ITALY * C * C * C*********************************************************************** C C C C REFERENCES C ========== C C MURLI A., M. RIZZARDI - "ALGORITHM XXX: TALBOT'S METHOD FOR THE C LAPLACE INVERSION PROBLEM". C ACM TRANS. MATH. SOFTWARE, VOL. ##, C NO. #, MONTH YEAR, PP. ##-##. C C C C DESCRIPTION C =========== C C A SET OF 8 LAPLACE TRANSFORMS HAS BEEN CHOSEN. C C THE SELECTED LAPLACE TRANSFORMS ARE DIVIDED INTO C C LAPFUN(S) WITH REAL SINGULARITIES: C C 1/S**2 C LOG(S)/S C EXP(-4*SQRT(S)) C EXP(-5*S)/S C EXP(-1/S)/SQRT(S) C C LAPFUN(S) WITH COMPLEX SINGULARITIES: C C ARCTAN(1/S) C LOG((S**2+1)/(S**2+4)) C S**2/(S**3+8) C C MOREOVER, THESE TEST FUNCTIONS CONTAIN AN EXAMPLE OF INVERSE LAPLACE C TRANSFORM TENDING TO INFINITY AS T INCREASES (OVERFLOW PROBLEM): C C INVLAP(T)=(EXP(-2*T)+2*EXP(T)*COS(T*SQRT(3)))/3 C C WHOSE TRANSFORM IS C C LAPFUN(S)=S**2/(S**3+8) C C C AND TWO PARTICULAR FUNCTIONS C C LAPFUN(S)=EXP(-5*S)/S C LAPFUN(S)=EXP(-1/S)/SQRT(S) C C WHICH ARE DISCUSSED IN THE PAPER. C C C FOR EACH OF THE EIGHT FUNCTIONS THIS PROGRAM CALCULATES INVLAP(T) C FOR SEVEN VALUES OF T, AND FOR EACH VALUE OF T THIS PROGRAM CALLS C THE ROUTINES IMPLEMENTING TALBOT'S METHOD FOR SEVERAL REQUIRED C ACCURACIES (MORE PRECISELY (C+2)-TIMES WHERE C IS THE DECIMAL C MACHINE ACCURACY). C C IN SUBPROGRAMS WE PROVIDE A GENERIC LAPLACE TRANSFORM COMPLEX C FUNCTION FLAP(S) AND ITS CORRESPONDING INVERSE FUNCTION FINV(T) C (THE SECOND IN DOUBLE PRECISION). BY MEANS OF A COMMON VARIABLE C (NAMED NF) WE ARE ABLE TO SELECT ONE OF THE 8 TEST FUNCTIONS. C C MOREOVER SINCE TALBOT'S METHOD NEEDS TO KNOW THE LOCATION AND NATURE C OF SINGULARITIES OF FLAP(S), WE PROVIDE A SUBROUTINE SINGS THAT C RETURNS THIS INFORMATION (USING THE SAME COMMON AREA TO ESTABLISH C WHICH FUNCTION ACTUALLY IS TESTING). C C IN THIS PROGRAM THE FOLLOWING STEPS ARE REPEATED FOR EACH TEST C FUNCTION: C C 1) INPUT LOCATION AND NATURE OF THE SINGULARITIES OF THE LAPLACE C TRANSFORM. C C 2) SUBSEQUENT STEPS ARE PERFORMED FOR EACH VALUE OF T: C C 2.A) THE INVERSE LAPLACE TRANSFORM IS COMPUTED IN T USING DOUBLE C PRECISION ARITHMETIC. THIS "REFERENCE VALUE" WILL BE C COMPARED WITH THE NUMERICAL APPROXIMATION TO ESTABLISH ITS C LOCAL ACCURACY. C C 2.B) A LOOP IS MADE VARYING THE REQUESTED ACCURACY FROM 1 CORRECT C DECIMAL DIGIT TILL 2 DIGITS MORE THAN THE "DECIMAL MACHINE C PRECISION" (DECDIG SHOWS THE REQUESTED ACCURACY IN OUTPUTS C OF THIS PROGRAM). C FOR EACH T AND FOR EACH REQUIRED ACCURACY, THE SUBROUTINE C TAPAR IS CALLED TO PROVIDE THE PARAMETERS OF TALBOT'S METHOD C SUBSEQUENTLY THE SUBROUTINE TSUM IS CALLED TO PROVIDE THE C NUMERICAL APPROXIMATION. C OUTPUTS SHOW FOR EACH VALUE OF T, THE VALUES OF THE GEOME- C TRICAL PARAMETERS LAMBDA, SIGMA, NU AND THE REFERENCE VALUE C OF THE INVERSE LAPLACE TRANSFORM. THEN VARYING THE REQUIRED C ACCURACY, THE OUTPUT SHOWS: C - THE REQUESTED ACCURACY IN THE RESULT C - THE NUMBER OF POINTS USED BY THE QUADRATURE RULE C - THE NUMERICAL APPROXIMATION PROVIDED BY THE METHOD C - TWO LOCAL ERRORS (*) C C C C======================================================================= C = C = C (*) ABOUT THE ERRORS, TWO MEASURES HAVE BEEN CHOSEN: = C = C I) WHAT WE CALLED "PSEUDO ERROR" IS THE RATIO = C PS.ERR = ABSOLUTE ERROR / EXP(SIGMA*T) = C WHERE SIGMA IS ONE OF THE GEOMETRICAL PARAMETERS OF THE = C METHOD. GENERALLY SIGMA IS ZERO, BUT WHEN IT IS POSITIVE = C THE METHOD WILL INVERT A TRANSFORM SUCH AS = C LAPFUN(S+SIGMA) INSTEAD OF LAPFUN(S) = C SO THAT THE FINAL ERROR HAS A FACTOR LIKE EXP(SIGMA*T). = C = C II) WHAT WE SIMPLY CALLED "ERROR" IS THE RELATIVE ERROR WHEN = C THE "REFERENCE VALUE" OF THE INVERSE IS GREATER THAN ONE, = C OTHERWISE IT IS THE ABSOLUTE ERROR. = C IN THE OUTPUT A CHARACTER "A" OR "R" SPECIFIES ABSOLUTE OR = C RELATIVE ERROR HAS BEEN COMPUTED RESPECTIVELY. = C = C = C======================================================================= C C C C REMARKS ON INPUT/OUTPUT STATEMENTS C ================================== C C <> NO INPUT IS REQUIRED C C <> OUTPUT IS DIRECTED TO A GENERIC DEVICE, INDIVIDUATED BY THE C INTEGER VARIABLE "OUTPUT" IN ALL THE "WRITE" STATEMENTS. C "OUTPUT" MUST BE PROPERLY DEFINED BEFORE USING THIS PROGRAM C CHANGING THE STATEMENT: C C DATA OUTPUT / 6 / C C THAT ACTUALLY DEFINES THE STANDARD OUTPUT DEVICE. C C C SUBROUTINES OR FUNCTIONS NEEDED C =============================== C C TAPAR, TSUM, C INVPC : SUBROUTINES IMPLEMENTING TALBOT'S METHOD. C C FLAP : COMPLEX FUNCTION THAT PROVIDES THE LAPLACE TRANSFORM C TO BE INVERTED. C C SINGS, FINV : AUXILIARY ROUTINES; THE FORMER PROVIDES THE SINGULA_ C RITIES AND THEIR MULTIPLICITIES OF THE LAPLACE C TRANSFORM SELECTED FOR INVERSION, THE LATTER PROVIDES C THE INVERSE LAPLACE TRANSFORM FUNCTION. C C R1MACH, D1MACH : FUNCTIONS RESPECTIVELY FOR SINGLE AND DOUBLE PRECI_ C SION MACHINE DEPENDENT CONSTANTS. C C ALOG, ALOG10, DABS, DBLE, DEXP, DLOG, IFIX, SQRT : FORTRAN INTRINSIC C FUNCTIONS. C C C C*********************************************************************** C C COMMON /INVLAP/ NF EXTERNAL FLAP COMPLEX FLAP DOUBLE PRECISION FINV,FV,TT,ERROR,PSERR,DINVF,DEMAX,D1MACH REAL INVF INTEGER OUTPUT CHARACTER*1 ALFA C C DIMENSION SINGRE(24),SINGIM(24),MULTSI(24),VT(8) DATA NDIM / 24 / DATA OUTPUT / 6 / C C C C T-VALUES SELECTED FOR TESTING C C DATA VT/.01,.1,1.,10.,100.,700.,1000.,10000./ C C C NT IS THE NUMBER OF T-VALUES USED HERE FOR THE FIRST SIX FUNCTIONS. C THE MAXIMUM ALLOWED IS THE DIMENSION OF THE ARRAY VT(I) C C NT=7 C C C IC IS THE MACHINE DECIMAL PRECISION USED BY TALBOT'S METHOD. HERE C IT IS COMPUTED, TO ASSURE PORTABILITY, USING R1MACH AS FOLLOWS C C IC=-IFIX(ALOG10(R1MACH(3))*3./4.) C C C EMAX IS THE LOGARITHM OF THE LARGEST MAGNITUDE FLOATING-POINT C NUMBER (FOR SINGLE PRECISION) C C EMAX = ALOG(R1MACH(2)) C C C DEMAX IS THE LOGARITHM OF THE LARGEST MAGNITUDE FLOATING-POINT C NUMBER (FOR DOUBLE PRECISION) C C DEMAX= DLOG(D1MACH(2)) C C WRITE(OUTPUT,100) WRITE(OUTPUT,101) WRITE(OUTPUT,102) C C DO 1 NF=1,6 IF(NF .EQ. 1) WRITE(OUTPUT,201) IF(NF .EQ. 2) WRITE(OUTPUT,202) IF(NF .EQ. 3) WRITE(OUTPUT,203) IF(NF .EQ. 4) WRITE(OUTPUT,204) IF(NF .EQ. 5) WRITE(OUTPUT,205) IF(NF .EQ. 6) WRITE(OUTPUT,206) C C C THIS IS THE START OF THE BLOCK FOR TESTING A PARTICULAR FUNCTION C C CALL SINGS(NDIM,SINGRE,SINGIM,MULTSI,NOSING) WRITE(OUTPUT,200) (SINGRE(K),SINGIM(K),MULTSI(K),K=1,NOSING) WRITE(OUTPUT,250) IC WRITE(OUTPUT,251) WRITE(OUTPUT,252) WRITE(OUTPUT,253) WRITE(OUTPUT,254) DO 2 IT=1,NT TVALUE=VT(IT) TT=DBLE(TVALUE) IF(TVALUE .GT. DEMAX .AND. NF .EQ. 6) THEN FV=0.0D0 ELSE FV=FINV(TVALUE) ENDIF IF(DABS(FV) .GT. 1.0D0) THEN ALFA='R' ELSE ALFA='A' ENDIF DO 3 ID=1,IC+2 CALL TAPAR(TVALUE,ID,NOSING,SINGRE,SINGIM,MULTSI,IC, + CONLAM,CONSIG,CONNU,NOPTS) IF(ID .EQ. 1) THEN IF(NF .EQ. 6 .AND. TVALUE .GT. DEMAX) THEN WRITE(OUTPUT,301) TVALUE,CONLAM,CONSIG,CONNU ELSE WRITE(OUTPUT,300) TVALUE,CONLAM,CONSIG,CONNU,FV ENDIF ENDIF CALL TSUM(FLAP,CONLAM,CONSIG,CONNU,NOPTS,TVALUE,INVF,IER) IF(IER .EQ. 0) THEN PSERR=(FV-INVF)/DEXP(TT*CONSIG) IF(ALFA .EQ. 'R') THEN ERROR=(FV-INVF)/FV ELSE ERROR= FV-INVF ENDIF WRITE(OUTPUT,350) ID,NOPTS,INVF,PSERR,ERROR,ALFA ELSE WRITE(OUTPUT,351) ID,NOPTS,INVF,IER IF(NF .EQ. 6 .AND.(TVALUE.GT.EMAX.AND.TVALUE.LT.DEMAX)) THEN DINVF=CONLAM/NOPTS*DEXP(TT*CONSIG)*INVF PSERR=(FV-DINVF)/DEXP(TT*CONSIG) ERROR=(FV-DINVF)/FV WRITE(OUTPUT,352) DINVF,PSERR,ERROR,ALFA ENDIF ENDIF 3 CONTINUE 2 CONTINUE C C C THIS IS THE END OF THE BLOCK FOR TESTING A PARTICULAR FUNCTION C C 1 CONTINUE C C C THE FOLLOWING BLOCKS REFER TO PARTICULAR SITUATIONS C DISCUSSED IN THE PAPER C C C C C THE FOLLOWING BLOCK REFERS TO AN EXAMPLE OF THE C NECESSITY OF CONDITION (2.2A) C C NF=7 WRITE(OUTPUT,207) CALL SINGS(NDIM,SINGRE,SINGIM,MULTSI,NOSING) WRITE(OUTPUT,200) (SINGRE(K),SINGIM(K),MULTSI(K),K=1,NOSING) WRITE(OUTPUT,250) IC WRITE(OUTPUT,251) WRITE(OUTPUT,252) WRITE(OUTPUT,253) WRITE(OUTPUT,254) DO 4 IT=4,NT TVALUE=VT(IT) TT=DBLE(TVALUE) FV=FINV(TVALUE) IF(DABS(FV) .GT. 1.0D0) THEN ALFA='R' ELSE ALFA='A' ENDIF DO 5 ID=1,IC+2 CALL TAPAR(TVALUE,ID,NOSING,SINGRE,SINGIM,MULTSI,IC, + CONLAM,CONSIG,CONNU,NOPTS) IF(ID .EQ. 1) WRITE(OUTPUT,300) TVALUE,CONLAM,CONSIG,CONNU,FV CALL TSUM(FLAP,CONLAM,CONSIG,CONNU,NOPTS,TVALUE,INVF,IER) IF(IER .EQ. 0) THEN PSERR=(FV-INVF)/DEXP(TT*CONSIG) IF(ALFA .EQ. 'R') THEN ERROR=(FV-INVF)/FV ELSE ERROR=FV-INVF ENDIF WRITE(OUTPUT,350) ID,NOPTS,INVF,PSERR,ERROR,ALFA ELSE WRITE(OUTPUT,351) ID,NOPTS,INVF,IER ENDIF 5 CONTINUE 4 CONTINUE C C C THE FOLLOWING THREE BLOCKS REFER TO AN EXAMPLE OF THE PROBLEM C OF ESSENTIAL SINGULARITIES C C C C C 1) NORMAL TALBOT'S METHOD C C NF=8 NT=6 WRITE(OUTPUT,208) CALL SINGS(NDIM,SINGRE,SINGIM,MULTSI,NOSING) WRITE(OUTPUT,200) (SINGRE(K),SINGIM(K),MULTSI(K),K=1,NOSING) WRITE(OUTPUT,250) IC WRITE(OUTPUT,251) WRITE(OUTPUT,252) WRITE(OUTPUT,253) WRITE(OUTPUT,254) DO 6 IT=1,NT TVALUE=VT(IT) TT=DBLE(TVALUE) FV=FINV(TVALUE) IF(DABS(FV) .GT. 1.0D0) THEN ALFA='R' ELSE ALFA='A' ENDIF DO 7 ID=1,IC+2 CALL TAPAR(TVALUE,ID,NOSING,SINGRE,SINGIM,MULTSI,IC, + CONLAM,CONSIG,CONNU,NOPTS) IF(ID .EQ. 1) WRITE(OUTPUT,300) TVALUE,CONLAM,CONSIG,CONNU,FV CALL TSUM(FLAP,CONLAM,CONSIG,CONNU,NOPTS,TVALUE,INVF,IER) IF(IER .EQ. 0) THEN PSERR=(FV-INVF)/DEXP(TT*CONSIG) IF(ALFA .EQ. 'R') THEN ERROR=(FV-INVF)/FV ELSE ERROR=FV-INVF ENDIF WRITE(OUTPUT,350) ID,NOPTS,INVF,PSERR,ERROR,ALFA ELSE WRITE(OUTPUT,351) ID,NOPTS,INVF,IER ENDIF 7 CONTINUE 6 CONTINUE C C C 2) TALBOT'S MODIFICATION C C NT=6 WRITE(OUTPUT,208) WRITE(OUTPUT,400) DO 8 IT=1,NT TVALUE=VT(IT) TT=DBLE(TVALUE) FV=FINV(TVALUE) IF(DABS(FV) .GT. 1.0D0) THEN ALFA='R' ELSE ALFA='A' ENDIF DO 9 ID=1,IC+2 CALL TAPAR(TVALUE,ID,NOSING,SINGRE,SINGIM,MULTSI,IC, + CONLAM,CONSIG,CONNU,NOPTS) OMEGA=CONLAM*TVALUE CONLAM=(OMEGA-1.)/TVALUE+1.0/30.0 IF(ID .EQ. 1) WRITE(OUTPUT,300) TVALUE,CONLAM,CONSIG,CONNU,FV CALL TSUM(FLAP,CONLAM,CONSIG,CONNU,NOPTS,TVALUE,INVF,IER) IF(IER .EQ. 0) THEN PSERR=(FV-INVF)/DEXP(TT*CONSIG) IF(ALFA .EQ. 'R') THEN ERROR=(FV-INVF)/FV ELSE ERROR=FV-INVF ENDIF WRITE(OUTPUT,350) ID,NOPTS,INVF,PSERR,ERROR,ALFA ELSE WRITE(OUTPUT,351) ID,NOPTS,INVF,IER ENDIF IF(TVALUE .LE. 10.0) GOTO 9 NOPTS=NOPTS*ALOG10(TVALUE) CALL TSUM(FLAP,CONLAM,CONSIG,CONNU,NOPTS,TVALUE,INVF,IER) IF(IER .EQ. 0) THEN PSERR=(FV-INVF)/DEXP(TT*CONSIG) IF(ALFA .EQ. 'R') THEN ERROR=(FV-INVF)/FV ELSE ERROR=FV-INVF ENDIF WRITE(OUTPUT,350) ID,NOPTS,INVF,PSERR,ERROR,ALFA ELSE WRITE(OUTPUT,351) ID,NOPTS,INVF,IER ENDIF 9 CONTINUE 8 CONTINUE C C C 3) OUR MODIFICATION (SEE THE PAPER) C C NT=6 WRITE(OUTPUT,208) WRITE(OUTPUT,401) DO 10 IT=1,NT TVALUE=VT(IT) TT=DBLE(TVALUE) FV=FINV(TVALUE) IF(DABS(FV) .GT. 1.0D0) THEN ALFA='R' ELSE ALFA='A' ENDIF DO 11 ID=1,IC+2 CALL TAPAR(TVALUE,ID,NOSING,SINGRE,SINGIM,MULTSI,IC, + CONLAM,CONSIG,CONNU,NOPTS) OMEGA=CONLAM*TVALUE CONLAM=(OMEGA-1.0)/TVALUE+1.0/SQRT(OMEGA*TVALUE) IF(ID .EQ. 1) WRITE(OUTPUT,300) TVALUE,CONLAM,CONSIG,CONNU,FV CALL TSUM(FLAP,CONLAM,CONSIG,CONNU,NOPTS,TVALUE,INVF,IER) IF(IER .EQ. 0) THEN PSERR=(FV-INVF)/DEXP(TT*CONSIG) IF(ALFA .EQ. 'R') THEN ERROR=(FV-INVF)/FV ELSE ERROR=FV-INVF ENDIF WRITE(OUTPUT,350) ID,NOPTS,INVF,PSERR,ERROR,ALFA ELSE WRITE(OUTPUT,351) ID,NOPTS,INVF,IER ENDIF IF(TVALUE .LE. 10.0) GOTO 11 NOPTS=NOPTS*ALOG10(TVALUE) CALL TSUM(FLAP,CONLAM,CONSIG,CONNU,NOPTS,TVALUE,INVF,IER) IF(IER .EQ. 0) THEN PSERR=(FV-INVF)/DEXP(TT*CONSIG) IF(ALFA .EQ. 'R') THEN ERROR=(FV-INVF)/FV ELSE ERROR=FV-INVF ENDIF WRITE(OUTPUT,350) ID,NOPTS,INVF,PSERR,ERROR,ALFA ELSE WRITE(OUTPUT,351) ID,NOPTS,INVF,IER ENDIF 11 CONTINUE 10 CONTINUE C C C LAST EXECUTABLE STATEMENT OF THE DRIVER PROGRAM C C STOP C C C C C 100 FORMAT(////' THIS IS AN EXAMPLE OF RUN FOR TESTING SOFTWARE', + '(TAPAR AND TSUM)'/' IMPLEMENTING TALBOT-METHOD FOR THE NUMERI', + 'CAL INVERSION OF LAPLACE TRANSFORMS'///) 101 FORMAT(' A SET OF 8 LAPLACE TRANSFORMS HAS BEEN CHOSEN. FOR', + ' EACH OF THEM,'/' OUTPUT SHOWS THE LOCAL ACCURACY IN APPRO', + 'XIMATING ITS INVERSE LAPLACE TRANSFORM'/' COMPUTED IN SEVERAL', + ' VALUES OF T.') 102 FORMAT(/////' FOR MORE DETAILS, SEE COMMENTS IN THE LISTING' + /////' ***** NOW THE RUN IS STARTING *****'/////) 200 FORMAT(' SINGULARITIES AND THEIR MULTIPLICITIES'/(2(2X,E12.6) + ,2X,I2)) 201 FORMAT(////' TEST FUNCTION IS LAPFUN(S)=1/S**2' + //22X,'INVLAP(T)=T'//) 202 FORMAT(////' TEST FUNCTION IS LAPFUN(S)=LOG(S)/S' + //22X,'INVLAP(T)=-LOG(T)-EULER.NUMBER'//) 203 FORMAT(////' TEST FUNCTION IS LAPFUN(S)=EXP(-4*SQRT(S))' + //22X,'INVLAP(T)=2*EXP(-4/T)/(T*SQRT(PI*T))'//) 204 FORMAT(////' TEST FUNCTION IS LAPFUN(S)=ARCTAN(1/S)' + //22X,'INVLAP(T)=SIN(T)/T'//) 205 FORMAT(////' TEST FUNCTION IS LAPFUN(S)=LOG((S**2+1)/(S**2+', + '4))'//22X,'INVLAP(T)=2*(COS(2*T)-COS(T))/T'//) 206 FORMAT(////' TEST FUNCTION IS LAPFUN(S)=S**2/(S**3+8)' + //22X,'INVLAP(T)=(EXP(-2*T)+2*EXP(T)*COS(T*SQRT(3)))/3'//) 207 FORMAT(////' TEST FUNCTION IS LAPFUN(S)=EXP(-5*S)/S' + //22X,'INVLAP(T)=HEAVISIDE(T-5)'//) 208 FORMAT(////' TEST FUNCTION IS LAPFUN(S)=EXP(-1/S)/SQRT(S)' + //22X,'INVLAP(T)=COS(2*SQRT(T))/SQRT(PI*T)'//) 250 FORMAT(//' MACHINE DECIMAL PRECISION IS ',I3//) 251 FORMAT(' INPUTS TO TAPAR ARE:'/25X,'TVALUE'/25X,'DECDIG'/25X, + 'SINGULARITIES AND THEIR MULTIPLICITIES'/25X,'MACHINE DECIMAL', + ' PRECISION'/) 252 FORMAT(' OUTPUTS FROM TAPAR ARE:'/25X,'LAMBDA'/25X,'SIGMA'/25X, + 'NU'/25X,'NO.PTS'///) 253 FORMAT(' INPUTS TO TSUM ARE:'/25X,'THE COMPLEX FUNCTION LAPFUN', + '(S)'/25X,'LAMBDA'/25X,'SIGMA'/25X,'NU'/25X,'NO.PTS'/25X, + 'TVALUE'/) 254 FORMAT(' OUTPUTS FROM TSUM ARE:'/25X,'THE APPROXIMATION TO THE', + ' INVERSE LAPLACE TRANSFORM'/25X,'AN ERROR INDICATOR'///) 300 FORMAT(///4X,'TVALUE',11X,'LAMBDA',7X,'SIGMA',8X,'NU',11X, + 'EXACT INV.TRANSFORM'/1X,E13.7,4X,3E13.7,2X,D20.14// + ' DECDIG',2X,'NO.PTS',6X,'APPROXIMATION',8X,'PSEUDO ERROR',12X, + 'ERROR') 301 FORMAT(///4X,'TVALUE',11X,'LAMBDA',7X,'SIGMA',8X,'NU',11X, + 'EXACT INV.TRANSFORM'/1X,E13.7,4X,3E13.7,2X,'DEXP(TVALUE).', + 'OVERFLOW'//' DECDIG',2X,'NO.PTS',6X,'APPROXIMATION',8X, + 'PSEUDO ERROR',12X,'ERROR') 350 FORMAT(3X,I2,3X,I7,2X,E20.14,2X,D18.12,2X,D18.12,A2) 351 FORMAT(3X,I2,3X,I7,2X,E20.14,' RETURNED VALUE, IER = ',I2) 352 FORMAT(17X,D20.14,2X,D18.12,2X,D18.12,A2) 400 FORMAT(' RERUN WITH TALBOT-MODIFICATION FOR LAMBDA AND NO.PTS'/ + ' LAMBDA = (OMEGA-1)/TVALUE+1/30'/ + ' NO.PTS = NO.PTS*LOG10(TVALUE) IF TVALUE .GT. 10'//) 401 FORMAT(' RERUN WITH OUR MODIFICATION FOR LAMBDA AND NO.PTS'/ + ' LAMBDA = (OMEGA-1)/TVALUE+1/SQRT(OMEGA*TVALUE)'/ + ' NO.PTS = NO.PTS*LOG10(TVALUE) IF TVALUE .GT. 10'//) C C C*********************************************************************** C * C * C END DRIVER * C * C * C*********************************************************************** C C END COMPLEX FUNCTION FLAP (S) C C C C********************************************************************** C C C LAPLACE TRANSFORM TEST FUNCTIONS C C C********************************************************************* C C PURPOSE C ======= C IT RETURNS THE VALUE OF ONE OF THE 8 LAPLACE TRANSFORMS PROVIDED C HERE FOR TESTING, COMPUTED IN S, DEPENDING ON THE VALUE ASSIGNED C TO THE VARIABLE NF PASSED BY A COMMON AREA. C C REMARKS C ======= C THIS FUNCTION REPRESENTS A GENERIC LAPLACE TRANSFORM SUPPLIED C BY THE USER. C C********************************************************************* C C COMMON /INVLAP/NF COMPLEX S COMPLEX A,B,I C ELIM=ALOG(R1MACH(1)) C IF(NF .EQ. 1) THEN C C F1(S)=1/S**2 C FLAP=1./(S*S) RETURN ENDIF C IF(NF .EQ. 2) THEN C C F2(S)=LOG(S)/S C FLAP=CLOG(S)/S RETURN ENDIF C IF(NF .EQ. 3) THEN C C F3(S)=EXP(-4*SQRT(S)) C A=4.*CSQRT(S) IF(REAL(A).LT.0.0) A=-A IF((REAL(A).EQ.0.0).AND.(AIMAG(A).LT.0.0)) A=-A IF (-REAL(A).LE.ELIM) THEN FLAP=(0.0,0.0) ELSE FLAP=CEXP(-A) ENDIF RETURN ENDIF C IF(NF .EQ. 4) THEN C C F4(S)=ARCTAN(1/S) C I=(0.0,1.0) A=1./S FLAP=(I-A)/(I+A) FLAP=CLOG(FLAP)/(2.*I) RETURN ENDIF C IF(NF .EQ. 5) THEN C C F5(S)=LOG((S**2+1)/(S**2+4)) C FLAP=S*S FLAP=(FLAP+1.)/(FLAP+4.) FLAP=CLOG(FLAP) RETURN ENDIF C IF(NF .EQ. 6) THEN C C F6(S)=S**2/(S**3+8) C FLAP=(1./(S+2.)+2*(S-1.)/(4+S*(S-2.)))/3. RETURN ENDIF C IF(NF .EQ. 7) THEN C C F7(S)=EXP(-5*S)/S C A=5.0*S B=S IF (-REAL(A).LE.ELIM) THEN FLAP=(0.0,0.0) ELSE FLAP=CEXP(-A)/B ENDIF RETURN ENDIF C IF(NF .EQ. 8) THEN C C F8(S)=EXP(-1/S)/SQRT(S) C A=1./S B=CSQRT(S) IF(REAL(B).LT.0.0) B=-B IF((REAL(B).EQ.0.0).AND.(AIMAG(B).LT.0.0)) B=-B IF (-REAL(A).LE.ELIM) THEN FLAP=(0.0,0.0) ELSE FLAP=CEXP(-A)/B ENDIF RETURN ENDIF C C WRITE(6,100) 100 FORMAT(//,4X,'NO TEST FUNCTION CORRESPONDS TO THIS NF-VALUE') STOP C******************************************************************** C C END FLAP C C******************************************************************** END SUBROUTINE SINGS(NDIM,SINGRE,SINGIM,MULTSI,NOSING) C C C C********************************************************************** C C C SINGULARITIES OF THE LAPLACE TRANSFORM TEST FUNCTIONS C C C********************************************************************** C C PURPOSE C ======= C FOR EACH TRANSFORM DEFINED BY NF (PASSED THROUGH A COMMON C AREA), IT RETURNS THE NUMBER OF SINGULARITIES, THEIR LOCATIONS C AND NATURE. C C********************************************************************* C C COMMON /INVLAP/ NF DIMENSION SINGRE(NDIM),SINGIM(NDIM),MULTSI(NDIM) C IF(NF .EQ. 1) THEN NOSING=1 SINGRE(1)=0. SINGIM(1)=0. MULTSI(1)=2 RETURN ENDIF C IF(NF .EQ. 2) THEN NOSING=1 SINGRE(1)=0. SINGIM(1)=0. MULTSI(1)=0 RETURN ENDIF C IF(NF .EQ. 3) THEN NOSING=1 SINGRE(1)=0. SINGIM(1)=0. MULTSI(1)=0 RETURN ENDIF C IF(NF .EQ. 4) THEN NOSING=2 SINGRE(1)=0. SINGIM(1)=0. MULTSI(1)=0 SINGRE(2)=0. SINGIM(2)=1. MULTSI(2)=0 RETURN ENDIF C IF(NF .EQ. 5) THEN NOSING=2 SINGRE(1)=0. SINGIM(1)=1. MULTSI(1)=0 SINGRE(2)=0. SINGIM(2)=2. MULTSI(2)=0 RETURN ENDIF C IF(NF .EQ. 6) THEN NOSING=2 SINGRE(1)=-2. SINGIM(1)=0. MULTSI(1)=1 SINGRE(2)=1. SINGIM(2)=SQRT(3.0) MULTSI(2)=1 RETURN ENDIF C IF(NF .EQ. 7) THEN NOSING=1 SINGRE(1)=0.0 SINGIM(1)=0.0 MULTSI(1)=1 RETURN ENDIF C IF(NF .EQ. 8) THEN NOSING=1 SINGRE(1)=0. SINGIM(1)=0. MULTSI(1)=0 RETURN ENDIF C C WRITE(6,100) 100 FORMAT(//,4X,'NO TEST FUNCTION CORRESPONDS TO THIS NF-VALUE') STOP C********************************************************************** C C END SINGS C C********************************************************************** END DOUBLE PRECISION FUNCTION FINV (T) C C C C********************************************************************** C C C INVERSE LAPLACE TRANSFORM TEST FUNCTIONS C C C********************************************************************* C C PURPOSE C ======= C IT RETURNS THE VALUE OF ONE OF THE 8 INVERSE LAPLACE TRANSFORMS C PROVIDED HERE FOR TESTING, COMPUTED IN T, DEPENDING ON THE VALUE C ASSIGNED TO THE VARIABLE NF PASSED BY A COMMON AREA. C C REMARKS C ======= C THIS FUNCTION ONLY NEEDS FOR TESTING THE LOCAL ACCURACY OF THE C TALBOT'S METHOD IMPLEMENTATION. C C********************************************************************** C C COMMON /INVLAP/ NF DOUBLE PRECISION A,B,RTRE,EULER,DELIM,D1MACH DATA EULER/0.5772156649015325D0/ C DELIM = DLOG(D1MACH(1)) C IF(NF .EQ. 1) THEN C C F1(T)=T C FINV=DBLE(T) RETURN ENDIF C IF(NF .EQ. 2) THEN C C F2(T)=-LOG(T)-EULER'S CONSTANT C A=DBLE(T) FINV=-DLOG(A)- EULER RETURN ENDIF C IF(NF .EQ. 3) THEN C C F3(T)=2*EXP(-4/T)/(T*SQRT(PI*T)) C A=4./DBLE(T) IF(-A .LE. DELIM) THEN FINV = 0.0D0 ELSE FINV = 2.*DEXP(-A) ENDIF A=4.0D0*DATAN(1.0D0)*T A=DSQRT(A) FINV=FINV/A/DBLE(T) RETURN ENDIF C IF(NF .EQ. 4) THEN C C F4(T)=SIN(T)/T C A=DBLE(T) FINV=DSIN(A)/A RETURN ENDIF C IF(NF .EQ. 5) THEN C C F5(T)=2*(COS(2*T)-COS(T))/T C A=DBLE(T) A=DCOS(A) FINV=A*(2.*A-1.)-1. FINV=2.*FINV/DBLE(T) RETURN ENDIF C IF(NF .EQ. 6) THEN C C F6(T)=(EXP(-2*T)+2*EXP(T)*COS(T*SQRT(3)))/3 C RTRE=DSQRT(3.0D0) A=DBLE(T) B=-2.*A IF(B .LE. DELIM) THEN B = 0.0D0 ELSE B = DEXP(B) ENDIF FINV=( B + 2.*DEXP(A)*DCOS(A*RTRE))/3. RETURN ENDIF C IF(NF .EQ. 7) THEN C C F7(T)=HEAVISIDE(T-5) C IF(T .EQ. 5.0) THEN FINV=0.5D0 RETURN ENDIF IF(T .GT. 5.0) THEN FINV=1.0D0 ELSE FINV=0.0D0 ENDIF RETURN ENDIF C IF(NF .EQ. 8) THEN C C F8(T)=COS(2*SQRT(T))/SQRT(PI*T) C A=DBLE(T) B=4.0*DATAN(1.0D0) A=DSQRT(A) B=DSQRT(B) FINV=DCOS(A*2.)/(A*B) RETURN ENDIF C C WRITE(6,100) 100 FORMAT(//,4X,'NO TEST FUNCTION CORRESPONDS TO THIS NF-VALUE') STOP C******************************************************************** C C END FINV C C******************************************************************** END SUBROUTINE TAPAR (TVALUE,DECDIG,NOSING,SINGRE,SINGIM,MULTSI,IC, + CONLAM,CONSIG,CONNU,NOPTS) C C C*********************************************************************** C * C * C TALBOT'S METHOD IMPLEMENTATION FOR THE NUMERICAL INVERSION OF * C * C LAPLACE TRANSFORMS * C * C (MODULE 1) * C * C * C * C * C>>>>>>>>>>>>> VERSION 3.0 FEBRUARY 1, 1989 <<<<<<<<<<* C * C * C * C * C AUTHORS: ALMERICO MURLI, MARIAROSARIA RIZZARDI * C * C * C DEPARTMENT OF MATHEMATIC AND APPLICATIONS * C UNIVERSITY OF NAPLES - ITALY * C * C * C*********************************************************************** C C C C REFERENCES C ========== C C MURLI A., M. RIZZARDI - "ALGORITHM XXX: TALBOT'S METHOD FOR THE C LAPLACE INVERSION PROBLEM". C ACM TRANS. MATH. SOFTWARE, VOL. ##, C NO. #, MONTH YEAR, PP. ##-##. C C C C PURPOSE C ======= C C THIS SUBROUTINE PROVIDES VALUES OF THE CONTOUR PARAMETERS C (LAMBDA, SIGMA, NU) C AND OF THE ACCURACY PARAMETER (N) ACCORDING TO THE TALBOT METHOD C FOR THE NUMERICAL INVERSION OF LAPLACE TRANSFORMS. C AFTER THIS ROUTINE HAS RETURNED ITS RESULTS, THESE WOULD NORMALLY C BE USED IN SUBROUTINE TSUM (MODULE 2 OF TALBOT'S METHOD IMPLEMENTA- C TION), WHICH CALCULATES THE INVERSE LAPLACE TRANSFORM F(TVALUE) C BY CONTOUR INTEGRATION. C C C C CALLING SEQUENCE C ================ C C CALL TAPAR (TVALUE, DECDIG, NOSING, SINGRE, SINGIM, MULTSI, IC, C + CONLAM, CONSIG, CONNU, NOPTS) C C C C INPUT PARAMETERS C ================ C C TVALUE - (REAL) THE REAL VALUE IN WHICH THE INVERSE LAPLACE C TRANSFORM WILL BE COMPUTED. C TVALUE MUST BE A POSITIVE NUMBER. C C DECDIG - (INTEGER) THE ACCURACY REQUIRED IN THE RESULT, IN C TERMS OF ABSOLUTE ERROR OF MAGNITUDE: C 10**(-DECDIG) IF ABS(F(TVALUE)) .LE. 1.0 C OR C DECDIG SIGNIFICANT (DECIMAL) DIGITS OTHERWISE C C NOSING - (INTEGER) THE NUMBER OF SINGULARITIES OF THE C LAPLACE TRANSFORM. C THE USER PROVIDES ONLY SINGULARITIES WITH NON C NEGATIVE IMAGINARY PARTS. THESE ARE LOCATED IN C SINGRE(J), SINGIM(J), MULTSI(J), J=1,...,NOSING. C C SINGRE - (REAL ARRAY) THIS ARRAY CONTAINS THE REAL PARTS C OF SINGULARITIES. IT MUST BE DIMENSIONED AT LEAST C "NOSING". C C SINGIM - (REAL ARRAY) THIS ARRAY CONTAINS THE IMAGINARY C PARTS OF SINGULARITIES. IT MUST BE DIMENSIONED AT C LEAST "NOSING". C C MULTSI - (INTEGER ARRAY) THIS ARRAY CONTAINS THE MULTIPLI_ C CITIES OF THOSE SINGULARITIES WHICH ARE POLES, ZERO C OTHERWISE; I.E. MULTSI(J) CONTAINS 0 IF THE J-TH C SINGULARITY IS NOT A POLE AND OTHERWISE THE MULTI_ C PLICITY OF THIS POLE. IT MUST BE DIMENSIONED LIKE C SINGRE AND SINGIM. C C IC - (INTEGER) DECIMAL MACHINE PRECISION. C IT MIGHT BE SIMPLY COMPUTED AS C IC = INT(P*LOG10(B)*3./4.) C WHERE C B IS THE RADIX OF THE COMPUTER FLOATING-POINT C SYSTEM AND P IS THE NUMBER OF SIGNIFICANT B-DIGITS C FOR MANTISSAS, OR EQUIVALENTLY, IC MAY BE COMPUTED C FOR SINGLE PRECISION COMPUTATIONS (USING THE C ROUTINE R1MACH) AS C IC = -INT(LOG10(R1MACH(3))*3./4.) C AND FOR DOUBLE PRECISION COMPUTATIONS (USING THE C ROUTINE D1MACH) AS C IC = -INT(LOG10(D1MACH(3))*3./4.) C C C C OUTPUT PARAMETERS C ================= C C CONLAM - (REAL) C C CONSIG - " C C CONNU - " GEOMETRICAL PARAMETERS FOR THE TALBOT INTE_ C GRATION CONTOUR (RESPECTIVELY LAMBDA, SIGMA AND NU) C C NOPTS - (INTEGER) THE NUMBER OF POINTS WHICH WOULD BE C REQUIRED BY THE QUADRATURE ROUTINE TO OBTAIN THE C ACCURACY INDICATED BY "DECDIG". C C C C SUBROUTINES OR FUNCTIONS NEEDED C =============================== C C INVPC : AUXILIARY SUBROUTINE OF THE ALGORITHM WHICH FINDS THE C "PRINCIPAL INVERSE" APPLYING A SUITABLE REAL NEWTON C PROCESS. C C AMAX1, AMIN1, ATAN2, DATAN, EXP, MAX0, MIN0, TAN : FORTRAN INTRINSIC C FUNCTIONS. C C C C REMARKS C ======= C C 0 IN THE FOLLOWING COMMENTS OF CODE, THE NUMBERS IN BRACKETS C REFER TO FORMULAS IN THE PAPER: C C A. TALBOT - "THE ACCURATE NUMERICAL INVERSION OF LAPLACE C TRANSFORMS". C J. INST. MATHS. APPLICS. (1979), N.23, PP.97-120. C C C 1 CHANGING THE VALUES PASSED FOR "IC", THIS ROUTINE IS COMMON C BOTH TO SINGLE AND DOUBLE PRECISION VERSION OF THE SOFTWARE C IMPLEMENTING TALBOT'S METHOD. C C C C*********************************************************************** C C LOGICAL CASO1 DOUBLE PRECISION PI INTEGER IC,DECDIG,DD DIMENSION SINGRE(NOSING),SINGIM(NOSING),MULTSI(NOSING) C DATA PI/0.0D0/ C IF(PI .EQ. 0.0D0) PI=4.*DATAN(1.0D0) CASO1=.TRUE. OMEGA=0.4*(IC+1) C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C THE FOLLOWING STATEMENTS COMPUTE C "PMAX" = MAXIMUM REAL PARTS OF SINGULARITIES (58) C "KMAX" = ARRAY-INDEX OF THE SINGULARITY HAVING THE LARGEST C MULTIPLICITY C "SIG0" = MAXIMUM BETWEEN "PMAX" AND ZERO (59) C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C SIG0=0.0 PMAX=SINGRE(1) KMAX=1 IF(NOSING .EQ. 1) GOTO 3 DO 2 J=1,NOSING IF(SINGRE(J) .GT. PMAX) PMAX=SINGRE(J) IF(MULTSI(J) .GT. MULTSI(KMAX)) KMAX=J 2 CONTINUE 3 IF(PMAX .GT. SIG0) SIG0=PMAX C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C THE FOLLOWING STATEMENTS COMPUTE THE "DOMINANT SINGULARITY" C C SD = SRD + I*SID ("I" IS THE IMAGINARY UNIT) C C SUCH THAT IM(SD)/ARG(SD) IS MAXIMUM OVER ALL THE COMPLEX C SINGULARITIES (SHIFTED BY "SIG0") WITH IM(SD) > 0. (63) C C OTHERWISE IM(SD) IS SET TO ZERO AND ARG(SD) IS SET TO "PI" C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C RD=0. SID=0. TETD=PI DO 7 J=1,NOSING IF(SINGIM(J) .LE. 0.0) GOTO 7 CASO1=.FALSE. ETAJ=ATAN2(-SINGRE(J)+SIG0, SINGIM(J)) TETJ=ETAJ + PI/2. RJ=SINGIM(J)/TETJ IF(RD .GE. RJ) GOTO 7 RD=RJ TETD=TETJ ID=J 7 CONTINUE C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C FOR AN APPROPRIATE CHOICE OF THE GEOMETRICAL PARAMETERS: C IN THE CASE OF REAL SINGULARITIES, ALWAYS "CASE-1" HOLDS C (LABEL 8). C OTHERWISE IT COMPUTES THE SHIFTED "DOMINANT SINGULARITY" C AND THE QUANTITIES "V" AND "OMEGA". (65) AND (66) C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C IF(CASO1) GOTO 8 SRD=SINGRE(ID) - SIG0 SID=SINGIM(ID) MD =MULTSI(ID) V =SID*TVALUE OMEGA=AMIN1(OMEGA + V/2., 2.*(IC+1)/3.) C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C IT CHECKS FOR "CASE-1" OR "CASE-2" C WHEN SINGULARITIES ARE COMPLEX C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C IF(1.8*V .GT. OMEGA*TETD) GOTO 9 C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C GEOMETRICAL PARAMETERS' CHOICE FOR "CASE-1" (68) C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C CASO1=.TRUE. 8 CONLAM=OMEGA/TVALUE CONSIG=SIG0 CONNU=1.0 GOTO 10 C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C GEOMETRICAL PARAMETERS' CHOICE FOR "CASE-2" (74) C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C 9 CASO1=.FALSE. PK=1.6 + 12./(V+25.) FI=1.05 + 1050./AMAX1(553.0, 800.-V) PMU=(OMEGA/TVALUE + SIG0 - PMAX)/(PK/FI - 1./TAN(FI)) CONLAM=PK*PMU/FI CONSIG=PMAX-PMU/TAN(FI) CONNU=SID/PMU C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C THE FOLLOWING STATEMENTS COMPUTE NOPTS = MAX (N0, N1, N2) (78) C C C C AT FIRST "N1" IS COMPUTED (87) C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C 10 N1=0 E=(2.3*DECDIG+OMEGA)/(CONLAM*TVALUE) IF(E .GT. 4.4) GOTO 11 UNRO=(16. + 4.3 * E)/(24.8 - 2.5 * E) GOTO 13 11 IF(E .GT. 10.) GOTO 12 UNRO=(50. + 3. * E)/(129. / E - 4.0) GOTO 13 12 UNRO=(44. + 19. * E)/(256. / E + 0.4) 13 N1=CONLAM*TVALUE*(UNRO + (CONNU-1.)/2.) N1=N1+1 N0=0 N2=0 IF(SID .EQ. 0.0) GOTO 15 C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C THE FOLLOWING STATEMENTS COMPUTE "N0" AND "N2" C FOR COMPLEX SINGULARITIES (81) AND (89) C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C DD=DECDIG+2*MIN0(MD-1, 1) + MD/4 IF(.NOT. CASO1 .OR. (MD .EQ. 0)) GOTO 14 P=SRD/CONLAM Q=SID/CONLAM CALL INVPC(P,Q,TETD,U,PI) IF(U .LE. 0.0) GOTO 14 N0=(2.3*DD+SRD*TVALUE)/U N0=N0+1 14 GAMM=(CONSIG-SIG0)/CONLAM Y=V/1000. ETA=AMIN1(1.78, 1.236+0.0064*(1.78**DD)) ETA=ETA*(1.09-Y*(0.92-0.8*Y)) N2=ETA*CONNU*(2.3*DD+OMEGA)/(3.+4.*GAMM+1./EXP(GAMM)) N2=N2+1 C C C NOPTS=MAX0(N0, N1) NOPTS=MAX0(N2, NOPTS) C C RETURN C C C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C THE FOLLOWING STATEMENTS COMPUTE "N0" AND "N2" C FOR REAL SINGULARITIES (82) AND (89) C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C 15 DD=DECDIG+2*MIN0(MULTSI(KMAX)-1, 1) + MULTSI(KMAX)/4 N2=(2.3*DD+OMEGA)/2. N2=N2+1 ICM1=IC-1 IF(DD .LT. ICM1) GOTO 17 Q=0.0 DO 16 L=1,NOSING IF(MULTSI(L) .EQ. 0) GOTO 16 DD=DECDIG+2*MIN0(MULTSI(L)-1, 1) + MULTSI(L)/4 SRD=SINGRE(L) - SIG0 IF(.NOT. (DD .GE. ICM1 .AND. SRD .LT. 0.0)) GOTO 16 P=SRD/CONLAM CALL INVPC(P,Q,TETD,U,PI) IF(U .LE. 0.0) GOTO 16 N0L=(2.3*DD+TVALUE*SRD)/U IF(N0L .GT. N0) N0=N0L 16 CONTINUE N0=N0+1 C C C 17 NOPTS=MAX0(N0, N1) NOPTS=MAX0(N2, NOPTS) C C RETURN C C C C*********************************************************************** C * C END TAPAR * C * C*********************************************************************** C END SUBROUTINE TSUM (FLAP,CONLAM,CONSIG,CONNU,NOPTS,TVALUE,FINV,IER) C C C*********************************************************************** C * C * C TALBOT'S METHOD IMPLEMENTATION FOR THE NUMERICAL INVERSION OF * C * C LAPLACE TRANSFORMS * C * C (MODULE 2) * C * C * C * C * C>>>>>>>>>>>>> VERSION 3.0 FEBRUARY 1, 1989 <<<<<<<<<<* C * C * C * C * C AUTHORS: ALMERICO MURLI, MARIAROSARIA RIZZARDI * C * C * C DEPARTMENT OF MATHEMATIC AND APPLICATIONS * C UNIVERSITY OF NAPLES - ITALY * C * C * C*********************************************************************** C C C C REFERENCES C ========== C C MURLI A., M. RIZZARDI - "ALGORITHM XXX: TALBOT'S METHOD FOR THE C LAPLACE INVERSION PROBLEM". C ACM TRANS. MATH. SOFTWARE, VOL. ##, C NO. #, MONTH YEAR, PP. ##-##. C C C C PURPOSE C ======= C C THIS SUBROUTINE PROVIDES A NUMERICAL APPROXIMATION, ACCORDING TO C TALBOT'S METHOD, OF THE INVERSE LAPLACE TRANSFORM. C THIS ROUTINE APPLIES THE TRAPEZOIDAL RULE FOR THE APPROXIMATION OF C THE CONTOUR INTEGRAL GIVING THE INVERSE LAPLACE TRANSFORM. C IN THE QUADRATURE RULE, THE SUM IS ACCOMPLISHED USING THE REINSCH C ALGORITHM. C SINGLE PRECISION ARITHMETIC IS USED. C THE GEOMETRICAL PARAMETERS (LAMBDA, SIGMA, NU) INDIVIDUATE A PARTI_ C CULAR CONTOUR AND THE ACCURACY PARAMETER (N) GIVES THE NUMBER OF C TERMS IN THE SUM. THESE PARAMETERS MAY BE PREVIOUSLY COMPUTED BY A C CALL TO THE SUBROUTINE "TAPAR" OR THEY MAY BE ASSIGNED BY AN EXPERT C USER. C C C C CALLING SEQUENCE C ================ C C CALL TSUM (FLAP, CONLAM, CONSIG, CONNU, NOPTS, TVALUE, FINV, IER) C C C C INPUT PARAMETERS C ================ C C FLAP - (COMPLEX FUNCTION) THIS FUNCTION, USER SUPPLIED, C REPRESENTS THE LAPLACE TRANSFORM TO BE INVERTED. C IT MUST BE DEFINED BY AN EXTERNAL STATEMENT IN THE C CALLING PROGRAM AND IT MUST BE SUPPLIED AS A C FUNCTION LIKE: C COMPLEX FUNCTION FLAP (S) C COMPLEX S C ..... C END C C CONLAM - (REAL) C CONSIG - " C CONNU - " THE GEOMETRICAL PARAMETERS OF THE INTEGRAL C CONTOUR (RESPECTIVELY LAMBDA, SIGMA AND NU). THEY C MAY BE PREVIOUSLY COMPUTED BY THE SUBROUTINE C "TAPAR" ACCORDING TO TALBOT'S METHOD OR SUPPLIED C BY THE USER. C C NOPTS - (INTEGER) THE ACCURACY PARAMETER (N), I.E. THE C NUMBER OF TERMS IN THE SUM THAT GIVES THE RESULT. C IT ALSO REPRESENTS THE NUMBER OF FUNCTION EVALUA_ C TIONS. IT MAY BE PREVIOUSLY COMPUTED BY THE C SUBROUTINE "TAPAR" ACCORDING TO TALBOT'S METHOD OR C SUPPLIED BY THE USER. C C TVALUE - (REAL) THE VALUE IN WHICH THE INVERSE TRANSFORM C WILL BE APPROXIMATED. C C C C OUTPUT PARAMETERS C ================= C C FINV - (REAL) THE NUMERICAL APPROXIMATION TO THE INVERSE C LAPLACE TRANSFORM COMPUTED IN "TVALUE". C C IER - (INTEGER) ERROR INDICATOR. C IER = 0 NORMAL TERMINATION. C IER = 1 THE OUTPUT VALUE OF "FINV" CAUSES OVER_ C FLOW. IN THIS CASE THE ROUTINE RETURNS A C SCALED VALUE IN "FINV" SUCH THAT THE C FINAL RESULT MAY BE COMPUTED AS: C C F(T) = EXP(SIGMA*T)*CONLAM/NOPTS*FINV C C C C SUBROUTINES OR FUNCTIONS NEEDED C =============================== C C FLAP : COMPLEX FUNCTION (INPUT PARAMETER) FOR THE LAPLACE C TRANSFORM TO BE INVERTED. C C R1MACH : FUNCTION FOR SINGLE PRECISION MACHINE DEPENDENT C CONSTANTS. C C XERRWV, XERROR : SUBROUTINES (FROM SLATEC LIBRARY) TO PROCESS C A DIAGNOSTIC MESSAGE. C IF THESE ROUTINES ARE NOT AVAILABLE, THE CODE ALSO C CONTAINS STATEMENTS TO PRINT THE OUTPUT MESSAGE. C IT SUFFICES TO REMOVE THE C FROM COLUMN 1 TO ACTIVATE C THE DESIRED CHOICE. C XERRWV CALL IS ACTUALLY ACTIVATED. C C ABS, AIMAG, ALOG, CMPLX, COS, DATAN, EXP, FLOAT, REAL, SIN : C FORTRAN INTRINSIC FUNCTIONS. C C C C*********************************************************************** C C EXTERNAL FLAP DOUBLE PRECISION PI COMPLEX FLAP, S, FF REAL CONLAM, CONSIG, CONNU, R1MACH C C DATA PI/0.0D0/ C C ELIM=ALOG(R1MACH(1)) IF(PI .EQ. 0.0D0) PI=4.*DATAN(1.0D0) PIOVN=PI/FLOAT(NOPTS) IER=0 TAU=CONLAM*TVALUE PSI=PIOVN*TAU*CONNU C=COS(PSI) BR=0. BI=0. DBR=0. DBI=0. NM1=NOPTS - 1 C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C COMPUTES ONLY T0 FOR NOPTS = 1 C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C FINV=0. IF(NM1 .EQ. 0) GOTO 5 C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C BEGIN REINSCH ALGORITHM C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C IF(C .GT. 0.) GOTO 2 C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C COS(PSI) <= 0 C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C U=4.*COS(PSI/2)**2 DO 1 K=NM1, 1, -1 TETA=K*PIOVN ALFA=TETA*COS(TETA)/SIN(TETA) BETA=TETA+ALFA*(ALFA-1.)/TETA SR =CONLAM*ALFA+CONSIG SI =CONLAM*CONNU*TETA S =CMPLX(SR,SI) FF =FLAP(S) G =REAL(FF) H =AIMAG(FF) C EAT=EXP(ALFA*TAU) ARG =ALFA*TAU IF (ARG.LE.ELIM) THEN EAT =0.0 ELSE EAT =EXP(ARG) ENDIF BR =DBR - BR BI =DBI - BI DBR =U*BR - DBR + EAT*(G*CONNU - H*BETA) DBI =U*BI - DBI + EAT*(H*CONNU + G*BETA) C 1 CONTINUE C C BR =DBR - BR BI =DBI - BI DBR =U*BR - DBR C C C GOTO 4 C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C COS(PSI) > 0 C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C 2 CONTINUE U=-4.*SIN(PSI/2)**2 DO 3 K=NM1, 1, -1 TETA=K*PIOVN ALFA=TETA*COS(TETA)/SIN(TETA) BETA=TETA+ALFA*(ALFA-1.)/TETA SR =CONLAM*ALFA+CONSIG SI =CONLAM*CONNU*TETA S =CMPLX(SR,SI) FF =FLAP(S) G =REAL(FF) H =AIMAG(FF) C EAT=EXP(ALFA*TAU) ARG =ALFA*TAU IF (ARG.LE.ELIM) THEN EAT=0.0 ELSE EAT =EXP(ARG) ENDIF BR =DBR+BR BI =DBI+BI DBR =U*BR + DBR + EAT*(G*CONNU - H*BETA) DBI =U*BI + DBI + EAT*(H*CONNU + G*BETA) C 3 CONTINUE C C BR =DBR + BR BI =DBI + BI DBR =U*BR + DBR C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C END REINSCH ALGORITHM C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C 4 FINV=DBR - BR*U/2 - BI*SIN(PSI) C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C COMPUTES T0 C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C 5 SR =CONLAM + CONSIG S =CMPLX(SR,0.0) FF =FLAP(S) G =REAL(FF) FINV=(FINV+CONNU*EXP(TAU)*G/2.) IF(FINV .EQ. 0.0) RETURN IF(FINV .GT. 0.0) THEN SIGNUM=1.0 ELSE SIGNUM=-1.0 ENDIF C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C OVERFLOW PROBLEM HANDLING C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C COST= ALOG ( R1MACH(2) ) OF =CONSIG*TVALUE + ALOG(CONLAM*ABS(FINV)/NOPTS) IF(OF .GT. COST) THEN IER = 1 C C DIAGNOSTIC MESSAGE C CALL XERRWV(' ABNORMAL RETURN FROM TSUM (POSSIBILITY OF OVERF +LOW), ERROR NUMBER = ',69,1,0,1,IER,0,0,0.,0.) C C PRINT 999, IER C999 FORMAT(/' ABNORMAL RETURN FROM TSUM'/' ERROR NUMBER = ', C + I2,' (POSSIBILITY OF OVERFLOW)'/) C ELSE FINV=EXP(OF)*SIGNUM ENDIF C C RETURN C C C C*********************************************************************** C * C END TSUM * C * C*********************************************************************** C END SUBROUTINE INVPC(P,Q,TETA,U,PI ) C C C C********************************************************************** C * C TALBOT'S METHOD IMPLEMENTATION FOR THE NUMERICAL INVERSION OF * C * C LAPLACE TRANSFORMS * C * C (MODULE 3 - AUXILIARY ROUTINE) * C * C * C * C * C * C>>>>>>>>>>>>> VERSION 3.0 FEBRUARY 1, 1989 <<<<<<<<<<* C * C * C * C * C AUTHORS: ALMERICO MURLI, MARIAROSARIA RIZZARDI * C * C * C DEPARTMENT OF MATHEMATIC AND APPLICATIONS * C UNIVERSITY OF NAPLES - ITALY * C * C * C********************************************************************** C C PURPOSE C ======= C FIND THE PRINCIPAL INVERSE Z=-U + I*Y OF S=P + I*Q=R*EXP(I*TETA) C IN S=Z/(1 - EXP(-Z)) C APPLYING NEWTON'S METHOD FOR REAL ROOTS C C*********************************************************************** C C C ROUTINE PARAMETERS DESCRIPTION C ============================== C C INPUT PARAMETERS C ================ C P - (REAL) C Q - C TETA - RESPECTIVELY REAL PART, IMAGINARY PART AND C ARGUMENT OF S. C IT IS SUPPOSED THAT C P .LE. 0.0 C Q .GE. 0.0 C PI/2 .LE. TETA .LE. PI C C PI - (DOUBLE PRECISION) THE CONSTANT PI ALREADY COMPUTE C IN THE CALLING PROGRAM. C C C OUTPUT PARAMETERS C ================= C U - (REAL) AN APPROXIMATION TO THE REAL PART (SIGN INVER C TED) OF THE PRINCIPAL INVERSE OF S. C C C C*********************************************************************** C DOUBLE PRECISION PI C R=P*P+Q*Q R=SQRT(R) C C SET THE INITIAL GUESS C X=13./(5.-2*P-Q-0.45*EXP(P)) Y=2*PI -X K=1 157 CONTINUE ANG=Y-TETA S=SIN(ANG) C=COS(ANG) B=Q-Y E=-B/S ARG=E/R IF(ARG .LE. 0.0) THEN U=-1.0 GOTO 200 ENDIF U=ALOG(ARG) G=P+B*C/S+U DY=B*G/(1.+E*(E-2.*C)) C C TEST FOR CONVERGENCE (EPS IS 1.E-4) C IF(ABS(DY) .LT. 0.0001) GOTO 200 C C TEST FOR NUMBER OF ITERATION (MAXIMUM IS 100) C IF(K .GE. 100) GOTO 200 K=K+1 Y=Y+DY GOTO 157 200 CONTINUE C RETURN C********************************************************************* C C END INVPC C C********************************************************************* END