home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 2 / RISC_DISC_2.iso / pd_share / program / fortran77_210 / acmtoms / f77 / invlap619 < prev    next >
Encoding:
Text File  |  1992-01-17  |  16.9 KB  |  528 lines
  1. ***********************  COPYRIGHT NOTICE  ********************************
  2. *  The material in this library is copyrighted by the ACM, which grants   *
  3. *  general permission to distribute provided the copies are not made for  *
  4. *  direct commercial advantage.  For details of the copyright and         *
  5. *  dissemination agreement, consult the current issue of TOMS.            *
  6. ***************************************************************************
  7. CDate:           Thu, 17 Oct 91  11:21 BST
  8. CFrom:           CHECAJ@UK.AC.HERIOT-WATT.VAXB
  9. CTo:             KMC@UK.AC.RUTHERFORD.DEC-E
  10. CSubject:        fortran - inverse laplace transform
  11. C
  12. C     ALGORITHM 619, COLLECTED ALGORITHMS FROM ACM.
  13. C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 3,
  14. C     SEP., 1984, P. 348-353.
  15. C   DRIVER PROGRAM FOR DLAINV
  16. C   THE INVERSE LAPLACE TRANSFORM OF 1/(P**2+1) IS COMPUTED
  17. C   FOR T=0.1,1,2,3,4,5,10,20,30,40,50,60,70,80,90 AND 100
  18. C   (THESE VALUES ARE STORED IN THE ARRAY TVAL).
  19. C   THE REQUESTED TOLERANCES ARE EPSAB=EPSRE=1.0D-4, 1.0D-8
  20. C   AND 1.0D-12 (THESE VALUES ARE STORED IN THE ARRAY E)
  21. C   THE EXACT INVERSE LAPLACE TRANSFORM IS SIN(T).
  22. C   ALSO THE EXACT ERROR IS COMPUTED
  23. C
  24.       DOUBLE PRECISION C,DIF,E(3),ESTERR,EXA,RESULT,T,TVAL(16)
  25.       DOUBLE PRECISION DABS,DSIN
  26.       EXTERNAL FUN
  27. C
  28. C   THE ARRAY E CONTAINS THE REQUESTED TOLERANCES
  29. C
  30.       DATA E(1),E(2),E(3)/1.0D-4,1.0D-8,1.0D-12/
  31. C
  32. C   THE ARRAY TVAL CONTAINS THE VALUES OF THE INDEPENDENT
  33. C   VARIABLE OF THE INVERSE LAPLACE TRANSFORM
  34. C
  35.       DATA TVAL(1),TVAL(2),TVAL(3),TVAL(4),TVAL(5),TVAL(6),
  36.      *  TVAL(7),TVAL(8),TVAL(9),TVAL(10),TVAL(11),TVAL(12),
  37.      *  TVAL(13),TVAL(14),TVAL(15),TVAL(16) /0.1D0,1.0D0,2.0D0,
  38.      *  3.0D0,4.0D0,5.0D0,1.0D1,2.0D1,3.0D1,4.0D1,5.0D1,6.0D1,
  39.      *  7.0D1,8.0D1,9.0D1,1.0D2/
  40. C
  41. C   C IS THE ABSCISSA OF CONVERGENCE
  42. C
  43.       C = 0.0D0
  44.       WRITE(6,99997)
  45.       DO 20 I=1,3
  46.         WRITE(6,99998) E(I)
  47.         DO 10 K=1,16
  48.           T = TVAL(K)
  49.         WRITE(*,*)' T= ',T
  50. C
  51. C   COMPUTATION OF THE APPROXIMATE INVERSE LAPLACE TRANSFORM
  52. C
  53.           CALL DLAINV(FUN,T,C,E(I),E(I),RESULT,ESTERR,NUM,IER)
  54. C
  55. C   COMPUTATION OF THE EXACT INVERSE LAPLACE TRANSFORM
  56. C
  57.           EXA = DSIN(T)
  58. C
  59. C   COMPUTATION OF THE EXACT ERROR
  60. C
  61.           DIF = DABS(EXA-RESULT)
  62.           WRITE(6,99999) T,EXA,RESULT,DIF,ESTERR,NUM,IER
  63.    10   CONTINUE
  64.    20 CONTINUE
  65. 99997 FORMAT(1H0,4X,1HT,5X,1H*,8X,11HEXACT VALUE,13X,8HCOMPUTED,
  66.      *  6H VALUE,9X,5HERROR,7X,6HESTERR,4X,3HNUM,4X,3HIER//1H ,
  67.      *  98(1H*)/)
  68. 99998 FORMAT(1H0,16HEPSAB = EPSRE = ,D9.2/)
  69. 99999 FORMAT(5X,F5.1,1X,1H*,1X,2(D23.16,3X),2(D9.2,3X),I3,I6)
  70.       STOP
  71.       END
  72.        SUBROUTINE DLAINV(FUN,T,C,EPSRE,EPSAB,RESULT,ESTERR,NUM,
  73.      *   IER)
  74. C
  75. C ......................................................................
  76. C
  77. C   1. DLAINV
  78. C        INVERSION OF LAPLACE TRANSFORM USING THE DURBIN FORMULA
  79. C        IN COMBINATION WITH THE EPSILON ALGORITHM
  80. C
  81. C   2. PURPOSE
  82. C        THE ROUTINE CALCULATES AN APPROXIMATE RESULT TO THE
  83. C        INVERSE LAPLACE TRANSFORM F(T) OF FUN, FOR THE VALUE
  84. C        OF THE INDEPENDENT VARIABLE EQUAL TO T, HOPEFULLY
  85. C        SATISFYING FOLLOWING CLAIM FOR ACCURACY :
  86. C        ABS(F(T)-RESULT).LE.MAX(EPSAB,EPSR*ABS(F(T)))
  87. C
  88. C   3. CALLING SEQUENCE
  89. C        CALL DLAINV(FUN,T,C,EPSRE,EPSAB,RESULT,ESTERR,NUM,IER)
  90. C
  91. C      INPUT PARAMETERS
  92. C        FUN    - DOUBLE PRECISION
  93. C                 SUBROUTINE DEFINING THE LAPLACE TRANSFORM AS
  94. C                 A COMPLEX FUNCTION.  THE CALLING SEQUENCE OF
  95. C                 FUN IS : CALL FUN(A,B,C,D) WHERE
  96. C                 A : DOUBLE PRECISION
  97. C                     REAL PART OF THE INDEPENDENT VARIABLE
  98. C                     OF THE LAPLACE TRANSFORM (INPUT)
  99. C                 B : DOUBLE PRECISION
  100. C                     IMAGINARY PART OF THE INDEPENDENT
  101. C                     VARIABLE OF THE LAPLACE TRANSFORM (INPUT)
  102. C                 C : DOUBLE PRECISION
  103. C                     REAL PART OF THE VALUE OF THE LAPLACE
  104. C                     TRANSFORM (OUTPUT)
  105. C                 D : DOUBLE PRECISION
  106. C                     IMAGINARY PART OF THE VALUE OF THE
  107. C                     LAPLACE TRANSFORM (OUTPUT)
  108. C                 THE ACTUAL NAME FOR FUN NEEDS TO BE DECLARED
  109. C                 EXTERNAL IN THE DRIVER PROGRAM.
  110. C
  111. C        T      - DOUBLE PRECISION
  112. C                 VALUE OF THE INDEPENDENT VARIABLE FOR WHICH THE
  113. C                 INVERSE LAPLACE TRANSFORM HAS TO BE COMPUTED.
  114. C                 T SHOULD BE GREATER THAN ZERO.
  115. C
  116. C        C      - DOUBLE PRECISION
  117. C                 ABSCISSA OF CONVERGENCE OF THE LAPLACE TRANSFORM
  118. C
  119. C        EPSRE  - DOUBLE PRECISION
  120. C                 RELATIVE ACCURACY REQUESTED
  121. C
  122. C        EPSAB  - DOUBLE PRECISION
  123. C                 ABSOLUTE ACCURACY REQUESTED.
  124. C                 THE ROUTINE TRIES TO SATISFY THE LEAST STRINGENT
  125. C                 OF BOTH ACCURACY REQUIREMENT.
  126. C
  127. C      OUTPUT PARAMETERS
  128. C        RESULT - DOUBLE PRECISION
  129. C                 INVERSE LAPLACE TRANSFORM
  130. C
  131. C        ESTERR - DOUBLE PRECISION
  132. C                 ESTIMATE OF THE ABSOLUTE ERROR ABS(F(T)-RESULT)
  133. C
  134. C        NUM    - INTEGER
  135. C                 NUMBER OF EVALUATIONS OF FUN
  136. C
  137. C        IER    - INTEGER
  138. C                 PARAMETER GIVING INFORMATION ON THE TERMINATION
  139. C                 OF THE ALGORITHM
  140. C                 IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  141. C                         ROUTINE
  142. C                 IER = 1 THE COMPUTATIONS ARE TERMINATED BECAUSE
  143. C                         THE BOUND ON THE NUMBER OF EVALUATIONS
  144. C                         OF FUN HAS BEEN ACHIEVED.  THIS BOUND
  145. C                         IS EQUAL TO 8*MAX+5 WHERE  MAX  IS A
  146. C                         NUMBER INITIALIZED IN A DATA
  147. C                         STATEMENT.  ONE CAN ALLOW MORE FUNCTION
  148. C                         EVALUATIONS BY INCREASING THE VALUE OF
  149. C                         MAX IN THE DATA-STATEMENT.
  150. C                 IER = 2 THE VALUE OF T IS LESS THAN OR EQUAL
  151. C                         TO ZERO.
  152. C
  153. C   4. SUBROUTINES OR FUNCTIONS NEEDED
  154. C               - FUN : USER PROVIDED SUBROUTINE
  155. C               - DQEXT : EPSILON ALGORITHM
  156. C               - D1MACH : THIS SUBPROGRAM IS CALLED BY
  157. C                          DQEXT, AND PROVIDES MACHINE CONSTANTS
  158. C               - DABS, DATAN, DEXP, DMAX1, DSIN : FORTRAN FUNCTIONS
  159. C
  160. C ......................................................................
  161. C
  162.       DOUBLE PRECISION AIM,AK,ARE,ARG,BB,C,DABS,DATAN,DEXP,DMAX1,
  163.      *  DSIN,EPSAB,EPSRE,ESTERR,FIM,FRE,PID16,R,RESULT,RES3LA(3),
  164.      *  REX(52),SI(32),T
  165.       INTEGER I,IER,K,KC,KK,KS,M,NEX,NRES,NUM
  166.       DATA REX /52*0.0D0/
  167. C
  168. C   THE ARRAY SI CONTAINS VALUES OF THE SINE AND COSINE FUNCTIONS
  169. C   REQUIRED IN THE DURBIN FORMULA. SI(8) AND SI(16) ARE GIVEN IN
  170. C   THE FOLLOWING DATA STATEMENT. THE OTHER VALUES ARE COMPUTED.
  171. C
  172.       DATA SI(8),SI(16)/ 1.0D+00,0.0D+00/
  173. C
  174. C   MAX IS A BOUND ON THE NUMBER OF TERMS USED IN THE DURBIN
  175. C   FORMULA.
  176. C
  177.       DATA MAX/500/
  178. C
  179. C   TEST ON VALIDITY OF THE INPUT PARAMETER T
  180. C
  181.       IER = 2
  182.       RESULT = 0.0D+00
  183.       ESTERR = 1.0D+00
  184.       NUM = 0
  185.       IF (T.LE.0.0D+00) GO TO 999
  186. C
  187. C   PID16 IS EQUAL TO PI/16
  188. C
  189.       PID16 = DATAN(1.0D+00)/4.0D+00
  190. C
  191. C   COMPUTATION OF THE ELEMENTS OF THE ARRAY SI
  192. C
  193.       AK = 1.0D+00
  194.       DO 10 K=1,7
  195.         SI(K) = DSIN(AK*PID16)
  196.         AK = AK+1.0D+00
  197.         KK = 16-K
  198.         SI(KK) = SI(K)
  199. 10    CONTINUE
  200.       IER = 0
  201.       NRES = 0
  202.       DO 20 K=17,32
  203.         SI(K) = -SI(K-16)
  204. 20    CONTINUE
  205. C
  206. C   INITIALIZATION OF THE SUMMATION OF THE DURBIN FORMULA.
  207. C
  208.       ARG = PID16/T
  209.       ARE = C+2.0D+00/T
  210.       AIM = 0.0D+00
  211.       BB = DEXP(ARE*T)/(1.6D+01*T)
  212.       CALL FUN (ARE,AIM,FRE,FIM)
  213.       NUM = 5
  214.       R = 5.0D-01*FRE
  215.       NEX = 0
  216.       KC = 8
  217.       KS = 0
  218. C
  219. C   MAIN LOOP FOR THE SUMMATION
  220. C
  221.       DO 40 I=1,MAX
  222.         M = 8
  223.         IF (I.EQ.1) M = 12
  224.         DO 30 K=1,M
  225.           AIM = AIM+ARG
  226.           KC = KC+1
  227.           KS = KS+1
  228.           IF (KC.GT.32) KC = 1
  229.           IF (KS.GT.32) KS = 1
  230.           CALL FUN(ARE,AIM,FRE,FIM)
  231.           R = R+FRE*SI(KC)-FIM*SI(KS)
  232. 30      CONTINUE
  233.         NUM = NUM+8
  234.         NEX = NEX+1
  235.         REX(NEX) = R
  236. C
  237. C   EXTRAPOLATION USING THE EPSILON ALGORITHM
  238. C
  239.         IF(NEX.GE.3) CALL DQEXT(NEX,REX,RESULT,ESTERR,RES3LA,NRES)
  240.         IF(NRES.LT.4) GO TO 40
  241. C
  242. C   COMPUTATION OF INTERMEDIATE RESULT AND ESTIMATE OF THE
  243. C   ABSOLUTE ERROR
  244. C
  245.         RESULT = RESULT * BB
  246.         ESTERR = ESTERR * BB
  247.         IF (ESTERR.LT.DMAX1(EPSAB,EPSRE*DABS(RESULT)).AND.DABS(R*BB-
  248.      *   RESULT).LT.5.0D-01*DABS(RESULT)) GO TO 999
  249. 40    CONTINUE
  250. C
  251. C   SET ERROR FLAG IN THE CASE THAT THE NUMBER OF TERMS IN THE
  252. C   SUMMATION IS EQUAL TO MAX
  253. C
  254.       IER = 1
  255. 999   RETURN
  256.       END
  257.       SUBROUTINE DQEXT(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES)
  258. C
  259. C 1.        DQEXT
  260. C           EPSILON ALGORITHM
  261. C              STANDARD FORTRAN SUBROUTINE
  262. C              DOUBLE PRECISION VERSION
  263. C
  264. C 2.        PURPOSE
  265. C              THE ROUTINE DETERMINES THE LIMIT OF A GIVEN SEQUENCE OF
  266. C              APPROXIMATIONS, BY MEANS OF THE EPSILON ALGORITHM
  267. C              OF P. WYNN.
  268. C              AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN.
  269. C              THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE
  270. C              ELEMENTS NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL
  271. C              ARE PRESERVED.
  272. C
  273. C 3.        CALLING SEQUENCE
  274. C              CALL DQEXT(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES)
  275. C
  276. C           PARAMETERS
  277. C              N      - INTEGER
  278. C                       EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE
  279. C                       FIRST COLUMN OF THE EPSILON TABLE.
  280. C
  281. C              EPSTAB - DOUBLE PRECISION
  282. C                       VECTOR OF DIMENSION 52 CONTAINING THE ELEMENTS
  283. C                       OF THE TWO LOWER DIAGONALS OF THE
  284. C                       TRIANGULAR EPSILON TABLE
  285. C                       THE ELEMENTS ARE NUMBERED STARTING AT THE
  286. C                       RIGHT-HAND CORNER OF THE TRIANGLE.
  287. C
  288. C              RESULT - DOUBLE PRECISION
  289. C                       RESULTING APPROXIMATION TO THE INTEGRAL
  290. C
  291. C              ABSERR - DOUBLE PRECISION
  292. C                       ESTIMATE OF THE ABSOLUTE ERROR COMPUTED FROM
  293. C                       RESULT AND THE 3 PREVIOUS RESULTS
  294. C
  295. C              RES3LA - DOUBLE PRECISION
  296. C                       VECTOR OF DIMENSION 3 CONTAINING THE LAST 3
  297. C                       RESULTS
  298. C
  299. C              NRES   - INTEGER
  300. C                       NUMBER OF CALLS TO THE ROUTINE
  301. C                       (SHOULD BE ZERO AT FIRST CALL)
  302. C
  303. C 4.        SUBROUTINES OR FUNCTIONS NEEDED
  304. C                     - D1MACH
  305. C                     - FORTRAN DABS, DMAX1
  306. C
  307. C     ..................................................................
  308. C
  309.       DOUBLE PRECISION ABSERR,D1MACH,DABS,DELTA1,DELTA2,DELTA3,DMAX1,
  310.      *  EPMACH,EPSINF,EPSTAB(52),ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3,
  311.      *  OFLOW,RES,RESULT,RES3LA(3),SS,TOL1,TOL2,TOL3
  312.       INTEGER I,IB,IB2,IE,INDX,K1,K2,K3,LIMEXP,N,NEWELM,NRES,NUM
  313. C
  314. C           LIST OF MAJOR VARIABLES
  315. C           -----------------------
  316. C
  317. C           E0     - THE 4 ELEMENTS ON WHICH THE
  318. C           E1       COMPUTATION OF A NEW ELEMENT IN
  319. C           E2       THE EPSILON TABLE IS BASED
  320. C           E3                 E0
  321. C                        E3    E1    NEW
  322. C                              E2
  323. C           NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW
  324. C                    DIAGONAL
  325. C           ERROR  - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2)
  326. C           RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE
  327. C                    OF ERROR
  328. C
  329. C           MACHINE DEPENDENT CONSTANTS
  330. C           ---------------------------
  331. C
  332. C           EPMACH IS THE LARGEST RELATIVE SPACING.
  333. C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
  334. C           LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON
  335. C           TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER
  336. C           DIAGONAL OF THE EPSILON TABLE IS DELETED.
  337. C
  338. C***FIRST EXECUTABLE STATEMENT
  339.       OFLOW = D1MACH(2)
  340.       EPMACH = D1MACH(4)
  341.       NRES = NRES+1
  342.       ABSERR = OFLOW
  343.  
  344.       RESULT = EPSTAB(N)
  345.       IF(N.LT.3) GO TO 100
  346.       LIMEXP = 50
  347.       EPSTAB(N+2) = EPSTAB(N)
  348.       NEWELM = (N-1)/2
  349.       EPSTAB(N) = OFLOW
  350.       NUM = N
  351.       K1 = N
  352.       DO 40 I = 1,NEWELM
  353.         K2 = K1-1
  354.         K3 = K1-2
  355.         RES = EPSTAB(K1+2)
  356.         E0 = EPSTAB(K3)
  357.         E1 = EPSTAB(K2)
  358.         E2 = RES
  359.         E1ABS = DABS(E1)
  360.         DELTA2 = E2-E1
  361.         ERR2 = DABS(DELTA2)
  362.         TOL2 = DMAX1(DABS(E2),E1ABS)*EPMACH
  363.         DELTA3 = E1-E0
  364.         ERR3 = DABS(DELTA3)
  365.         TOL3 = DMAX1(E1ABS,DABS(E0))*EPMACH
  366.         IF(ERR2.GT.TOL2.OR.ERR3.GT.TOL3) GO TO 10
  367. C
  368. C           IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
  369. C           ACCURACY, CONVERGENCE IS ASSUMED.
  370. C           RESULT = E2
  371. C           ABSERR = ABS(E1-E0)+ABS(E2-E1)
  372. C
  373.         RESULT = RES
  374.         ABSERR = ERR2+ERR3
  375. C***JUMP OUT OF DO-LOOP
  376.         GO TO 100
  377.    10   CONTINUE
  378.         E3 = EPSTAB(K1)
  379.         EPSTAB(K1) = E1
  380.         DELTA1 = E1-E3
  381.         ERR1 = DABS(DELTA1)
  382.         TOL1 = DMAX1(E1ABS,DABS(E3))*EPMACH
  383. C
  384. C           IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
  385. C           A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
  386. C
  387.         IF(ERR1.LE.TOL1.OR.ERR2.LE.TOL2.OR.ERR3.LE.TOL3) GO TO 20
  388.         SS = 1.0D+00/DELTA1+1.0D+00/DELTA2-1.0D+00/DELTA3
  389.         EPSINF = DABS(SS*E1)
  390. C
  391. C           TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND
  392. C           EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE
  393. C           OF N.
  394. C
  395.         IF(EPSINF.GT.1.0D-04) GO TO 30
  396.    20   N = I+I-1
  397. C***JUMP OUT OF DO-LOOP
  398.         GO TO 50
  399. C
  400. C           COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST
  401. C           THE VALUE OF RESULT.
  402. C
  403.    30   RES = E1+1.0D+00/SS
  404.         EPSTAB(K1) = RES
  405.         K1 = K1-2
  406.         ERROR = ERR2+DABS(RES-E2)+ERR3
  407.         IF(ERROR.GT.ABSERR) GO TO 40
  408.         ABSERR = ERROR
  409.         RESULT = RES
  410.    40 CONTINUE
  411. C
  412. C           SHIFT THE TABLE.
  413. C
  414.    50 IF(N.EQ.LIMEXP) N = 2*(LIMEXP/2)-1
  415.       IB = 1
  416.       IF((NUM/2)*2.EQ.NUM) IB = 2
  417.       IE = NEWELM+1
  418.       DO 60 I=1,IE
  419.         IB2 = IB+2
  420.         EPSTAB(IB) = EPSTAB(IB2)
  421.         IB = IB2
  422.    60 CONTINUE
  423.       IF(NUM.EQ.N) GO TO 80
  424.       INDX = NUM-N+1
  425.       DO 70 I = 1,N
  426.         EPSTAB(I)= EPSTAB(INDX)
  427.         INDX = INDX+1
  428.    70 CONTINUE
  429.    80 IF(NRES.GE.4) GO TO 90
  430.       RES3LA(NRES) = RESULT
  431.       ABSERR = OFLOW
  432.       GO TO 100
  433. C
  434. C           COMPUTE ERROR ESTIMATE
  435. C
  436.    90 ABSERR = DABS(RESULT-RES3LA(3))+DABS(RESULT-RES3LA(2))
  437.      *  +DABS(RESULT-RES3LA(1))
  438.       RES3LA(1) = RES3LA(2)
  439.       RES3LA(2) = RES3LA(3)
  440.       RES3LA(3) = RESULT
  441.   100 ABSERR = DMAX1(ABSERR,5.0D+00*EPMACH*DABS(RESULT))
  442.       RETURN
  443.       END
  444.       SUBROUTINE FUN(X,Y,A,B)
  445. C
  446. C   IN THIS SUBROUTINE THE FUNCTION 1/(P**2+1) IS EVALUATED
  447. C   WHERE P=X+iY. THIS FUNCTION IS THE LAPLACE TRANSFORM OF
  448. C   SIN(T)
  449. C
  450.       DOUBLE PRECISION X,Y,A,B,C,D
  451. C  A+IB = 1/((X+IY)**2+1)
  452.       C = X*X-Y*Y+1.0D0
  453.       D = C*C + 4.0D0*X*X*Y*Y
  454.       A = C/D
  455.       B = -2.0D0*X*Y/D
  456.       RETURN
  457.       END
  458.       DOUBLE PRECISION FUNCTION D1MACH(I)
  459. C***BEGIN PROLOGUE  D1MACH
  460. C***DATE WRITTEN   750101   (YYMMDD)
  461. C***REVISION DATE  820801   (YYMMDD)
  462. C***CATEGORY NO.  R1
  463. C***KEYWORDS  MACHINE CONSTANTS
  464. C***AUTHOR  FOX, P. A., (BELL LABS)
  465. C           HALL, A. D., (BELL LABS)
  466. C           SCHRYER, N. L., (BELL LABS)
  467. C***PURPOSE  RETURNS DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS
  468. C***DESCRIPTION
  469. C     D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
  470. C     FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
  471. C     SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
  472. C     AS FOLLOWS, FOR EXAMPLE
  473. C
  474. C          D = D1MACH(I)
  475. C
  476. C     WHERE I=1,...,5.  THE (OUTPUT) VALUE OF D ABOVE IS
  477. C     DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
  478. C     VARIOUS VALUES OF I ARE DISCUSSED BELOW.
  479. C
  480. C  DOUBLE-PRECISION MACHINE CONSTANTS
  481. C  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
  482. C  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
  483. C  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
  484. C  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
  485. C  D1MACH( 5) = LOG10(B)
  486. C***REFERENCES  FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A
  487. C                 PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL
  488. C                 SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.
  489. C***ROUTINES CALLED  XERROR
  490. C***END PROLOGUE  D1MACH
  491. C
  492.       INTEGER SMALL(4)
  493.       INTEGER LARGE(4)
  494.       INTEGER RIGHT(4)
  495.       INTEGER DIVER(4)
  496.       INTEGER LOG10(4)
  497. C
  498.       DOUBLE PRECISION DMACH(5)
  499. C
  500.       EQUIVALENCE (DMACH(1),SMALL(1))
  501.       EQUIVALENCE (DMACH(2),LARGE(1))
  502.       EQUIVALENCE (DMACH(3),RIGHT(1))
  503.       EQUIVALENCE (DMACH(4),DIVER(1))
  504.       EQUIVALENCE (DMACH(5),LOG10(1))
  505. C
  506. C     Machine Constants for Acorn Cambridge Co-Processor
  507. C    (expressed in Hexadecimal)
  508. C
  509.       DATA SMALL(1), SMALL(2) / ?I00000000, ?I00100000 /
  510.       DATA LARGE(1), LARGE(2) / ?IFFFFFFFF, ?I7FEFFFFF /
  511.       DATA RIGHT(1), RIGHT(2) / ?I00000000, ?I3CA00000 /
  512.       DATA DIVER(1), DIVER(2) / ?I00000000, ?I3CB00000 /
  513.       DATA LOG10(1), LOG10(2) / ?I509F79FF, ?I3FD34413 /
  514. C
  515. C     DMACH(1)=B**(-1022)
  516. C     DMACH(2)=B**(1023)*(B-B**(-52))
  517. C     DMACH(3)=B**(-53)
  518. C     DMACH(4)=B**(-52)
  519. C     DMACH(5)=DLOG10(B)
  520. C
  521. C     where B=2.0D00
  522. C
  523. C***FIRST EXECUTABLE STATEMENT  D1MACH
  524. C
  525.       D1MACH=DMACH(I)
  526.       RETURN
  527.       END
  528.