home *** CD-ROM | disk | FTP | other *** search
/ Best Objectech Shareware Selections / UNTITLED.iso / boss / educ / math / 011 / ders4.for < prev    next >
Encoding:
Text File  |  1993-01-03  |  20.0 KB  |  552 lines
  1. C ===  Derivating with respect to: 
  2. C     A(1,1) A(1,2) A(1,3) A(1,4) A(2,1) A(2,2) A(2,3) A(2,4) A44 B(1) B(2)
  3. C     C(1,1) C(1,2) C(1,3) C(1,4) C(2,1) C(2,2) C(2,3) C44                 
  4. C
  5.       SUBROUTINE VALUE(A,B,C,A44,C44,FN,INF)
  6.       IMPLICIT REAL*8(A-H,O-Z)
  7.       REAL*4 X
  8.       DIMENSION A(2,4),B(2),C(2,4)
  9. C ****** SUBST. MAX. NUMBER OF OBSERVATIONS ******
  10.       COMMON X(13,1800),NOBS,B3
  11.       COMMON/TRAPCODE/ITR
  12.       COMMON /VRNCS/ W1,W2,W3
  13.       ITR=0
  14. C -- COMPUTE EIGENVALUES AND EIGENVECTORS -
  15.       DET_1=A(2,2)
  16.       DET_2=-A(2,1)
  17.       DET_5=-A(1,2)
  18.       DET_6=A(1,1)
  19.       DET=A(1,1)*A(2,2)-A(1,2)*A(2,1)
  20.       HFTRA_1=0.5
  21. C*** WARNING: New identifier HFTRA_1 too long ***
  22.       HFTRA_6=0.5
  23. C*** WARNING: New identifier HFTRA_6 too long ***
  24.       HFTRA=0.5*(A(1,1)+A(2,2))
  25.       DISCR_1=2*HFTRA**(2-1)*HFTRA_1-DET_1
  26. C*** WARNING: New identifier DISCR_1 too long ***
  27.       DISCR_2=-DET_2
  28. C*** WARNING: New identifier DISCR_2 too long ***
  29.       DISCR_5=-DET_5
  30. C*** WARNING: New identifier DISCR_5 too long ***
  31.       DISCR_6=2*HFTRA**(2-1)*HFTRA_6-DET_6
  32. C*** WARNING: New identifier DISCR_6 too long ***
  33.       DISCR=HFTRA**2-DET
  34.       IF(DISCR.GT.0.) GO TO 1
  35.          WRITE(3,11)
  36.    11    FORMAT(' COMPLEX EIGENVALUES')
  37.          INF=1
  38.          RETURN
  39.     1 IF(HFTRA.lt.0.) goto 91
  40.       XL1_1=HFTRA_1+DISCR_1/2./DSQRT(DISCR)
  41.       XL1_2=DISCR_2/2./DSQRT(DISCR)
  42.       XL1_5=DISCR_5/2./DSQRT(DISCR)
  43.       XL1_6=HFTRA_6+DISCR_6/2./DSQRT(DISCR)
  44.       XL1=HFTRA+DSQRT(DISCR)
  45.  91   continue      
  46.       IF(HFTRA.ge.0.) goto 92
  47.       XL1_1=HFTRA_1-DISCR_1/2./DSQRT(DISCR)
  48.       XL1_2=-DISCR_2/2./DSQRT(DISCR)
  49.       XL1_5=-DISCR_5/2./DSQRT(DISCR)
  50.       XL1_6=HFTRA_6-DISCR_6/2./DSQRT(DISCR)
  51.       XL1=HFTRA-DSQRT(DISCR)
  52.  92   continue
  53.       XL2_1=(DET_1-DET/XL1*XL1_1)/XL1
  54.       XL2_2=(DET_2-DET/XL1*XL1_2)/XL1
  55.       XL2_5=(DET_5-DET/XL1*XL1_5)/XL1
  56.       XL2_6=(DET_6-DET/XL1*XL1_6)/XL1
  57.       XL2=DET/XL1
  58.       P11_1=XL1_1
  59.       P11_2=XL1_2
  60.       P11_5=XL1_5
  61.       P11_6=XL1_6-1.
  62.       P11=XL1-A(2,2)
  63.       P12_2=1.
  64.       P12=A(1,2)
  65.       P13_1=(A(1,3)*P11_1-(A(1,3)*P11+A(2,3)*P12)/(XL1+1.)*XL1_1)/(XL1+1
  66.      :.)
  67.       P13_2=(A(1,3)*P11_2+A(2,3)*P12_2-(A(1,3)*P11+A(2,3)*P12)/(XL1+1.)*
  68.      :XL1_2)/(XL1+1.)
  69.       P13_3=P11/(XL1+1.)
  70.       P13_5=(A(1,3)*P11_5-(A(1,3)*P11+A(2,3)*P12)/(XL1+1.)*XL1_5)/(XL1+1
  71.      :.)
  72.       P13_6=(A(1,3)*P11_6-(A(1,3)*P11+A(2,3)*P12)/(XL1+1.)*XL1_6)/(XL1+1
  73.      :.)
  74.       P13_7=P12/(XL1+1.)
  75.       P13=(A(1,3)*P11+A(2,3)*P12)/(XL1+1.)
  76.       P14_1=(A(1,4)*P11_1-(A(1,4)*P11+A(2,4)*P12)/(XL1-A44)*XL1_1)/(XL1-
  77.      :A44)
  78.       P14_2=(A(1,4)*P11_2+A(2,4)*P12_2-(A(1,4)*P11+A(2,4)*P12)/(XL1-A44)
  79.      :*XL1_2)/(XL1-A44)
  80.       P14_4=P11/(XL1-A44)
  81.       P14_5=(A(1,4)*P11_5-(A(1,4)*P11+A(2,4)*P12)/(XL1-A44)*XL1_5)/(XL1-
  82.      :A44)
  83.       P14_6=(A(1,4)*P11_6-(A(1,4)*P11+A(2,4)*P12)/(XL1-A44)*XL1_6)/(XL1-
  84.      :A44)
  85.       P14_8=P12/(XL1-A44)
  86.       P14_9=(-(A(1,4)*P11+A(2,4)*P12)/(XL1-A44)*(-1.))/(XL1-A44)
  87.       P14=(A(1,4)*P11+A(2,4)*P12)/(XL1-A44)
  88.       P21_5=1.
  89.       P21=A(2,1)
  90.       P22_1=XL2_1-1.
  91.       P22_2=XL2_2
  92.       P22_5=XL2_5
  93.       P22_6=XL2_6
  94.       P22=XL2-A(1,1)
  95.       P23_1=(A(2,3)*P22_1-(A(1,3)*P21+A(2,3)*P22)/(XL2+1.)*XL2_1)/(XL2+1
  96.      :.)
  97.       P23_2=(A(2,3)*P22_2-(A(1,3)*P21+A(2,3)*P22)/(XL2+1.)*XL2_2)/(XL2+1
  98.      :.)
  99.       P23_3=P21/(XL2+1.)
  100.       P23_5=(A(1,3)*P21_5+A(2,3)*P22_5-(A(1,3)*P21+A(2,3)*P22)/(XL2+1.)*
  101.      :XL2_5)/(XL2+1.)
  102.       P23_6=(A(2,3)*P22_6-(A(1,3)*P21+A(2,3)*P22)/(XL2+1.)*XL2_6)/(XL2+1
  103.      :.)
  104.       P23_7=P22/(XL2+1.)
  105.       P23=(A(1,3)*P21+A(2,3)*P22)/(XL2+1.)
  106.       P24_1=(A(2,4)*P22_1-(A(1,4)*P21+A(2,4)*P22)/(XL2-A44)*XL2_1)/(XL2-
  107.      :A44)
  108.       P24_2=(A(2,4)*P22_2-(A(1,4)*P21+A(2,4)*P22)/(XL2-A44)*XL2_2)/(XL2-
  109.      :A44)
  110.       P24_4=P21/(XL2-A44)
  111.       P24_5=(A(1,4)*P21_5+A(2,4)*P22_5-(A(1,4)*P21+A(2,4)*P22)/(XL2-A44)
  112.      :*XL2_5)/(XL2-A44)
  113.       P24_6=(A(2,4)*P22_6-(A(1,4)*P21+A(2,4)*P22)/(XL2-A44)*XL2_6)/(XL2-
  114.      :A44)
  115.       P24_8=P22/(XL2-A44)
  116.       P24_9=(-(A(1,4)*P21+A(2,4)*P22)/(XL2-A44)*(-1.))/(XL2-A44)
  117.       P24=(A(1,4)*P21+A(2,4)*P22)/(XL2-A44)
  118. C -- COMPUTE LOG-LIKELIHOOD -
  119.       DETC_12=C(2,2)
  120. C*** WARNING: New identifier DETC_12 too long ***
  121.       DETC_13=-C(2,1)
  122. C*** WARNING: New identifier DETC_13 too long ***
  123.       DETC_16=-C(1,2)
  124. C*** WARNING: New identifier DETC_16 too long ***
  125.       DETC_17=C(1,1)
  126. C*** WARNING: New identifier DETC_17 too long ***
  127.       DETC=C(1,1)*C(2,2)-C(1,2)*C(2,1)
  128.       IF(ABS(DETC).GT.0.) GO TO 3
  129.          WRITE(3,33)
  130.    33    FORMAT(' C IS SINGULAR')
  131.          INF=1
  132.          RETURN
  133.  *** Use CONTINUE for label ***
  134.       DETP_1=P11_1*P22+P11*P22_1
  135.       DETP_2=P11_2*P22+P11*P22_2-P12_2*P21
  136.       DETP_5=P11_5*P22+P11*P22_5-P12*P21_5
  137.       DETP_6=P11_6*P22+P11*P22_6
  138.     3 DETP=P11*P22-P12*P21
  139.       IF(ABS(DETP).GT.0.) GO TO 4
  140.          WRITE(3,44)
  141.    44    FORMAT(' P IS SINGULAR')
  142.          INF=1
  143.          RETURN
  144.     4 continue
  145.       FN_1=-DSIGN(1.,DETC*DETP)*DETC*DETP_1/DABS(DETC*DETP)
  146.       FN_2=-DSIGN(1.,DETC*DETP)*DETC*DETP_2/DABS(DETC*DETP)
  147.       FN_3=0.
  148.       FN_4=0.
  149.       FN_5=-DSIGN(1.,DETC*DETP)*DETC*DETP_5/DABS(DETC*DETP)
  150.       FN_6=-DSIGN(1.,DETC*DETP)*DETC*DETP_6/DABS(DETC*DETP)
  151.       FN_7=0.
  152.       FN_8=0.
  153.       FN_9=0.
  154.       FN_10=0.
  155.       FN_11=0.
  156.       FN_12=-DSIGN(1.,DETC*DETP)*DETC_12*DETP/DABS(DETC*DETP)
  157.       FN_13=-DSIGN(1.,DETC*DETP)*DETC_13*DETP/DABS(DETC*DETP)
  158.       FN_14=0.
  159.       FN_15=0.
  160.       FN_16=-DSIGN(1.,DETC*DETP)*DETC_16*DETP/DABS(DETC*DETP)
  161.       FN_17=-DSIGN(1.,DETC*DETP)*DETC_17*DETP/DABS(DETC*DETP)
  162.       FN_18=0.
  163.       FN_19=0.
  164. C*** WARNING: Statement interpreted as
  165. C***  FN=-DLOG(DABS(DETC*DETP))
  166.       FN=-DLOG(ABS(DETC*DETP))
  167.       W1_1=0.
  168.       W1_2=0.
  169.       W1_3=0.
  170.       W1_4=0.
  171.       W1_5=0.
  172.       W1_6=0.
  173.       W1_7=0.
  174.       W1_8=0.
  175.       W1_9=0.
  176.       W1_10=0.
  177.       W1_11=0.
  178.       W1_12=0.
  179.       W1_13=0.
  180.       W1_14=0.
  181.       W1_15=0.
  182.       W1_16=0.
  183.       W1_17=0.
  184.       W1_18=0.
  185.       W1_19=0.
  186.       W1=0.
  187.       W2_1=0.
  188.       W2_2=0.
  189.       W2_3=0.
  190.       W2_4=0.
  191.       W2_5=0.
  192.       W2_6=0.
  193.       W2_7=0.
  194.       W2_8=0.
  195.       W2_9=0.
  196.       W2_10=0.
  197.       W2_11=0.
  198.       W2_12=0.
  199.       W2_13=0.
  200.       W2_14=0.
  201.       W2_15=0.
  202.       W2_16=0.
  203.       W2_17=0.
  204.       W2_18=0.
  205.       W2_19=0.
  206.       W2=0.
  207.       W3=0.
  208.       PB1_1=P11_1*B(1)+P13_1*B3-P14_1*A44
  209.       PB1_2=P11_2*B(1)+P12_2*B(2)+P13_2*B3-P14_2*A44
  210.       PB1_3=P13_3*B3
  211.       PB1_4=-P14_4*A44
  212.       PB1_5=P11_5*B(1)+P13_5*B3-P14_5*A44
  213.       PB1_6=P11_6*B(1)+P13_6*B3-P14_6*A44
  214.       PB1_7=P13_7*B3
  215.       PB1_8=-P14_8*A44
  216.       PB1_9=-P14_9*A44-P14
  217.       PB1_10=P11
  218.       PB1_11=P12
  219.       PB1=P11*B(1)+P12*B(2)+P13*B3-P14*A44
  220.       PB2_1=P22_1*B(2)+P23_1*B3-P24_1*A44
  221.       PB2_2=P22_2*B(2)+P23_2*B3-P24_2*A44
  222.       PB2_3=P23_3*B3
  223.       PB2_4=-P24_4*A44
  224.       PB2_5=P21_5*B(1)+P22_5*B(2)+P23_5*B3-P24_5*A44
  225.       PB2_6=P22_6*B(2)+P23_6*B3-P24_6*A44
  226.       PB2_7=P23_7*B3
  227.       PB2_8=-P24_8*A44
  228.       PB2_9=-P24_9*A44-P24
  229.       PB2_10=P21
  230.       PB2_11=P22
  231.       PB2=P21*B(1)+P22*B(2)+P23*B3-P24*A44
  232. C -- LOOP BEGINS -
  233.       DO 10 I=1,NOBS
  234.       IF(ITR.EQ.1) GO TO 2
  235.       Y04_9=-(1.0-DEXP(X(11,I)*C44))*DEXP(X(12,I)*(-A44))*X(12,I)*(-1.)
  236.       Y04_19=-(-DEXP(X(11,I)*C44)*X(11,I))*DEXP(X(12,I)*(-A44))
  237.       Y04=1.0-(1.0-DEXP(X(11,I)*C44))*DEXP(X(12,I)*(-A44))
  238.       Y4_9=-(1.0-DEXP(X(11,I)*C44))*DEXP(X(13,I)*(-A44))*X(13,I)*(-1.)
  239.       Y4_19=-(-DEXP(X(11,I)*C44)*X(11,I))*DEXP(X(13,I)*(-A44))
  240.       Y4=1.0-(1.0-DEXP(X(11,I)*C44))*DEXP(X(13,I)*(-A44))
  241.       ALX04_9=Y04_9/Y04/C44
  242. C*** WARNING: New identifier ALX04_9 too long ***
  243.       ALX04_19=(Y04_19/Y04-DLOG(Y04)/C44)/C44
  244. C*** WARNING: New identifier ALX04_19 too long ***
  245.       ALX04=DLOG(Y04)/C44
  246.       ALX4_9=Y4_9/Y4/C44
  247.       ALX4_19=(Y4_19/Y4-DLOG(Y4)/C44)/C44
  248. C*** WARNING: New identifier ALX4_19 too long ***
  249.       ALX4=DLOG(Y4)/C44
  250.       ALXC1_9=C(1,4)*ALX4_9
  251. C*** WARNING: New identifier ALXC1_9 too long ***
  252.       ALXC1_12=X(1,I)
  253. C*** WARNING: New identifier ALXC1_12 too long ***
  254.       ALXC1_13=X(2,I)
  255. C*** WARNING: New identifier ALXC1_13 too long ***
  256.       ALXC1_14=X(3,I)
  257. C*** WARNING: New identifier ALXC1_14 too long ***
  258.       ALXC1_15=ALX4
  259. C*** WARNING: New identifier ALXC1_15 too long ***
  260.       ALXC1_19=C(1,4)*ALX4_19
  261. C*** WARNING: New identifier ALXC1_19 too long ***
  262.       ALXC1=C(1,1)*X(1,I)+C(1,2)*X(2,I)+C(1,3)*X(3,I)+C(1,4)*ALX4
  263.       ALXC2_9=C(2,4)*ALX4_9
  264. C*** WARNING: New identifier ALXC2_9 too long ***
  265.       ALXC2_16=X(1,I)
  266. C*** WARNING: New identifier ALXC2_16 too long ***
  267.       ALXC2_17=X(2,I)
  268. C*** WARNING: New identifier ALXC2_17 too long ***
  269.       ALXC2_18=X(3,I)
  270. C*** WARNING: New identifier ALXC2_18 too long ***
  271.       ALXC2_19=C(2,4)*ALX4_19
  272. C*** WARNING: New identifier ALXC2_19 too long ***
  273.       ALXC2=C(2,1)*X(1,I)+C(2,2)*X(2,I)+C(2,3)*X(3,I)+C(2,4)*ALX4
  274.       Y1_9=DEXP(ALXC1)*ALXC1_9
  275.       Y1_12=DEXP(ALXC1)*ALXC1_12
  276.       Y1_13=DEXP(ALXC1)*ALXC1_13
  277.       Y1_14=DEXP(ALXC1)*ALXC1_14
  278.       Y1_15=DEXP(ALXC1)*ALXC1_15
  279.       Y1_19=DEXP(ALXC1)*ALXC1_19
  280.       Y1=DEXP(ALXC1)
  281.       Y2_9=DEXP(ALXC2)*ALXC2_9
  282.       Y2_16=DEXP(ALXC2)*ALXC2_16
  283.       Y2_17=DEXP(ALXC2)*ALXC2_17
  284.       Y2_18=DEXP(ALXC2)*ALXC2_18
  285.       Y2_19=DEXP(ALXC2)*ALXC2_19
  286.       Y2=DEXP(ALXC2)
  287.       Y01_9=DEXP(C(1,1)*X(4,I)+C(1,2)*X(5,I)+C(1,3)*X(6,I)+C(1,4)*ALX04)
  288.      :*C(1,4)*ALX04_9
  289.       Y01_12=DEXP(C(1,1)*X(4,I)+C(1,2)*X(5,I)+C(1,3)*X(6,I)+C(1,4)*ALX04
  290.      :)*X(4,I)
  291.       Y01_13=DEXP(C(1,1)*X(4,I)+C(1,2)*X(5,I)+C(1,3)*X(6,I)+C(1,4)*ALX04
  292.      :)*X(5,I)
  293.       Y01_14=DEXP(C(1,1)*X(4,I)+C(1,2)*X(5,I)+C(1,3)*X(6,I)+C(1,4)*ALX04
  294.      :)*X(6,I)
  295.       Y01_15=DEXP(C(1,1)*X(4,I)+C(1,2)*X(5,I)+C(1,3)*X(6,I)+C(1,4)*ALX04
  296.      :)*ALX04
  297.       Y01_19=DEXP(C(1,1)*X(4,I)+C(1,2)*X(5,I)+C(1,3)*X(6,I)+C(1,4)*ALX04
  298.      :)*C(1,4)*ALX04_19
  299.       Y01=DEXP(C(1,1)*X(4,I)+C(1,2)*X(5,I)+C(1,3)*X(6,I)+C(1,4)*ALX04)
  300.       Y02_9=DEXP(C(2,1)*X(4,I)+C(2,2)*X(5,I)+C(2,3)*X(6,I)+C(2,4)*ALX04)
  301.      :*C(2,4)*ALX04_9
  302.       Y02_16=DEXP(C(2,1)*X(4,I)+C(2,2)*X(5,I)+C(2,3)*X(6,I)+C(2,4)*ALX04
  303.      :)*X(4,I)
  304.       Y02_17=DEXP(C(2,1)*X(4,I)+C(2,2)*X(5,I)+C(2,3)*X(6,I)+C(2,4)*ALX04
  305.      :)*X(5,I)
  306.       Y02_18=DEXP(C(2,1)*X(4,I)+C(2,2)*X(5,I)+C(2,3)*X(6,I)+C(2,4)*ALX04
  307.      :)*X(6,I)
  308.       Y02_19=DEXP(C(2,1)*X(4,I)+C(2,2)*X(5,I)+C(2,3)*X(6,I)+C(2,4)*ALX04
  309.      :)*C(2,4)*ALX04_19
  310.       Y02=DEXP(C(2,1)*X(4,I)+C(2,2)*X(5,I)+C(2,3)*X(6,I)+C(2,4)*ALX04)
  311.       EL1T_1=DEXP(X(7,I)*(-XL1))*X(7,I)*(-XL1_1)
  312.       EL1T_2=DEXP(X(7,I)*(-XL1))*X(7,I)*(-XL1_2)
  313.       EL1T_5=DEXP(X(7,I)*(-XL1))*X(7,I)*(-XL1_5)
  314.       EL1T_6=DEXP(X(7,I)*(-XL1))*X(7,I)*(-XL1_6)
  315.       EL1T=DEXP(X(7,I)*(-XL1))
  316.       EL2T_1=DEXP(X(7,I)*(-XL2))*X(7,I)*(-XL2_1)
  317.       EL2T_2=DEXP(X(7,I)*(-XL2))*X(7,I)*(-XL2_2)
  318.       EL2T_5=DEXP(X(7,I)*(-XL2))*X(7,I)*(-XL2_5)
  319.       EL2T_6=DEXP(X(7,I)*(-XL2))*X(7,I)*(-XL2_6)
  320.       EL2T=DEXP(X(7,I)*(-XL2))
  321.       E1_1=P11_1*Y1+P13_1*X(8,I)+P14_1*Y4-EL1T_1*(P11*Y01+P12*Y02+P13*X(
  322.      :9,I)+P14*Y04)-EL1T*(P11_1*Y01+P13_1*X(9,I)+P14_1*Y04)
  323.       E1_2=P11_2*Y1+P12_2*Y2+P13_2*X(8,I)+P14_2*Y4-EL1T_2*(P11*Y01+P12*Y
  324.      :02+P13*X(9,I)+P14*Y04)-EL1T*(P11_2*Y01+P12_2*Y02+P13_2*X(9,I)+P14_
  325.      :2*Y04)
  326.       E1_3=P13_3*X(8,I)-EL1T*P13_3*X(9,I)
  327.       E1_4=P14_4*Y4-EL1T*P14_4*Y04
  328.       E1_5=P11_5*Y1+P13_5*X(8,I)+P14_5*Y4-EL1T_5*(P11*Y01+P12*Y02+P13*X(
  329.      :9,I)+P14*Y04)-EL1T*(P11_5*Y01+P13_5*X(9,I)+P14_5*Y04)
  330.       E1_6=P11_6*Y1+P13_6*X(8,I)+P14_6*Y4-EL1T_6*(P11*Y01+P12*Y02+P13*X(
  331.      :9,I)+P14*Y04)-EL1T*(P11_6*Y01+P13_6*X(9,I)+P14_6*Y04)
  332.       E1_7=P13_7*X(8,I)-EL1T*P13_7*X(9,I)
  333.       E1_8=P14_8*Y4-EL1T*P14_8*Y04
  334.       E1_9=P11*Y1_9+P12*Y2_9+P14_9*Y4+P14*Y4_9-EL1T*(P11*Y01_9+P12*Y02_9
  335.      :+P14_9*Y04+P14*Y04_9)
  336.       E1_10=0.
  337.       E1_11=0.
  338.       E1_12=P11*Y1_12-EL1T*P11*Y01_12
  339.       E1_13=P11*Y1_13-EL1T*P11*Y01_13
  340.       E1_14=P11*Y1_14-EL1T*P11*Y01_14
  341.       E1_15=P11*Y1_15-EL1T*P11*Y01_15
  342.       E1_16=P12*Y2_16-EL1T*P12*Y02_16
  343.       E1_17=P12*Y2_17-EL1T*P12*Y02_17
  344.       E1_18=P12*Y2_18-EL1T*P12*Y02_18
  345.       E1_19=P11*Y1_19+P12*Y2_19+P14*Y4_19-EL1T*(P11*Y01_19+P12*Y02_19+P1
  346.      :4*Y04_19)
  347.       E1=P11*Y1+P12*Y2+P13*X(8,I)+P14*Y4
  348.      1     -EL1T*(P11*Y01+P12*Y02+P13*X(9,I)+P14*Y04)
  349.       E2_1=P22_1*Y2+P23_1*X(8,I)+P24_1*Y4-EL2T_1*(P21*Y01+P22*Y02+P23*X(
  350.      :9,I)+P24*Y04)-EL2T*(P22_1*Y02+P23_1*X(9,I)+P24_1*Y04)
  351.       E2_2=P22_2*Y2+P23_2*X(8,I)+P24_2*Y4-EL2T_2*(P21*Y01+P22*Y02+P23*X(
  352.      :9,I)+P24*Y04)-EL2T*(P22_2*Y02+P23_2*X(9,I)+P24_2*Y04)
  353.       E2_3=P23_3*X(8,I)-EL2T*P23_3*X(9,I)
  354.       E2_4=P24_4*Y4-EL2T*P24_4*Y04
  355.       E2_5=P21_5*Y1+P22_5*Y2+P23_5*X(8,I)+P24_5*Y4-EL2T_5*(P21*Y01+P22*Y
  356.      :02+P23*X(9,I)+P24*Y04)-EL2T*(P21_5*Y01+P22_5*Y02+P23_5*X(9,I)+P24_
  357.      :5*Y04)
  358.       E2_6=P22_6*Y2+P23_6*X(8,I)+P24_6*Y4-EL2T_6*(P21*Y01+P22*Y02+P23*X(
  359.      :9,I)+P24*Y04)-EL2T*(P22_6*Y02+P23_6*X(9,I)+P24_6*Y04)
  360.       E2_7=P23_7*X(8,I)-EL2T*P23_7*X(9,I)
  361.       E2_8=P24_8*Y4-EL2T*P24_8*Y04
  362.       E2_9=P21*Y1_9+P22*Y2_9+P24_9*Y4+P24*Y4_9-EL2T*(P21*Y01_9+P22*Y02_9
  363.      :+P24_9*Y04+P24*Y04_9)
  364.       E2_10=0.
  365.       E2_11=0.
  366.       E2_12=P21*Y1_12-EL2T*P21*Y01_12
  367.       E2_13=P21*Y1_13-EL2T*P21*Y01_13
  368.       E2_14=P21*Y1_14-EL2T*P21*Y01_14
  369.       E2_15=P21*Y1_15-EL2T*P21*Y01_15
  370.       E2_16=P22*Y2_16-EL2T*P22*Y02_16
  371.       E2_17=P22*Y2_17-EL2T*P22*Y02_17
  372.       E2_18=P22*Y2_18-EL2T*P22*Y02_18
  373.       E2_19=P21*Y1_19+P22*Y2_19+P24*Y4_19-EL2T*(P21*Y01_19+P22*Y02_19+P2
  374.      :4*Y04_19)
  375.       E2=P21*Y1+P22*Y2+P23*X(8,I)+P24*Y4
  376.      1     -EL2T*(P21*Y01+P22*Y02+P23*X(9,I)+P24*Y04)
  377.       IF(ABS(XL1).GT.1.D-22) GO TO 222
  378.       R1_1=0.
  379.       R1_2=0.
  380.       R1_5=0.
  381.       R1_6=0.
  382.          R1=-X(7,I)
  383.       E1_1=E1_1-R1_1*PB1-R1*PB1_1
  384.       E1_2=E1_2-R1_2*PB1-R1*PB1_2
  385.       E1_3=E1_3-R1*PB1_3
  386.       E1_4=E1_4-R1*PB1_4
  387.       E1_5=E1_5-R1_5*PB1-R1*PB1_5
  388.       E1_6=E1_6-R1_6*PB1-R1*PB1_6
  389.       E1_7=E1_7-R1*PB1_7
  390.       E1_8=E1_8-R1*PB1_8
  391.       E1_9=E1_9-R1*PB1_9
  392.       E1_10=E1_10-R1*PB1_10
  393.       E1_11=E1_11-R1*PB1_11
  394.          E1=E1-R1*PB1
  395.          GO TO 333
  396.  *** Use CONTINUE for label ***
  397.       R1_1=(EL1T_1*EL1T+EL1T*EL1T_1-(EL1T*EL1T-1.)/(XL1+XL1)*(XL1_1+XL1_
  398.      :1))/(XL1+XL1)
  399.       R1_2=(EL1T_2*EL1T+EL1T*EL1T_2-(EL1T*EL1T-1.)/(XL1+XL1)*(XL1_2+XL1_
  400.      :2))/(XL1+XL1)
  401.       R1_5=(EL1T_5*EL1T+EL1T*EL1T_5-(EL1T*EL1T-1.)/(XL1+XL1)*(XL1_5+XL1_
  402.      :5))/(XL1+XL1)
  403.       R1_6=(EL1T_6*EL1T+EL1T*EL1T_6-(EL1T*EL1T-1.)/(XL1+XL1)*(XL1_6+XL1_
  404.      :6))/(XL1+XL1)
  405.  222  R1=(EL1T*EL1T-1.)/(XL1+XL1)
  406.       E1_1=E1_1+((-EL1T_1)-(1.-EL1T)/XL1*XL1_1)/XL1*PB1+(1.-EL1T)/XL1*PB
  407.      :1_1
  408.       E1_2=E1_2+((-EL1T_2)-(1.-EL1T)/XL1*XL1_2)/XL1*PB1+(1.-EL1T)/XL1*PB
  409.      :1_2
  410.       E1_3=E1_3+(1.-EL1T)/XL1*PB1_3
  411.       E1_4=E1_4+(1.-EL1T)/XL1*PB1_4
  412.       E1_5=E1_5+((-EL1T_5)-(1.-EL1T)/XL1*XL1_5)/XL1*PB1+(1.-EL1T)/XL1*PB
  413.      :1_5
  414.       E1_6=E1_6+((-EL1T_6)-(1.-EL1T)/XL1*XL1_6)/XL1*PB1+(1.-EL1T)/XL1*PB
  415.      :1_6
  416.       E1_7=E1_7+(1.-EL1T)/XL1*PB1_7
  417.       E1_8=E1_8+(1.-EL1T)/XL1*PB1_8
  418.       E1_9=E1_9+(1.-EL1T)/XL1*PB1_9
  419.       E1_10=E1_10+(1.-EL1T)/XL1*PB1_10
  420.       E1_11=E1_11+(1.-EL1T)/XL1*PB1_11
  421.       E1=E1+(1.-EL1T)/XL1*PB1
  422.  333  IF(ABS(XL2).GT.1.D-22) GO TO 444
  423.       R2_1=0.
  424.       R2_2=0.
  425.       R2_5=0.
  426.       R2_6=0.
  427.          R2=-X(7,I)
  428.       E2_1=E2_1-R2_1*PB2-R2*PB2_1
  429.       E2_2=E2_2-R2_2*PB2-R2*PB2_2
  430.       E2_3=E2_3-R2*PB2_3
  431.       E2_4=E2_4-R2*PB2_4
  432.       E2_5=E2_5-R2_5*PB2-R2*PB2_5
  433.       E2_6=E2_6-R2_6*PB2-R2*PB2_6
  434.       E2_7=E2_7-R2*PB2_7
  435.       E2_8=E2_8-R2*PB2_8
  436.       E2_9=E2_9-R2*PB2_9
  437.       E2_10=E2_10-R2*PB2_10
  438.       E2_11=E2_11-R2*PB2_11
  439.          E2=E2-R2*PB2
  440.          GO TO 555
  441.  *** Use CONTINUE for label ***
  442.       R2_1=(EL2T_1*EL2T+EL2T*EL2T_1-(EL2T*EL2T-1.)/(XL2+XL2)*(XL2_1+XL2_
  443.      :1))/(XL2+XL2)
  444.       R2_2=(EL2T_2*EL2T+EL2T*EL2T_2-(EL2T*EL2T-1.)/(XL2+XL2)*(XL2_2+XL2_
  445.      :2))/(XL2+XL2)
  446.       R2_5=(EL2T_5*EL2T+EL2T*EL2T_5-(EL2T*EL2T-1.)/(XL2+XL2)*(XL2_5+XL2_
  447.      :5))/(XL2+XL2)
  448.       R2_6=(EL2T_6*EL2T+EL2T*EL2T_6-(EL2T*EL2T-1.)/(XL2+XL2)*(XL2_6+XL2_
  449.      :6))/(XL2+XL2)
  450.  444  R2=(EL2T*EL2T-1.)/(XL2+XL2)
  451.       E2_1=E2_1+((-EL2T_1)-(1.-EL2T)/XL2*XL2_1)/XL2*PB2+(1.-EL2T)/XL2*PB
  452.      :2_1
  453.       E2_2=E2_2+((-EL2T_2)-(1.-EL2T)/XL2*XL2_2)/XL2*PB2+(1.-EL2T)/XL2*PB
  454.      :2_2
  455.       E2_3=E2_3+(1.-EL2T)/XL2*PB2_3
  456.       E2_4=E2_4+(1.-EL2T)/XL2*PB2_4
  457.       E2_5=E2_5+((-EL2T_5)-(1.-EL2T)/XL2*XL2_5)/XL2*PB2+(1.-EL2T)/XL2*PB
  458.      :2_5
  459.       E2_6=E2_6+((-EL2T_6)-(1.-EL2T)/XL2*XL2_6)/XL2*PB2+(1.-EL2T)/XL2*PB
  460.      :2_6
  461.       E2_7=E2_7+(1.-EL2T)/XL2*PB2_7
  462.       E2_8=E2_8+(1.-EL2T)/XL2*PB2_8
  463.       E2_9=E2_9+(1.-EL2T)/XL2*PB2_9
  464.       E2_10=E2_10+(1.-EL2T)/XL2*PB2_10
  465.       E2_11=E2_11+(1.-EL2T)/XL2*PB2_11
  466.       E2=E2+(1.-EL2T)/XL2*PB2
  467.  *** Use CONTINUE for label ***
  468.  555  R3=(1.-DEXP(2.D0*X(7,I)))*0.5
  469.       W1_1=W1_1+(E1_1*E1+E1*E1_1-E1*E1/R1*R1_1)/R1
  470.       W1_2=W1_2+(E1_2*E1+E1*E1_2-E1*E1/R1*R1_2)/R1
  471.       W1_3=W1_3+(E1_3*E1+E1*E1_3)/R1
  472.       W1_4=W1_4+(E1_4*E1+E1*E1_4)/R1
  473.       W1_5=W1_5+(E1_5*E1+E1*E1_5-E1*E1/R1*R1_5)/R1
  474.       W1_6=W1_6+(E1_6*E1+E1*E1_6-E1*E1/R1*R1_6)/R1
  475.       W1_7=W1_7+(E1_7*E1+E1*E1_7)/R1
  476.       W1_8=W1_8+(E1_8*E1+E1*E1_8)/R1
  477.       W1_9=W1_9+(E1_9*E1+E1*E1_9)/R1
  478.       W1_10=W1_10+(E1_10*E1+E1*E1_10)/R1
  479.       W1_11=W1_11+(E1_11*E1+E1*E1_11)/R1
  480.       W1_12=W1_12+(E1_12*E1+E1*E1_12)/R1
  481.       W1_13=W1_13+(E1_13*E1+E1*E1_13)/R1
  482.       W1_14=W1_14+(E1_14*E1+E1*E1_14)/R1
  483.       W1_15=W1_15+(E1_15*E1+E1*E1_15)/R1
  484.       W1_16=W1_16+(E1_16*E1+E1*E1_16)/R1
  485.       W1_17=W1_17+(E1_17*E1+E1*E1_17)/R1
  486.       W1_18=W1_18+(E1_18*E1+E1*E1_18)/R1
  487.       W1_19=W1_19+(E1_19*E1+E1*E1_19)/R1
  488.       W1=W1+E1*E1/R1
  489.       W2_1=W2_1+(E2_1*E2+E2*E2_1-E2*E2/R2*R2_1)/R2
  490.       W2_2=W2_2+(E2_2*E2+E2*E2_2-E2*E2/R2*R2_2)/R2
  491.       W2_3=W2_3+(E2_3*E2+E2*E2_3)/R2
  492.       W2_4=W2_4+(E2_4*E2+E2*E2_4)/R2
  493.       W2_5=W2_5+(E2_5*E2+E2*E2_5-E2*E2/R2*R2_5)/R2
  494.       W2_6=W2_6+(E2_6*E2+E2*E2_6-E2*E2/R2*R2_6)/R2
  495.       W2_7=W2_7+(E2_7*E2+E2*E2_7)/R2
  496.       W2_8=W2_8+(E2_8*E2+E2*E2_8)/R2
  497.       W2_9=W2_9+(E2_9*E2+E2*E2_9)/R2
  498.       W2_10=W2_10+(E2_10*E2+E2*E2_10)/R2
  499.       W2_11=W2_11+(E2_11*E2+E2*E2_11)/R2
  500.       W2_12=W2_12+(E2_12*E2+E2*E2_12)/R2
  501.       W2_13=W2_13+(E2_13*E2+E2*E2_13)/R2
  502.       W2_14=W2_14+(E2_14*E2+E2*E2_14)/R2
  503.       W2_15=W2_15+(E2_15*E2+E2*E2_15)/R2
  504.       W2_16=W2_16+(E2_16*E2+E2*E2_16)/R2
  505.       W2_17=W2_17+(E2_17*E2+E2*E2_17)/R2
  506.       W2_18=W2_18+(E2_18*E2+E2*E2_18)/R2
  507.       W2_19=W2_19+(E2_19*E2+E2*E2_19)/R2
  508.       W2=W2+E2*E2/R2
  509.       W3=W3+X(10,I)*X(10,I)/R3
  510.       FN_1=FN_1+0.5*(R1_1*R2+R1*R2_1)*R3/R1/R2/R3/NOBS
  511.       FN_2=FN_2+0.5*(R1_2*R2+R1*R2_2)*R3/R1/R2/R3/NOBS
  512.       FN_5=FN_5+0.5*(R1_5*R2+R1*R2_5)*R3/R1/R2/R3/NOBS
  513.       FN_6=FN_6+0.5*(R1_6*R2+R1*R2_6)*R3/R1/R2/R3/NOBS
  514.       FN_9=FN_9+((-ALXC1_9)-ALXC2_9)/NOBS
  515.       FN_12=FN_12+(-ALXC1_12)/NOBS
  516.       FN_13=FN_13+(-ALXC1_13)/NOBS
  517.       FN_14=FN_14+(-ALXC1_14)/NOBS
  518.       FN_15=FN_15+(-ALXC1_15)/NOBS
  519.       FN_16=FN_16+(-ALXC2_16)/NOBS
  520.       FN_17=FN_17+(-ALXC2_17)/NOBS
  521.       FN_18=FN_18+(-ALXC2_18)/NOBS
  522.       FN_19=FN_19+((-ALXC1_19)-ALXC2_19)/NOBS
  523.       FN=FN+(0.5*DLOG(R1*R2*R3)-ALXC1-ALXC2)/NOBS
  524.    10 continue
  525. C -- LOOP ENDS -
  526.       FN_1=FN_1+0.5*(W1_1*W2+W1*W2_1)*W3/W1/W2/W3
  527.       FN_2=FN_2+0.5*(W1_2*W2+W1*W2_2)*W3/W1/W2/W3
  528.       FN_3=FN_3+0.5*(W1_3*W2+W1*W2_3)*W3/W1/W2/W3
  529.       FN_4=FN_4+0.5*(W1_4*W2+W1*W2_4)*W3/W1/W2/W3
  530.       FN_5=FN_5+0.5*(W1_5*W2+W1*W2_5)*W3/W1/W2/W3
  531.       FN_6=FN_6+0.5*(W1_6*W2+W1*W2_6)*W3/W1/W2/W3
  532.       FN_7=FN_7+0.5*(W1_7*W2+W1*W2_7)*W3/W1/W2/W3
  533.       FN_8=FN_8+0.5*(W1_8*W2+W1*W2_8)*W3/W1/W2/W3
  534.       FN_9=FN_9+0.5*(W1_9*W2+W1*W2_9)*W3/W1/W2/W3
  535.       FN_10=FN_10+0.5*(W1_10*W2+W1*W2_10)*W3/W1/W2/W3
  536.       FN_11=FN_11+0.5*(W1_11*W2+W1*W2_11)*W3/W1/W2/W3
  537.       FN_12=FN_12+0.5*(W1_12*W2+W1*W2_12)*W3/W1/W2/W3
  538.       FN_13=FN_13+0.5*(W1_13*W2+W1*W2_13)*W3/W1/W2/W3
  539.       FN_14=FN_14+0.5*(W1_14*W2+W1*W2_14)*W3/W1/W2/W3
  540.       FN_15=FN_15+0.5*(W1_15*W2+W1*W2_15)*W3/W1/W2/W3
  541.       FN_16=FN_16+0.5*(W1_16*W2+W1*W2_16)*W3/W1/W2/W3
  542.       FN_17=FN_17+0.5*(W1_17*W2+W1*W2_17)*W3/W1/W2/W3
  543.       FN_18=FN_18+0.5*(W1_18*W2+W1*W2_18)*W3/W1/W2/W3
  544.       FN_19=FN_19+0.5*(W1_19*W2+W1*W2_19)*W3/W1/W2/W3
  545. C*** WARNING: Statement interpreted as
  546. C***  FN=FN+0.5*DLOG(W1*W2*W3)
  547.       FN=FN+0.5*(DLOG(W1*W2*W3))
  548.     2 CONTINUE
  549.       INF=ITR
  550.       RETURN
  551.       END
  552.