**************************************************** * * 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(*,*) 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 