****************************************************
*
*    Program Name : XMR
*    To compute the average run lengths for combined
*    individual and moving range charts .
*
****************************************************
*
       DOUBLE PRECISION A(161,161),B(161),X(161),M,R,D,
     1 D1(161),ARG,P,AI,BI,H,U,ARL
       INTEGER IPIVOT(161),N1,N2,IFLAG
*
*
*    User inputs:
*
*     H=subinterval width for trapezoidal
*       rule (.05 adequate) .
*
*     M=control limit for chart for individuals,
*       expressed in units of sigma from nominal.
*       (choose M which can be expressed as k*H,
*       where k is a positive integer)
*
*     R=control limit for moving range chart,
*       expressed in units of sigma. (choose R
*       which can be expressed as j*H, where j
*       is a positive integer)
*
*     D=process mean shift from nominal, expressed
*       in units of sigma.
*
*     N=2*k
*
*     NN=2*k+1  (NN is the necessary dimension for
*        the matrix A and the vectors B, X, and the
*        work vectors D1 and IPIVOT)
*
*
       WRITE(*,1)
   1   FORMAT(10X,'Compute ARLs for Individual-Moving Range Chart')
       WRITE(*,2)
   2   FORMAT(/,10X,'Enter M and R : ',/)
       READ(6,*) M,R
       WRITE(*,3)
   3   FORMAT(/)
       N=160
       NN=N+1
       H=2.0D0*DBLE(M)/DBLE(N)
       DINC=0.25D0
       ND=17
       DO 300 ID=1,ND
        D=0.0D0+DINC*(ID-1)
*
*     Apply trapezoidal rule to replace equation (1)
*     with a system of linear algebraic
*     equations A*X=B .
*
       DO 10 I=1,NN
        U=-M+H*DBLE(I-1)
        BI=DMIN1(M,U+R)
        AI=DMAX1(-M,U-R)
        N1=IFIX(SNGL((AI+M)/H)) + 1
        N2=IFIX(SNGL((BI+M)/H)) + 1
        B(I)=-1.0D0
        DO 20 J=1,NN
         ARG=AI+H*DBLE(J-N1)-D
         A(I,J)=H*P(ARG)
         IF (J .LT. N1) A(I,J)=0.0D0
         IF (J .EQ. N1) A(I,J)=H/2.0D0*P(AI-D)
         IF (J .EQ. N2) A(I,J)=H/2.0D0*P(BI-D)
         IF (J .GT. N2) A(I,J)=0.0D0
         IF (J .EQ. I) A(I,J)=A(I,J)-1.0D0
  20    CONTINUE
  10   CONTINUE
*
*    Calls to subroutines FACTOR and SUBST: Triangular
*    factorization and substitution to solve the
*    linear system A*X=B. The solution X is
*    returned by subroutine SUBST.
*
       CALL FACTOR(A,NN,D1,IPIVOT,IFLAG)
*
*    If IFLAG .NE. 0 then the linear system A*X=B can
*    be solved for X by CALL SUBST(A,IPIVOT,B,NN,X) .
*
       CALL SUBST(A,IPIVOT,B,NN,X)
*
*    Approximate ARL (equation (2)) using trapezoidal
*    rule and solution to linear system A*X=B.
*
       ARL=H/2.0D0*(P(-M-D)*X(1) + P(M-D)*X(NN))
       DO 30 I=2,N
        ARG=-M+H*DBLE(I-1)-D
        ARL=ARL + H*X(I)*P(ARG)
  30   CONTINUE
       ARL=ARL + 1.0D0
       WRITE(*,100) M,R,D,ARL
  100  FORMAT(10X,'M=',F5.2,2X,'R=',F5.2,2X,'D=',F5.2,
     1 2X,'ARL=',F10.4)
  300  CONTINUE
       STOP
       END
*
*    Evaluate the standard Normal density function .
*
       DOUBLE PRECISION FUNCTION P(X)
       DOUBLE PRECISION X,C0
       C0=6.283185307179586D0
       P=(1.0D0/DSQRT(C0))*DEXP(-.50D0*X**2)
       RETURN
       END
*
       SUBROUTINE SUBST(W,IPIVOT,B,N,X)
*
*   Source: Conte,S. D. and de Boor,C. (1980).
*   Elementary Numerical Analysis. McGraw-Hill,
*   New York, NY.
*
       INTEGER IPIVOT(N),I,IP,J
       REAL*8 B(N),W(N,N),X(N),SUM
       IF (N .LE. 1) THEN
        X(1)=B(1)/W(1,1)
        RETURN
       END IF
       IP=IPIVOT(1)
       X(1)=B(IP)
       DO 15 I=2,N
        SUM=0.0D0
        I1=I-1
        DO 14 J=1,I1
         SUM=W(I,J)*X(J) + SUM
  14    CONTINUE
        IP=IPIVOT(I)
        X(I)=B(IP) - SUM
  15   CONTINUE
       X(N)=X(N)/W(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=W(I,J)*X(J) + SUM
  19    CONTINUE
        X(I)=(X(I)-SUM)/W(I,I)
  20   CONTINUE
       RETURN
       END
*
       SUBROUTINE FACTOR(W,N,D,IPIVOT,IFLAG)
*
*   Source: Conte,S.D., and de Boor,C. (1980).
*   Elementary Numerical Analysis. McGraw-Hill,
*   New York, NY.
*
       REAL*8 D(N),W(N,N),AWIKOD,COLMAX,RATIO,
     1 ROWMAX,TEMP
       INTEGER IFLAG,IPIVOT(N),I,ISTAR,J,K
       IFLAG=1
       DO 9 I=1,N
        IPIVOT(I)=I
        ROWMAX=0.0D0
        DO 5 J=1,N
         ROWMAX=DMAX1(ROWMAX,DABS(W(I,J)))
    5   CONTINUE
        IF (ROWMAX .EQ. 0.0D0) THEN
         IFLAG=0
         ROWMAX=1.0D0
        END IF
        D(I)=ROWMAX
    9  CONTINUE
       IF (N .LE. 1) RETURN
       N1=N-1
       DO 20 K=1,N1
        COLMAX=DABS(W(K,K)/D(K))
        ISTAR=K
        K1=K+1
        DO 13 I=K1,N
         AWIKOD=DABS(W(I,K))/D(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=D(ISTAR)
          D(ISTAR)=D(K)
          D(K)=TEMP
          DO 15 J=1,N
           TEMP=W(ISTAR,J)
           W(ISTAR,J)=W(K,J)
           W(K,J)=TEMP
   15     CONTINUE
         END IF
         K2=K+1
         DO 19 I=K2,N
          W(I,K)=W(I,K)/W(K,K)
          RATIO=W(I,K)
          K3=K+1
          DO 18 J=K3,N
           W(I,J)=W(I,J)-RATIO*W(K,J)
   18     CONTINUE
   19    CONTINUE
         END IF
   20   CONTINUE
        IF (W(N,N) .EQ. 0.0D0) IFLAG=0
        RETURN
        END