real a(301,300),x(301),b(301),mflops real emax integer ipvt(300) c lda = 301 c write(6,40) 40 format(' lu decomposition timing') do 200 n = 50,300,50 call matgen(a,lda,n,b) write(6,50)lda,n 50 format(/' sp size of the arrays',i5,' and order is ',i5/) c t1 = second(t) call lu1(a,lda,n,ipvt,ierr) call lus(a,lda,n,ipvt,x,b) t2 = second(t) - t1 mflops = (float(n*(5+n*(-3+n*4)))/(6.0*t2))/1.0e6 emax = 0.0 do 73 i = 1,n emax = amax1(emax,abs(x(i))) 73 continue write(6,80) t2,mflops,emax 80 format(' unrolling depth 1 time = ',e12.3,' MFLOPS = ',f9.4, $ ' check = ',e12.3) c call matgen(a,lda,n,b) t1 = second(t) call lu2(a,lda,n,ipvt,ierr) call lus(a,lda,n,ipvt,x,b) t2 = second(t) - t1 mflops = (float(n*(5+n*(-3+n*4)))/(6.0*t2))/1.0e6 emax = 0.0 do 81 i = 1,n emax = amax1(emax,abs(x(i))) 81 continue write(6,90) t2,mflops,emax 90 format(' unrolling depth 2 time = ',e12.3,' MFLOPS = ',f9.4, $ ' check = ',e12.3) c call matgen(a,lda,n,b) t1 = second(t) call lu4(a,lda,n,ipvt,ierr) call lus(a,lda,n,ipvt,x,b) t2 = second(t) - t1 mflops = (float(n*(5+n*(-3+n*4)))/(6.0*t2))/1.0e6 emax = 0.0 do 91 i = 1,n emax = amax1(emax,abs(x(i))) 91 continue write(6,100) t2,mflops,emax 100 format(' unrolling depth 4 time = ',e12.3,' MFLOPS = ',f9.4, $ ' check = ',e12.3) c call matgen(a,lda,n,b) t1 = second(t) call lu8(a,lda,n,ipvt,ierr) call lus(a,lda,n,ipvt,x,b) t2 = second(t) - t1 mflops = (float(n*(5+n*(-3+n*4)))/(6.0*t2))/1.0e6 emax = 0.0 do 101 i = 1,n emax = amax1(emax,abs(x(i))) 101 continue write(6,110) t2,mflops,emax 110 format(' unrolling depth 8 time = ',e12.3,' MFLOPS = ',f9.4, $ ' check = ',e12.3) c call matgen(a,lda,n,b) t1 = second(t) call lux(a,lda,n,ipvt,ierr) call lus(a,lda,n,ipvt,x,b) t2 = second(t) - t1 mflops = (float(n*(5+n*(-3+n*4)))/(6.0*t2))/1.0e6 emax = 0.0 do 111 i = 1,n emax = amax1(emax,abs(x(i))) 111 continue write(6,120) t2,mflops,emax 120 format(' unrolling depth 16 time = ',e12.3,' MFLOPS = ',f9.4, $ ' check = ',e12.3) c 200 continue stop end subroutine matgen(a,lda,n,b) real a(lda,*),b(*) c init = 1325 do 30 j = 1,n do 20 i = 1,n init = mod(3125*init,65536) a(i,j) = (init - 32768.0)/16384.0 20 continue 30 continue do 35 i = 1,n b(i) = 0.0 35 continue do 50 j = 1,n do 40 i = 1,n b(i) = b(i) + a(i,j) 40 continue 50 continue return end SUBROUTINE LUS (A,LDA,N,IPVT,X,B) INTEGER IPVT(*) REAL A(LDA,*), X(*), B(*), XK C C PURPOSE: C Solve the linear system Ax = b given the LU factorization of A (as C computed in LU). C C ADDITIONAL PARAMETERS (NOT PARAMETERS OF LU): C C X REAL(N), solution of linear system C C B REAL(N), right-hand-side of linear system C C ---------------------------------------------------------------------- C DO 10 K = 1, N X(K) = B(K) 10 CONTINUE C DO 20 K = 1, N L = IPVT(K) XK = X(L) X(L) = X(K) X(K) = XK 20 CONTINUE C DO 40 K = 1, N XK = X(K)*A(K,K) DO 30 I = K+1, N X(I) = X(I) - A(I,K)*XK 30 CONTINUE X(K) = XK 40 CONTINUE C DO 60 K = N, 1, -1 XK = X(K) DO 50 I = 1, K-1 X(I) = X(I) + A(I,K)*XK 50 CONTINUE 60 CONTINUE C RETURN END SUBROUTINE LU1(A,LDA,N,IPVT,INFO) C ** LU.F -- LU DECOMPOSITION C C INTEGER LDA,N,IPVT(N),INFO REAL A(LDA,*),t C INFO = 0 DO 40 J = 1, N C C FORM J-TH COLUMN OF L C CALL SMXPY1(N-J+1,A(J,J),J-1,LDA,A(1,J),A(J,1)) C C SEARCH FOR PIVOT C T = ABS(A(J,J)) K = J DO 10 I = J+1, N IF (ABS(A(I,J)) .GT. T) THEN T = ABS(A(I,J)) K = I END IF 10 CONTINUE IPVT(J) = K C C TEST FOR ZERO PIVOT C IF (T .EQ. 0) THEN INFO = J GO TO 50 ENDIF C C INTERCHANGE ROWS C DO 20 I = 1, N T = A(J,I) A(J,I) = A(K,I) A(K,I) = T 20 CONTINUE C C FORM J-TH ROW OF U C A(J,J) = 1/A(J,J) CALL SXMPY1(N-J,LDA,A(J,J+1),J-1,LDA,A(J,1),LDA,A(1,J+1)) T = -A(J,J) DO 30 I = J+1, N A(J,I) = T*A(J,I) 30 CONTINUE 40 CONTINUE 50 RETURN END SUBROUTINE LU2(A,LDA,N,IPVT,INFO) C ** LU.F -- LU DECOMPOSITION C C INTEGER LDA,N,IPVT(N),INFO REAL A(LDA,*),t C INFO = 0 DO 40 J = 1, N C C FORM J-TH COLUMN OF L C CALL SMXPY2(N-J+1,A(J,J),J-1,LDA,A(1,J),A(J,1)) C C SEARCH FOR PIVOT C T = ABS(A(J,J)) K = J DO 10 I = J+1, N IF (ABS(A(I,J)) .GT. T) THEN T = ABS(A(I,J)) K = I END IF 10 CONTINUE IPVT(J) = K C C TEST FOR ZERO PIVOT C IF (T .EQ. 0) THEN INFO = J GO TO 50 ENDIF C C INTERCHANGE ROWS C DO 20 I = 1, N T = A(J,I) A(J,I) = A(K,I) A(K,I) = T 20 CONTINUE C C FORM J-TH ROW OF U C A(J,J) = 1/A(J,J) CALL SXMPY2(N-J,LDA,A(J,J+1),J-1,LDA,A(J,1),LDA,A(1,J+1)) T = -A(J,J) DO 30 I = J+1, N A(J,I) = T*A(J,I) 30 CONTINUE 40 CONTINUE 50 RETURN END SUBROUTINE LU4(A,LDA,N,IPVT,INFO) C ** LU.F -- LU DECOMPOSITION C C INTEGER LDA,N,IPVT(N),INFO REAL A(LDA,*),t C INFO = 0 DO 40 J = 1, N C C FORM J-TH COLUMN OF L C CALL SMXPY4(N-J+1,A(J,J),J-1,LDA,A(1,J),A(J,1)) C C SEARCH FOR PIVOT C T = ABS(A(J,J)) K = J DO 10 I = J+1, N IF (ABS(A(I,J)) .GT. T) THEN T = ABS(A(I,J)) K = I END IF 10 CONTINUE IPVT(J) = K C C TEST FOR ZERO PIVOT C IF (T .EQ. 0) THEN INFO = J GO TO 50 ENDIF C C INTERCHANGE ROWS C DO 20 I = 1, N T = A(J,I) A(J,I) = A(K,I) A(K,I) = T 20 CONTINUE C C FORM J-TH ROW OF U C A(J,J) = 1/A(J,J) CALL SXMPY4(N-J,LDA,A(J,J+1),J-1,LDA,A(J,1),LDA,A(1,J+1)) T = -A(J,J) DO 30 I = J+1, N A(J,I) = T*A(J,I) 30 CONTINUE 40 CONTINUE 50 RETURN END SUBROUTINE LU8(A,LDA,N,IPVT,INFO) C ** LU.F -- LU DECOMPOSITION C C INTEGER LDA,N,IPVT(N),INFO REAL A(LDA,*),t C INFO = 0 DO 40 J = 1, N C C FORM J-TH COLUMN OF L C CALL SMXPY8(N-J+1,A(J,J),J-1,LDA,A(1,J),A(J,1)) C C SEARCH FOR PIVOT C T = ABS(A(J,J)) K = J DO 10 I = J+1, N IF (ABS(A(I,J)) .GT. T) THEN T = ABS(A(I,J)) K = I END IF 10 CONTINUE IPVT(J) = K C C TEST FOR ZERO PIVOT C IF (T .EQ. 0) THEN INFO = J GO TO 50 ENDIF C C INTERCHANGE ROWS C DO 20 I = 1, N T = A(J,I) A(J,I) = A(K,I) A(K,I) = T 20 CONTINUE C C FORM J-TH ROW OF U C A(J,J) = 1/A(J,J) CALL SXMPY8(N-J,LDA,A(J,J+1),J-1,LDA,A(J,1),LDA,A(1,J+1)) T = -A(J,J) DO 30 I = J+1, N A(J,I) = T*A(J,I) 30 CONTINUE 40 CONTINUE 50 RETURN END SUBROUTINE LUX(A,LDA,N,IPVT,INFO) C ** LU.F -- LU DECOMPOSITION C C INTEGER LDA,N,IPVT(N),INFO REAL A(LDA,*),t C INFO = 0 DO 40 J = 1, N C C FORM J-TH COLUMN OF L C CALL SMXPYX(N-J+1,A(J,J),J-1,LDA,A(1,J),A(J,1)) C C SEARCH FOR PIVOT C T = ABS(A(J,J)) K = J DO 10 I = J+1, N IF (ABS(A(I,J)) .GT. T) THEN T = ABS(A(I,J)) K = I END IF 10 CONTINUE IPVT(J) = K C C TEST FOR ZERO PIVOT C IF (T .EQ. 0) THEN INFO = J GO TO 50 ENDIF C C INTERCHANGE ROWS C DO 20 I = 1, N T = A(J,I) A(J,I) = A(K,I) A(K,I) = T 20 CONTINUE C C FORM J-TH ROW OF U C A(J,J) = 1/A(J,J) CALL SXMPYX(N-J,LDA,A(J,J+1),J-1,LDA,A(J,1),LDA,A(1,J+1)) T = -A(J,J) DO 30 I = J+1, N A(J,I) = T*A(J,I) 30 CONTINUE 40 CONTINUE 50 RETURN END SUBROUTINE SMXPY1 (N1,Y,N2,LDM,X,M) C ** SMXPY01.F -- SMXPY UNROLLED TO A DEPTH OF 1 C INTEGER LDM,N1,N2 REAL Y(*),X(*),M(LDM,*) C REAL Y(N1),X(N2),M(LDM,N2) C DO 20 J = 1, N2 DO 10 I = 1, N1 Y(I) = (Y(I)) + X(J)*M(I,J) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE SMXPY2(N1,Y,N2,LDM,X,M) C ** SMXPY02.F -- SMXPY UNROLLED TO A DEPTH OF 2 C INTEGER LDM,N1,N2 REAL Y(*),X(*),M(LDM,*) C REAL Y(N1),X(N2),M(LDM,N2) C C Cleanup odd vector C J = MOD(N2,2) IF (J .GE. 1) THEN DO 10 I = 1, N1 Y(I) = (Y(I)) + X(J)*M(I,J) 10 CONTINUE ENDIF C C Main loop - groups of two vectors C JMIN = J+2 DO 30 J = JMIN, N2, 2 DO 20 I = 1, N1 Y(I) = ( (Y(I)) $ + X(J-1)*M(I,J-1)) + X(J)*M(I,J) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE SMXPY4(N1,Y,N2,LDM,X,M) C ** SMXPY04.F -- SMXPY UNROLLED TO A DEPTH OF 4 C INTEGER LDM,N1,N2 REAL Y(*),X(*),M(LDM,*) C REAL Y(N1),X(N2),M(LDM,N2) C C Cleanup odd vector C J = MOD(N2,2) IF (J .GE. 1) THEN DO 10 I = 1, N1 Y(I) = (Y(I)) + X(J)*M(I,J) 10 CONTINUE ENDIF C C Cleanup odd group of two vectors C J = MOD(N2,4) IF (J .GE. 2) THEN DO 20 I = 1, N1 Y(I) = ( (Y(I)) $ + X(J-1)*M(I,J-1)) + X(J)*M(I,J) 20 CONTINUE ENDIF C C Main loop - groups of four vectors C JMIN = J+4 DO 40 J = JMIN, N2, 4 DO 30 I = 1, N1 Y(I) = ((( (Y(I)) $ + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2)) $ + X(J-1)*M(I,J-1)) + X(J) *M(I,J) 30 CONTINUE 40 CONTINUE RETURN END SUBROUTINE SMXPY8(N1,Y,N2,LDM,X,M) C ** SMXPY08.F -- SMXPY UNROLLED TO A DEPTH OF 8 C INTEGER LDM,N1,N2 REAL Y(*),X(*),M(LDM,*) C REAL Y(N1),X(N2),M(LDM,N2) C C Cleanup odd vector C J = MOD(N2,2) IF (J .GE. 1) THEN DO 10 I = 1, N1 Y(I) = (Y(I)) + X(J)*M(I,J) 10 CONTINUE ENDIF C C Cleanup odd group of two vectors C J = MOD(N2,4) IF (J .GE. 2) THEN DO 20 I = 1, N1 Y(I) = ( (Y(I)) $ + X(J-1)*M(I,J-1)) + X(J)*M(I,J) 20 CONTINUE ENDIF C C Cleanup odd group of four vectors C J = MOD(N2,8) IF (J .GE. 4) THEN DO 30 I = 1, N1 Y(I) = ((( (Y(I)) $ + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2)) $ + X(J-1)*M(I,J-1)) + X(J) *M(I,J) 30 CONTINUE ENDIF C C Main loop - groups of eight vectors C JMIN = J+8 DO 50 J = JMIN, N2, 8 DO 40 I = 1, N1 Y(I) = ((((((( (Y(I)) $ + X(J-7)*M(I,J-7)) + X(J-6)*M(I,J-6)) $ + X(J-5)*M(I,J-5)) + X(J-4)*M(I,J-4)) $ + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2)) $ + X(J-1)*M(I,J-1)) + X(J) *M(I,J) 40 CONTINUE 50 CONTINUE RETURN END SUBROUTINE SMXPYX(N1,Y,N2,LDM,X,M) C ** SMXPY16.F -- SMXPY UNROLLED TO A DEPTH OF 16 C INTEGER LDM,N1,N2 REAL Y(*),X(*),M(LDM,*) C REAL Y(N1),X(N2),M(LDM,N2) C C Cleanup odd vector C J = MOD(N2,2) IF (J .GE. 1) THEN DO 10 I = 1, N1 Y(I) = (Y(I)) + X(J)*M(I,J) 10 CONTINUE ENDIF C C Cleanup odd group of two vectors C J = MOD(N2,4) IF (J .GE. 2) THEN DO 20 I = 1, N1 Y(I) = ( (Y(I)) $ + X(J-1)*M(I,J-1)) + X(J)*M(I,J) 20 CONTINUE ENDIF C C Cleanup odd group of four vectors C J = MOD(N2,8) IF (J .GE. 4) THEN DO 30 I = 1, N1 Y(I) = ((( (Y(I)) $ + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2)) $ + X(J-1)*M(I,J-1)) + X(J) *M(I,J) 30 CONTINUE ENDIF C C Cleanup odd group of eight vectors C J = MOD(N2,16) IF (J .GE. 8) THEN DO 40 I = 1, N1 Y(I) = ((((((( (Y(I)) $ + X(J-7)*M(I,J-7)) + X(J-6)*M(I,J-6)) $ + X(J-5)*M(I,J-5)) + X(J-4)*M(I,J-4)) $ + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2)) $ + X(J-1)*M(I,J-1)) + X(J) *M(I,J) 40 CONTINUE ENDIF C C Main loop - groups of sixteen vectors C JMIN = J+16 DO 60 J = JMIN, N2, 16 DO 50 I = 1, N1 Y(I) = ((((((((((((((( (Y(I)) $ + X(J-15)*M(I,J-15)) + X(J-14)*M(I,J-14)) $ + X(J-13)*M(I,J-13)) + X(J-12)*M(I,J-12)) $ + X(J-11)*M(I,J-11)) + X(J-10)*M(I,J-10)) $ + X(J- 9)*M(I,J- 9)) + X(J- 8)*M(I,J- 8)) $ + X(J- 7)*M(I,J- 7)) + X(J- 6)*M(I,J- 6)) $ + X(J- 5)*M(I,J- 5)) + X(J- 4)*M(I,J- 4)) $ + X(J- 3)*M(I,J- 3)) + X(J- 2)*M(I,J- 2)) $ + X(J- 1)*M(I,J- 1)) + X(J) *M(I,J) 50 CONTINUE 60 CONTINUE RETURN END SUBROUTINE SXMPY1(N1,LDY,Y,N2,LDX,X,LDM,M) C ** SXMPY01.F -- SXMPY UNROLLED TO A DEPTH OF 1 C INTEGER N1,LDY,N2,LDX,LDM REAL Y(LDY,*),X(LDX,*),M(LDM,*) C REAL Y(LDY,N1),X(LDX,N2),M(LDM,N1) C DO 20 J = 1, N2 DO 10 I = 1, N1 Y(1,I) = (Y(1,I)) + X(1,J)*M(J,I) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE SXMPY2(N1,LDY,Y,N2,LDX,X,LDM,M) C ** SXMPY02.F -- SXMPY UNROLLED TO A DEPTH OF 2 C INTEGER N1,LDY,N2,LDX,LDM REAL Y(LDY,*),X(LDX,*),M(LDM,*) C REAL Y(LDY,N1),X(LDX,N2),M(LDM,N1) C C Cleanup odd vector C J = MOD(N2,2) IF (J .GE. 1) THEN DO 10 I = 1, N1 Y(1,I) = (Y(1,I)) + X(1,J)*M(J,I) 10 CONTINUE ENDIF C C Main loop - groups of two vectors C JMIN = J+2 DO 30 J = JMIN, N2, 2 DO 20 I = 1, N1 Y(1,I) = ( (Y(1,I)) $ + X(1,J-1)*M(J-1,I)) + X(1,J)*M(J,I) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE SXMPY4(N1,LDY,Y,N2,LDX,X,LDM,M) C ** SXMPY04.F -- SXMPY UNROLLED TO A DEPTH OF 4 C INTEGER N1,LDY,N2,LDX,LDM REAL Y(LDY,*),X(LDX,*),M(LDM,*) C REAL Y(LDY,N1),X(LDX,N2),M(LDM,N1) C C Cleanup odd vector C J = MOD(N2,2) IF (J .GE. 1) THEN DO 10 I = 1, N1 Y(1,I) = (Y(1,I)) + X(1,J)*M(J,I) 10 CONTINUE ENDIF C C Cleanup odd group of two vectors C J = MOD(N2,4) IF (J .GE. 2) THEN DO 20 I = 1, N1 Y(1,I) = ( (Y(1,I)) $ + X(1,J-1)*M(J-1,I)) + X(1,J)*M(J,I) 20 CONTINUE ENDIF C C Main loop - groups of four vectors C JMIN = J+4 DO 40 J = JMIN, N2, 4 DO 30 I = 1, N1 Y(1,I) = ((( (Y(1,I)) $ + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I)) $ + X(1,J-1)*M(J-1,I)) + X(1,J) *M(J,I) 30 CONTINUE 40 CONTINUE RETURN END SUBROUTINE SXMPY8(N1,LDY,Y,N2,LDX,X,LDM,M) C ** SXMPY08.F -- SXMPY UNROLLED TO A DEPTH OF 8 C INTEGER N1,LDY,N2,LDX,LDM REAL Y(LDY,*),X(LDX,*),M(LDM,*) C REAL Y(LDY,N1),X(LDX,N2),M(LDM,N1) C C Cleanup odd vector C J = MOD(N2,2) IF (J .GE. 1) THEN DO 10 I = 1, N1 Y(1,I) = (Y(1,I)) + X(1,J)*M(J,I) 10 CONTINUE ENDIF C C Cleanup odd group of two vectors C J = MOD(N2,4) IF (J .GE. 2) THEN DO 20 I = 1, N1 Y(1,I) = ( (Y(1,I)) $ + X(1,J-1)*M(J-1,I)) + X(1,J)*M(J,I) 20 CONTINUE ENDIF C C Cleanup odd group of four vectors C J = MOD(N2,8) IF (J .GE. 4) THEN DO 30 I = 1, N1 Y(1,I) = ((( (Y(1,I)) $ + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I)) $ + X(1,J-1)*M(J-1,I)) + X(1,J) *M(J,I) 30 CONTINUE ENDIF C C Main loop - groups of eight vectors C JMIN = J+8 DO 50 J = JMIN, N2, 8 DO 40 I = 1, N1 Y(1,I) = ((((((( (Y(1,I)) $ + X(1,J-7)*M(J-7,I)) + X(1,J-6)*M(J-6,I)) $ + X(1,J-5)*M(J-5,I)) + X(1,J-4)*M(J-4,I)) $ + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I)) $ + X(1,J-1)*M(J-1,I)) + X(1,J) *M(J,I) 40 CONTINUE 50 CONTINUE RETURN END SUBROUTINE SXMPYX(N1,LDY,Y,N2,LDX,X,LDM,M) C ** SXMPY16.F -- SXMPY UNROLLED TO A DEPTH OF 16 C INTEGER N1,LDY,N2,LDX,LDM REAL Y(LDY,*),X(LDX,*),M(LDM,*) C REAL Y(LDY,N1),X(LDX,N2),M(LDM,N1) C C Cleanup odd vector C J = MOD(N2,2) IF (J .GE. 1) THEN DO 10 I = 1, N1 Y(1,I) = (Y(1,I)) + X(1,J)*M(J,I) 10 CONTINUE ENDIF C C Cleanup odd group of two vectors C J = MOD(N2,4) IF (J .GE. 2) THEN DO 20 I = 1, N1 Y(1,I) = ( (Y(1,I)) $ + X(1,J-1)*M(J-1,I)) + X(1,J)*M(J,I) 20 CONTINUE ENDIF C C Cleanup odd group of four vectors C J = MOD(N2,8) IF (J .GE. 4) THEN DO 30 I = 1, N1 Y(1,I) = ((( (Y(1,I)) $ + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I)) $ + X(1,J-1)*M(J-1,I)) + X(1,J) *M(J,I) 30 CONTINUE ENDIF C C Cleanup odd group of eight vectors C J = MOD(N2,16) IF (J .GE. 8) THEN DO 40 I = 1, N1 Y(1,I) = ((((((( (Y(1,I)) $ + X(1,J-7)*M(J-7,I)) + X(1,J-6)*M(J-6,I)) $ + X(1,J-5)*M(J-5,I)) + X(1,J-4)*M(J-4,I)) $ + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I)) $ + X(1,J-1)*M(J-1,I)) + X(1,J) *M(J,I) 40 CONTINUE ENDIF C C Main loop - groups of sixteen vectors C JMIN = J+16 DO 60 J = JMIN, N2, 16 DO 50 I = 1, N1 Y(1,I) = ((((((((((((((( (Y(1,I)) $ + X(1,J-15)*M(J-15,I)) + X(1,J-14)*M(J-14,I)) $ + X(1,J-13)*M(J-13,I)) + X(1,J-12)*M(J-12,I)) $ + X(1,J-11)*M(J-11,I)) + X(1,J-10)*M(J-10,I)) $ + X(1,J- 9)*M(J- 9,I)) + X(1,J- 8)*M(J- 8,I)) $ + X(1,J- 7)*M(J- 7,I)) + X(1,J- 6)*M(J- 6,I)) $ + X(1,J- 5)*M(J- 5,I)) + X(1,J- 4)*M(J- 4,I)) $ + X(1,J- 3)*M(J- 3,I)) + X(1,J- 2)*M(J- 2,I)) $ + X(1,J- 1)*M(J- 1,I)) + X(1,J) *M(J,I) 50 CONTINUE 60 CONTINUE RETURN END