C
C     PROGRAM NAME :  EWMA
C     Calling Program to Evaluate the Average Run Lengths
C     of an EWMA with parameters Lambda and K, for shifts
C     in the process mean ranging from 0 to 4 sigma
C
C     Stephen V. Crowder
C     Sandia National Laboratories
C
C     February 13, 1995
C
       REAL*8 AL,ARL,ARG,L,A(24,24),B(24),W(24),P(24),
     1 X(24),F,D,DINC,H,WK(24)
       INTEGER IPIVOT(24),IFLAG,ND
       WRITE(*,2)
  2    FORMAT(/)
       WRITE(*,3)
  3    FORMAT(10X,'Evaluate Average Run Lengths of an EWMA ',
     1 'With Parameters Lambda and K',//,10X,'Enter Lambda and K: '
     1 ,/)
       READ(*,*) AL,L
       H=DSQRT(AL/(2.0D0-AL))*L
       DINC=0.25D0
       ND=17
       P(1)=.9951872199970213D0
       P(2)=.9747285559713095D0
       P(3)=.9382745520027327D0
       P(4)=.8864155270044010D0
       P(5)=.8200019859739029D0
       P(6)=.7401241915785543D0
       P(7)=.6480936519369755D0
       P(8)=.5454214713888395D0
       P(9)=.4337935076260451D0
       P(10)=.3150426796961634D0
       P(11)=.1911188674736163D0
       P(12)=.0640568928626056D0
       W(1)=.0123412297999872D0
       W(2)=.0285313886289337D0
       W(3)=.0442774388174198D0
       W(4)=.0592985849154368D0
       W(5)=.0733464814110803D0
       W(6)=.0861901615319533D0
       W(7)=.0976186521041139D0
       W(8)=.1074442701159656D0
       W(9)=.1155056680537256D0
       W(10)=.1216704729278034D0
       W(11)=.1258374563468283D0
       W(12)=.1279381953467521D0
       DO 1 I=1,12
        P(25-I)=-P(I)
        W(25-I)=W(I)
  1    CONTINUE
       DO 22 I=1,24
        W(I)=H*W(I)
        P(I)=H*P(I)
        B(I)=-1.0D0
  22    CONTINUE
       DO 100 K=1,ND
        D=DINC*DBLE(K-1)
       DO 10 I=1,24
        DO 20 J=1,24
         ARG=(P(J)-(1.0D0-AL)*P(I))/AL
         IF (I .EQ. J) THEN
         A(I,J)=(1.0D0/AL)*W(I)*F(ARG-D) - 1.0D0
         ELSE
         A(I,J)=(1.0D0/AL)*W(J)*F(ARG-D)
         END IF
 20     CONTINUE
 10    CONTINUE
       CALL FACTOR(A,24,WK,IPIVOT,IFLAG)
       IF (IFLAG .EQ. 0) THEN
        WRITE(*,200)
        STOP
       END IF
       CALL SUBST(A,IPIVOT,B,24,X)
       ARL=0.0D0
       DO 30 I=1,24
        ARG=P(I)/AL
        ARL=ARL + W(I)*X(I)*F(ARG-D)
 30    CONTINUE
       ARL=1.0D0 + ARL/AL
       WRITE(*,60) AL,L,D,ARL
 60    FORMAT(10X,'Lambda= ',F4.2,2X,'K= ',F4.2,2X,'Shift= ',
     1 F4.2,2X,'ARL= ',F10.3)
200    FORMAT(10X,'ZERO DETERMINANT FOR LINEAR SYSTEM')
100    CONTINUE
       STOP
       END
       DOUBLE PRECISION FUNCTION F(X)
       DOUBLE PRECISION X
       F=3.989422804014327D-1*DEXP(-.50D0*X*X)
       RETURN
       END
       SUBROUTINE SUBST(W1,IPIVOT,B,N,X2)
       INTEGER IPIVOT(24),I,IP,J
       REAL*8 B(24),W1(24,24),X2(24),SUM
       IF (N .LE. 1) THEN
        X2(1)=B(1)/W1(1,1)
        RETURN
       END IF
       IP=IPIVOT(1)
       X2(1)=B(IP)
       DO 15 I=2,N
        SUM=0.0D0
        I1=I-1
        DO 14 J=1,I1
         SUM=W1(I,J)*X2(J) + SUM
  14    CONTINUE
        IP=IPIVOT(I)
        X2(I)=B(IP) - SUM
  15   CONTINUE
       X2(N)=X2(N)/W1(N,N)
       I2=N-1
       DO 20 II=1,I2
        I=N-II
        SUM=0.0D0
        I3=I+1
        DO 19 J=I3,N
         SUM=W1(I,J)*X2(J) + SUM
  19    CONTINUE
        X2(I)=(X2(I)-SUM)/W1(I,I)
  20   CONTINUE
       RETURN
       END
       SUBROUTINE FACTOR(W1,N,D1,IPIVOT,IFLAG)
       INTEGER IFLAG,IPIVOT(24),I,ISTAR,J,K
       REAL*8 D1(24),W1(24,24),AWIKOD,COLMAX,RATIO,ROWMAX,TEMP
       IFLAG=1
       DO 9 I=1,N
        IPIVOT(I)=I
        ROWMAX=0.0D0
        DO 5 J=1,N
         ROWMAX=DMAX1(ROWMAX,DABS(W1(I,J)))
    5   CONTINUE
        IF (ROWMAX .EQ. 0.0D0) THEN
         IFLAG=0
         ROWMAX=1.0D0
        END IF
        D1(I)=ROWMAX
    9  CONTINUE
       IF (N .LE. 1) RETURN
       N1=N-1
       DO 20 K=1,N1
        COLMAX=DABS(W1(K,K))/D1(K)
        ISTAR=K
        K1=K+1
        DO 13 I=K1,N
         AWIKOD=DABS(W1(I,K))/D1(K)
         IF (AWIKOD .GT. COLMAX) THEN
          COLMAX=AWIKOD
          ISTAR=I
         END IF
   13   CONTINUE
        IF (COLMAX .EQ. 0.0D0) THEN
         IFLAG=0
        ELSE
         IF (ISTAR .GT. K) THEN
          IFLAG=-IFLAG
          I=IPIVOT(ISTAR)
          IPIVOT(ISTAR)=IPIVOT(K)
          IPIVOT(K)=I
          TEMP=D1(ISTAR)
          D1(ISTAR)=D1(K)
          D1(K)=TEMP
          DO 15 J=1,N
           TEMP=W1(ISTAR,J)
           W1(ISTAR,J)=W1(K,J)
           W1(K,J)=TEMP
   15     CONTINUE
         END IF
         K2=K+1
         DO 19 I=K2,N
          W1(I,K)=W1(I,K)/W1(K,K)
          RATIO=W1(I,K)
          K3=K+1
          DO 18 J=K3,N
           W1(I,J)=W1(I,J)-RATIO*W1(K,J)
   18     CONTINUE
   19    CONTINUE
         END IF
   20   CONTINUE
        IF (W1(N,N) .EQ. 0.0D0) IFLAG=0
        RETURN
        END