c Date: Fri, 24 Jun 88 11:16:54 PDT c From: David Bailey PROGRAM NASKER C C NAS KERNEL BENCHMARK PROGRAM C 12/17/84 DAVID H BAILEY C CHARACTER*8 PN(8) REAL ER(8), FP(8), TM(8), RT(8) COMMON /ARRAYS/ DATA(360000) DATA PN/'MXM', 'CFFT2D', 'CHOLSKY', 'BTRIX', 'GMTRY', 'EMIT', $ 'VPENTA', 'TOTAL'/ C WRITE (6, 1) 1 FORMAT (/16X, 'THE NAS KERNEL BENCHMARK PROGRAM'//) C CALL MXMTST (ER(1), FP(1), TM(1)) CALL FFTTST (ER(2), FP(2), TM(2)) CALL CHOTST (ER(3), FP(3), TM(3)) CALL BTRTST (ER(4), FP(4), TM(4)) CALL GMTTST (ER(5), FP(5), TM(5)) CALL EMITST (ER(6), FP(6), TM(6)) CALL VPETST (ER(7), FP(7), TM(7)) C TE = 0. TF = 0. TT = 0. DO 100 I = 1, 7 TE = TE + ER(I) TF = TF + FP(I) TT = TT + TM(I) RT(I) = 1E-6 * FP(I) / TM(I) 100 CONTINUE ER(8) = TE FP(8) = TF TM(8) = TT RT(8) = 1E-6 * TF / TT C WRITE (6, 2) (PN(I), ER(I), FP(I), TM(I), RT(I), I = 1, 8) 2 FORMAT (' PROGRAM', 8X, 'ERROR', 10X, 'FP OPS', 7X, 'SECONDS', $ 6X, 'MFLOPS'// 7(1X, A8, 1P2E15.4, 0PF12.4, F12.2/)/ $ 1X, A8, 1P2E15.4, 0PF12.4, F12.2//) C STOP END C FUNCTION CPTIME () C C RETURNS THE CPU TIME SINCE THE LAST CALL TO CPTIME. C THIS SUBPROGRAM MAY BE CHANGED AS NEEDED FOR A PARTICULAR COMPUTER C SYSTEM WITHOUT PENALTY, PROVIDED IT PERFORMS THIS FUNCTION. C DATA TX/0./ T = SECOND () CPTIME = T - TX TX = T RETURN END C SUBROUTINE COPY (N, A, B) C C ARRAY COPY ROUTINE C REAL A(N), B(N) DO 100 I = 1, N B(I) = A(I) 100 CONTINUE RETURN END C SUBROUTINE MXMTST (ER, FP, TM) C C FLOATING-POINT MATRIX MULTIPLY TEST C PARAMETER (L=256, M=128, N=64, F7=78125., T30=1073741824.) COMMON /ARRAYS/ A(L,M), S1, B(M,N), S2, C(L,N) DATA IT/100/, ANS/35.2026179738722/ C C INITIALIZATION C C THE ARRAYS A AND B ARE FILLED WITH PSEUDO-RANDOM (0., 1.) DATA C USING A RANDOM NUMBER GENERATOR BASED ON THE RECURSION C X(N+1) = 5**7 * X(N) (MOD 2**30) C THIS RECURSION WILL GENERATE 2**28 (APPROX. 268 MILLION) NUMBERS C BEFORE REPEATING. FOR THIS SCHEME TO WORK PROPERLY, THE HARDWARE C MULTIPLY OPERATION MUST BE CORRECT TO 47 BITS OF PRECISION. C THIS SAME SCHEME IS USED TO INITIALIZE DATA ARRAYS FOR ALL TESTS. C T = F7 / T30 DO 100 J = 1, M DO 100 I = 1, L T = MOD (F7 * T, 1.) A(I,J) = T 100 CONTINUE DO 110 J = 1, N DO 110 I = 1, M T = MOD (F7 * T, 1.) B(I,J) = T 110 CONTINUE TM = CPTIME () C C TIMING TEST C DO 120 II = 1, IT CALL MXM (A, B, C, L, M, N) 120 CONTINUE C TM = CPTIME () ER = ABS ((C(19,19) - ANS) / ANS) FP = 2. * IT * L * M * N C RETURN END C SUBROUTINE MXM (A, B, C, L, M, N) DIMENSION A(L,M), B(M,N), C(L,N) C C 4-WAY UNROLLED MATRIX MULTIPLY ROUTINE FOR VECTOR COMPUTERS. C M MUST BE A MULTIPLE OF 4. CONTIGUOUS DATA ASSUMED. C D H BAILEY 11/15/84 C DO 100 K = 1, N DO 100 I = 1, L C(I,K) = 0. 100 CONTINUE DO 110 J = 1, M, 4 DO 110 K = 1, N DO 110 I = 1, L C(I,K) = C(I,K) + A(I,J) * B(J,K) $ + A(I,J+1) * B(J+1,K) + A(I,J+2) * B(J+2,K) $ + A(I,J+3) * B(J+3,K) 110 CONTINUE C RETURN END C SUBROUTINE FFTTST (ER, FP, TM) C C 2-D FFT TEST PROGRAM C PARAMETER (M=128, N=256, M1=128, F7=78125., T30=1073741824.) COMPLEX X, Y, CT COMMON /ARRAYS/ X(M1,N), W1(M), W2(N), IP(2*N) DATA IT/100/, ANS/0.894799941219277/ C C INITIALIZE C AMN = M * N RMN = 1. / AMN T2 = F7 / T30 DO 100 J = 1, N DO 100 I = 1, M T1 = MOD (F7 * T2, 1.) T2 = MOD (F7 * T1, 1.) X(I,J) = CMPLX (T1, T2) 100 CONTINUE CALL CFFT2D1 (0, M, M1, N, X, W1, IP) CALL CFFT2D2 (0, M, M1, N, X, W2, IP) TM = CPTIME () C C TEST ITERATIONS C DO 120 K = 1, IT DO 110 J = 1, N DO 110 I = 1, M X(I,J) = RMN * X(I,J) 110 CONTINUE C CALL CFFT2D1 (1, M, M1, N, X, W1, IP) CALL CFFT2D2 (1, M, M1, N, X, W2, IP) CALL CFFT2D2 (-1, M, M1, N, X, W2, IP) CALL CFFT2D1 (-1, M, M1, N, X, W1, IP) 120 CONTINUE C TM = CPTIME () ER = ABS ((REAL(X(19,19)) - ANS) / ANS) FP = IT * AMN * (2. + 10. * LOG (AMN)/LOG (2.)) C RETURN END C SUBROUTINE CFFT2D1 (IS, M, M1, N, X, W, IP) C C PERFORMS COMPLEX RADIX 2 FFTS ON THE FIRST DIMENSION OF THE 2-D ARRAY X C D H BAILEY 11/15/84 C COMPLEX X(M1,N), W(M), CT, CX INTEGER IP(2,M) DATA PI/3.141592653589793/ C C IF IS = 0 THEN INITIALIZE ONLY C M2 = M / 2 IF (IS .EQ. 0) THEN DO 100 I = 1, M2 T = 2. * PI * (I-1) / M W(I) = CMPLX (COS (T), SIN (T)) 100 CONTINUE RETURN ENDIF C C PERFORM FORWARD OR BACKWARD FFTS ACCORDING TO IS = 1 OR -1 C DO 110 I = 1, M IP(1,I) = I 110 CONTINUE L = 1 I1 = 1 C 120 I2 = 3 - I1 DO 130 J = L, M2, L CX = W(J-L+1) IF (IS .LT. 0) CX = CONJG (CX) DO 130 I = J-L+1, J II = IP(I1,I) IP(I2,I+J-L) = II IM = IP(I1,I+M2) IP(I2,I+J) = IM DO 130 K = 1, N CT = X(II,K) - X(IM,K) X(II,K) = X(II,K) + X(IM,K) X(IM,K) = CT * CX 130 CONTINUE L = 2 * L I1 = I2 IF (L .LE. M2) GOTO 120 C DO 150 I = 1, M II = IP(I1,I) IF (II .GT. I) THEN DO 140 K = 1, N CT = X(I,K) X(I,K) = X(II,K) X(II,K) = CT 140 CONTINUE ENDIF 150 CONTINUE C RETURN END C SUBROUTINE CFFT2D2 (IS, M, M1, N, X, W, IP) C C PERFORMS COMPLEX RADIX 2 FFTS ON THE SECOND DIMENSION OF THE 2-D ARRAY X C D H BAILEY 11/15/84 C COMPLEX X(M1,N), W(N), CT, CX INTEGER IP(2,N) DATA PI/3.141592653589793/ C C IF IS = 0 THEN INITIALIZE ONLY C N2 = N / 2 IF (IS .EQ. 0) THEN DO 100 I = 1, N2 T = 2. * PI * (I-1) / N W(I) = CMPLX (COS (T), SIN (T)) 100 CONTINUE RETURN ENDIF C C PEFORM FORWARD OR BACKWARD FFTS ACCORDING TO IS = 1 OR -1 C DO 110 I = 1, N IP(1,I) = I 110 CONTINUE L = 1 I1 = 1 C 120 I2 = 3 - I1 DO 130 J = L, N2, L CX = W(J-L+1) IF (IS .LT. 0) CX = CONJG (CX) DO 130 I = J-L+1, J II = IP(I1,I) IP(I2,I+J-L) = II IM = IP(I1,I+N2) IP(I2,I+J) = IM DO 130 K = 1, M CT = X(K,II) - X(K,IM) X(K,II) = X(K,II) + X(K,IM) X(K,IM) = CT * CX 130 CONTINUE L = 2 * L I1 = I2 IF (L .LE. N2) GOTO 120 C DO 150 I = 1, N II = IP(I1,I) IF (II .GT. I) THEN DO 140 K = 1, M CT = X(K,I) X(K,I) = X(K,II) X(K,II) = CT 140 CONTINUE ENDIF 150 CONTINUE C RETURN END C SUBROUTINE CHOTST (ER, FP, TM) C C CHOLSKY TEST PROGRAM C PARAMETER (IDA=250, NMAT=250, M=4, N=40, NRHS=3, F7=78125., $ T30=1073741824.) COMMON /ARRAYS/ A(0:IDA, -M:0, 0:N), B(0:NRHS, 0:NMAT, 0:N), $ AX(0:IDA, -M:0, 0:N), BX(0:NRHS, 0:NMAT, 0:N) DATA IT/200/, ANS/5177.88531774562/ C C INITIALIZE C LA = (IDA+1) * (M+1) * (N+1) LB = (NRHS+1) * (NMAT+1) * (N+1) T = F7 / T30 DO 100 K = 0, N DO 100 J = -M, 0 DO 100 I = 0, IDA T = MOD (F7 * T, 1.) AX(I,J,K) = T 100 CONTINUE DO 110 K = 0, N DO 110 J = 0, NMAT DO 110 I = 0, NRHS T = MOD (F7 * T, 1.) BX(I,J,K) = T 110 CONTINUE TM = CPTIME () C C BEGIN TIMING TEST C DO 120 J = 1, IT CALL COPY (LA, AX, A) CALL COPY (LB, BX, B) CALL CHOLSKY (IDA, NMAT, M, N, A, NRHS, IDA, B) 120 CONTINUE C TM = CPTIME () ER = ABS ((B(1,19,19) - ANS) / ANS) FP = IT * (NMAT + 1) * 4403. C RETURN END C SUBROUTINE CHOLSKY (IDA, NMAT, M, N, A, NRHS, IDB, B) C C CHOLESKY DECOMPOSITION/SUBSTITUTION SUBROUTINE. C C 11/28/84 D H BAILEY MODIFIED FOR NAS KERNEL TEST C REAL A(0:IDA, -M:0, 0:N), B(0:NRHS, 0:IDB, 0:N), EPSS(0:256) DATA EPS/1E-13/ C C CHOLESKY DECOMPOSITION C DO 1 J = 0, N I0 = MAX ( -M, -J ) C C OFF DIAGONAL ELEMENTS C DO 2 I = I0, -1 DO 3 JJ = I0 - I, -1 DO 3 L = 0, NMAT 3 A(L,I,J) = A(L,I,J) - A(L,JJ,I+J) * A(L,I+JJ,J) DO 2 L = 0, NMAT 2 A(L,I,J) = A(L,I,J) * A(L,0,I+J) C C STORE INVERSE OF DIAGONAL ELEMENTS C DO 4 L = 0, NMAT 4 EPSS(L) = EPS * A(L,0,J) DO 5 JJ = I0, -1 DO 5 L = 0, NMAT 5 A(L,0,J) = A(L,0,J) - A(L,JJ,J) ** 2 DO 1 L = 0, NMAT 1 A(L,0,J) = 1. / SQRT ( ABS (EPSS(L) + A(L,0,J)) ) C C SOLUTION C DO 6 I = 0, NRHS DO 7 K = 0, N DO 8 L = 0, NMAT 8 B(I,L,K) = B(I,L,K) * A(L,0,K) DO 7 JJ = 1, MIN (M, N-K) DO 7 L = 0, NMAT 7 B(I,L,K+JJ) = B(I,L,K+JJ) - A(L,-JJ,K+JJ) * B(I,L,K) C DO 6 K = N, 0, -1 DO 9 L = 0, NMAT 9 B(I,L,K) = B(I,L,K) * A(L,0,K) DO 6 JJ = 1, MIN (M, K) DO 6 L = 0, NMAT 6 B(I,L,K-JJ) = B(I,L,K-JJ) - A(L,-JJ,K) * B(I,L,K) C RETURN END C SUBROUTINE BTRTST (ER, FP, TM) C C BTRIX TEST PROGRAM C PARAMETER (JD=30, KD=30, LD=30, MD=30, F7=78125., T30=1073741824.) COMMON /ARRAYS/ S(JD,KD,LD,5), A(5,5,MD,MD), B(5,5,MD,MD), $ C(5,5,MD,MD), SX(JD,KD,LD,5), BX(5,5,MD,MD) DATA JS/2/, JE/29/, LS/2/, LE/29/, IT/20/, ANS/-0.286282658663962/ C C INITIALIZATION C NB = 25 * MD * MD NS = JD * KD * LD * 5 T = F7 / T30 DO 100 L = 1, MD DO 100 K = 1, MD DO 100 J = 1, 5 DO 100 I = 1, 5 T = MOD (F7 * T, 1.) A(I,J,K,L) = T T = MOD (F7 * T, 1.) BX(I,J,K,L) = T T = MOD (F7 * T, 1.) C(I,J,K,L) = T 100 CONTINUE DO 110 L = 1, 5 DO 110 K = 1, LD DO 110 J = 1, KD DO 110 I = 1, JD T = MOD (F7 * T, 1.) SX(I,J,K,L) = T 110 CONTINUE TM = CPTIME () C C TIMING TEST C DO 120 II = 1, IT CALL COPY (NS, SX, S) DO 120 K = 1, KD CALL COPY (NB, BX, B) CALL BTRIX (JS, JE, LS, LE, K) 120 CONTINUE C TM = CPTIME () ER = ABS ((S(19,19,19,1) - ANS) / ANS) FP = IT * MD * (LE - 1) * 19165. C RETURN END C SUBROUTINE BTRIX (JS, JE, LS, LE, K) C C VECTORIZED BLOCK TRI-DIAGONAL SOLVER IN THE J DIRECTION C FOR K = CONSTANT PLANES C C 11/15/84 D H BAILEY MODIFIED FOR NAS KERNEL TEST C PARAMETER (JD=30, KD=30, LD=30, MD=30) COMMON /ARRAYS/ S(JD,KD,LD,5), A(5,5,MD,MD), B(5,5,MD,MD), $ C(5,5,MD,MD) C DIMENSION U12(MD), U13(MD), U14(MD), U15(MD), U23(MD), $ U24(MD), U25(MD), U34(MD), U35(MD), U45(MD) C REAL L11(MD), L21(MD), L31(MD), L41(MD), L51(MD), $ L22(MD), L32(MD), L42(MD), L52(MD), L33(MD), $ L43(MD), L53(MD), L44(MD), L54(MD), L55(MD) C C PART 1. FORWARD BLOCK SWEEP C C DO 100 J = JS,JE C C********** STEP 1. CONSTRUCT L(I) IN B ************************** C IF(J.EQ.JS) GO TO 4 DO 3 M = 1,5 DO 3 N = 1,5 DO 3 L = LS,LE B(M,N,J,L) = B(M,N,J,L) - A(M,1,J,L)*B(1,N,J-1,L) $ - A(M,2,J,L)*B(2,N,J-1,L) - A(M,3,J,L)*B(3,N,J-1,L) $ - A(M,4,J,L)*B(4,N,J-1,L) - A(M,5,J,L)*B(5,N,J-1,L) 3 CONTINUE C 4 CONTINUE C C********** STEP 2. CONPUTE L INVERSE *************************** C C C A. DECOMPOSE L(I) INTO L AND U C DO 20 L = LS,LE L11(L) = 1. / B(1,1,J,L) U12(L) = B(1,2,J,L)*L11(L) U13(L) = B(1,3,J,L)*L11(L) U14(L) = B(1,4,J,L)*L11(L) U15(L) = B(1,5,J,L)*L11(L) L21(L) = B(2,1,J,L) L22(L) = 1. / (B(2,2,J,L) - L21(L)*U12(L)) U23(L) = (B(2,3,J,L) - L21(L)*U13(L)) * L22(L) U24(L) = (B(2,4,J,L) - L21(L)*U14(L)) * L22(L) U25(L) = (B(2,5,J,L) - L21(L)*U15(L)) * L22(L) L31(L) = B(3,1,J,L) L32(L) = B(3,2,J,L) - L31(L)*U12(L) L33(L) = 1. / (B(3,3,J,L) - L31(L)*U13(L) - L32(L)*U23(L)) U34(L) = (B(3,4,J,L) - L31(L)*U14(L) - L32(L)*U24(L)) * L33(L) U35(L) = (B(3,5,J,L) - L31(L)*U15(L) - L32(L)*U25(L)) * L33(L) 20 CONTINUE C DO 25 L = LS,LE L41(L) = B(4,1,J,L) L42(L) = B(4,2,J,L) - L41(L)*U12(L) L43(L) = B(4,3,J,L) - L41(L)*U13(L) - L42(L)*U23(L) L44(L) = 1. / (B(4,4,J,L) - L41(L)*U14(L) - L42(L)*U24(L) $ - L43(L)*U34(L)) U45(L) = (B(4,5,J,L) - L41(L)*U15(L) - L42(L)*U25(L) $ - L43(L)*U35(L)) * L44(L) L51(L) = B(5,1,J,L) L52(L) = B(5,2,J,L) - L51(L)*U12(L) L53(L) = B(5,3,J,L) - L51(L)*U13(L) - L52(L)*U23(L) L54(L) = B(5,4,J,L) - L51(L)*U14(L) - L52(L)*U24(L) $ - L53(L)*U34(L) L55(L) = 1. / (B(5,5,J,L) - L51(L)*U15(L) - L52(L)*U25(L) $ - L53(L)*U35(L) - L54(L)*U45(L)) 25 CONTINUE C C********** STEP 3. SOLVE FOR INTERMEDIATE VECTOR *************** C C A. CONSTRUCT RHS C IF(J.EQ.JS) GO TO 34 DO 33 M = 1,5 DO 33 L = LS,LE S(J,K,L,M) = S(J,K,L,M) - A(M,1,J,L)*S(J-1,K,L,1) $ - A(M,2,J,L)*S(J-1,K,L,2) - A(M,3,J,L)*S(J-1,K,L,3) $ - A(M,4,J,L)*S(J-1,K,L,4) - A(M,5,J,L)*S(J-1,K,L,5) 33 CONTINUE C C B. INTERMEDIATE VECTOR C 34 CONTINUE C C FWD SUBSTITUTION C DO 35 L = LS,LE D1 = S(J,K,L,1)*L11(L) D2 = (S(J,K,L,2) - L21(L)*D1) * L22(L) D3 = (S(J,K,L,3) - L31(L)*D1 - L32(L)*D2) * L33(L) D4 = (S(J,K,L,4) - L41(L)*D1 - L42(L)*D2 - L43(L)*D3) * L44(L) D5 = (S(J,K,L,5) - L51(L)*D1 - L52(L)*D2 - L53(L)*D3 $ - L54(L)*D4) * L55(L) C C BWD SUBSTITUTION C S(J,K,L,5) = D5 S(J,K,L,4) = D4 - U45(L)*D5 S(J,K,L,3) = D3 - U34(L)*S(J,K,L,4) - U35(L)*D5 S(J,K,L,2) = D2 - U23(L)*S(J,K,L,3) - U24(L)*S(J,K,L,4) $ - U25(L)*D5 S(J,K,L,1) = D1 - U12(L)*S(J,K,L,2) - U13(L)*S(J,K,L,3) $ - U14(L)*S(J,K,L,4) - U15(L)*D5 35 CONTINUE C C********** STEP 4. CONSTRUCT U(I) = L(I)**(-1)*C(I+1) ********** C********** BY COLUMNS AND STORE IN B ********** C IF(J.EQ.JE) GO TO 100 DO 40 N = 1,5 DO 40 L = LS,LE C C FWD SUBSTITUTION C C1 = C(1,N,J,L)*L11(L) C2 = (C(2,N,J,L) - L21(L)*C1) * L22(L) C3 = (C(3,N,J,L) - L31(L)*C1 - L32(L)*C2) * L33(L) C4 = (C(4,N,J,L) - L41(L)*C1 - L42(L)*C2 - L43(L)*C3) $ * L44(L) C5 = (C(5,N,J,L) - L51(L)*C1 - L52(L)*C2 - L53(L)*C3 $ - L54(L)*C4) * L55(L) C C BWD SUBSTITUTION C B(5,N,J,L) = C5 B(4,N,J,L) = C4 - U45(L)*C5 B(3,N,J,L) = C3 - U34(L)*B(4,N,J,L) - U35(L)*C5 B(2,N,J,L) = C2 - U23(L)*B(3,N,J,L) - U24(L)*B(4,N,J,L) $ - U25(L)*C5 B(1,N,J,L) = C1 - U12(L)*B(2,N,J,L) - U13(L)*B(3,N,J,L) $ - U14(L)*B(4,N,J,L) - U15(L)*C5 40 CONTINUE C C 100 CONTINUE C C PART 2. BACKWARD BLOCK SWEEP C JEM1 = JE - 1 C DO 200 J = JEM1,JS,-1 DO 200 M = 1,5 DO 200 L = LS,LE S(J,K,L,M) = S(J,K,L,M) - B(M,1,J,L)*S(J+1,K,L,1) $ - B(M,2,J,L)*S(J+1,K,L,2) - B(M,3,J,L)*S(J+1,K,L,3) $ - B(M,4,J,L)*S(J+1,K,L,4) - B(M,5,J,L)*S(J+1,K,L,5) 200 CONTINUE C RETURN END C SUBROUTINE GMTTST (ER, FP, TM) C PARAMETER (NW=100, NB=5, F7=78125., T30=1073741824.) COMPLEX WALL, ZCR, PROJ, ZI, Z1, ZZ COMMON /ARRAYS/ NWALL(NB), WALL(NW,NB), RMATRX(NW*NB,NW*NB), $ ZCR(NW,NB), PROJ(NW,NB), XMAX(NB) DATA IT/2/, ANS/-2.57754233214174/ C C INITIALIZATION C LW = 2 * NW * NB T2 = F7 / T30 DO 100 J = 1, NB NWALL(J) = NW 100 CONTINUE DO 110 J = 1, NB DO 110 I = 1, NW T1 = MOD (F7 * T2, 1.) T2 = MOD (F7 * T1, 1.) WALL(I,J) = CMPLX (T1, T2) 110 CONTINUE TM = CPTIME () C C TIMING TEST C DO 120 I = 1, IT CALL GMTRY 120 CONTINUE C TM = CPTIME () ER = ABS ((RMATRX(19,19) - ANS) / ANS) FP = IT * (120. * (NB*NW) ** 2 + 0.666 * (NB*NW) ** 3) C RETURN END C SUBROUTINE GMTRY C C COMPUTE SOLID-RELATED ARRAYS, GAUSS ELIMINATE THE MATRIX OF WALL C INFLUENCE COEFFICIENTS. C C 11/30/84 D H BAILEY REVISED CODE FOR NAS KERNEL TEST C PARAMETER (NW=100, NB=5) COMPLEX WALL, ZCR, PROJ, ZI, Z1, ZZ COMMON /ARRAYS/ NWALL(NB), WALL(NW,NB), RMATRX(NW*NB,NW*NB), $ ZCR(NW,NB), PROJ(NW,NB), XMAX(NB) C DATA ARCL /5./, PI /3.141592653589793/, PERIOD/3./ C C COMPUTE ARCLENGTH. C MATDIM = 0 ARCL = 0. YMIN = 1.E+50 YMAX = -1.E+50 PIDP = PI / PERIOD C DO 9 L = 1, NB MATDIM = MATDIM + NWALL(L) DO 9 K = 1,NWALL(L) ARCL = ARCL + ABS(WALL(K,L) - WALL(1+MOD(K,NWALL(L)), L)) 9 CONTINUE C C COMPUTE CORE RADIUS. C R0 = ARCL / (MATDIM*2.) SIGMA = R0 / 2. C C DEFINE CREATION POINTS. C DO 6 L = 1,NB DO 5 K = 1,NWALL(L) ZZ = WALL(1+MOD(K+NWALL(L)-2,NWALL(L)), L) & - WALL(1+MOD(K,NWALL(L)), L) ZCR(K,L) = WALL(K,L) + CMPLX(0., R0/ABS(ZZ)) * ZZ 5 CONTINUE C C CHECK THAT WALL AND CREATION POINTS ARE NOT CROSSED DUE TO C TOO SHARP A CONCAVE KINK OR AN ERROR IN DEFINING THE BODY. C ALSO FIND HIGHEST, LOWEST AND RIGHT-MOST POINT. C XMAX(L) = REAL(ZCR(1,L)) LS = 0 DO 6 K = 1,NWALL(L) YMIN = MIN (YMIN, AIMAG(ZCR(K,L))) YMAX = MAX (YMAX, AIMAG(ZCR(K,L))) XMAX(L) = MAX (XMAX(L), REAL(ZCR(K,L))) KP = 1 + MOD(K, NWALL(L)) IF (REAL((ZCR(KP,L) - ZCR(K,L)) * & CONJG(WALL(KP,L) - WALL(K,L))).GT.0.) THEN LS = L KS = K ENDIF 6 CONTINUE C C IF (LS .NE. 0) THEN C WRITE (6, 102) LS, KS C102 FORMAT(" ON BODY NUMBER ", I3, " YOU HAVE TOO SHARP A", C & " KINK NEAR POINT ", I4) C STOP C ENDIF C C THE "MAIN PERIOD" WILL BE BETWEEN YLIMIT AND YLIMIT + PERIOD. C YLIMIT = (YMIN - PERIOD + YMAX)/2 C C PROJECT CREATION POINTS INTO MAIN PERIOD. THIS IS TECHNICAL. C DO 1 L = 1,NB DO 1 K = 1,NWALL(L) PROJ(K,L) = ZCR(K,L) - CMPLX(0., PERIOD* & (INT(5. + (AIMAG(ZCR(K,L)) - YLIMIT) / PERIOD) - 5.)) 1 CONTINUE C C COMPUTE MATRIX. C SIG2 = (2. * PIDP * SIGMA) ** 2 I0 = 0 DO 2 L1 = 1,NB J0 = 0 DO 4 L2 = 1,NB KRON = 0 IF (L1 .EQ. L2) KRON = 1 DO 3 J = 1,NWALL(L2) RMATRX(I0+1,J0+J) = KRON Z1 = EXP ((WALL(1,L1) - ZCR(J,L2)) * PIDP) Z1 = Z1 - 1. / Z1 DUM = SIG2 + REAL(Z1)**2 + AIMAG(Z1)**2 DO 3 I = 2,NWALL(L1) ZI = EXP ((WALL(I,L1) - ZCR(J,L2)) * PIDP) ZZ = ZI - 1. / ZI RMATRX(I0+I,J0+J) = -0.25 / PI * LOG (DUM / & (SIG2 + REAL(ZZ) ** 2 + AIMAG(ZZ) ** 2)) 3 CONTINUE J0 = J0 + NWALL(L2) 4 CONTINUE I0 = I0 + NWALL(L1) 2 CONTINUE C C GAUSS ELIMINATION C DO 8 I = 1, MATDIM RMATRX(I,I) = 1. / RMATRX(I,I) DO 8 J = I+1, MATDIM RMATRX(J,I) = RMATRX(J,I) * RMATRX(I,I) DO 8 K = I+1, MATDIM RMATRX(J,K) = RMATRX(J,K) - RMATRX(J,I) * RMATRX(I,K) 8 CONTINUE C RETURN END C SUBROUTINE EMITST (ER, FP, TM) C C EMIT TEST SUBROUTINE C PARAMETER (NW=100, NB=5, NV=1000, NVM=1500, F7=78125., $ T30=1073741824.) COMPLEX Z, WALL, ZCR, REFPT, EXPWKL, EXPMWK, FORCE, & UUPSTR, DUM3, EXPZ, EXPMZ COMMON /ARRAYS/ NWALL(NB), WALL(NW,NB), RMATRX(NW*NB,NW*NB), $ ZCR(NW,NB), Z(NVM), GAMMA(NVM), REFPT(NB), RHS(NW*NB), $ FORCE(NB), RMOM(NB), CP(NW,NB), DPDS(NW,NB), EXPZ(NVM), $ EXPMZ(NVM), PSI(NW), PS(NVM) DATA IT/10/, ANS/6.0088546832072/ C T2 = F7 / T30 DO 100 J = 1, NB NWALL(J) = NW REFPT(J) = 0. FORCE(J) = 0. RMOM(J) = 0. DO 100 I = 1, NW T1 = MOD (F7 * T2, 1.) T2 = MOD (F7 * T1, 1.) WALL(I,J) = CMPLX (T1, T2) T1 = MOD (F7 * T2, 1.) T2 = MOD (F7 * T1, 1.) ZCR(I,J) = CMPLX (T1, T2) DPDS(I,J) = 0. 100 CONTINUE DO 110 J = 1, NW*NB RMATRX(J,J) = 1. DO 110 I = 1, J-1 T2 = MOD (F7 * T2, 1.) RMATRX(I,J) = .001 * T2 RMATRX(J,I) = .001 * T2 110 CONTINUE DO 120 I = 1, NVM T1 = MOD (F7 * T2, 1.) T2 = MOD (F7 * T1, 1.) Z(I) = CMPLX (T1, T2) T2 = MOD (F7 * T2, 1.) GAMMA(I) = T2 120 CONTINUE TM = CPTIME () C C TIMING TEST C DO 130 I = 1, IT CALL EMIT 130 CONTINUE C TM = CPTIME () ER = ABS ((RHS(19) - ANS) / ANS) FP = IT * (56.*NV + NB*NW * (97. + 44.*NV + 2.*NB*NW)) C RETURN END C SUBROUTINE EMIT C C EMIT NEW VORTICES TO SATISFY BOUNDARY CONDITION. C FINISH COMPUTING PRESSURE, FORCES, ETC. C C 11/28/84 D H BAILEY MODIFIED FOR NAS KERNEL TEST C PARAMETER (NW=100, NB=5, NVM=1500) COMPLEX Z, WALL, ZCR, REFPT, EXPWKL, EXPMWK, FORCE, & UUPSTR, DUM3, EXPZ, EXPMZ, ZZ COMMON /ARRAYS/ NWALL(NB), WALL(NW,NB), RMATRX(NW*NB,NW*NB), $ ZCR(NW,NB), Z(NVM), GAMMA(NVM), REFPT(NB), RHS(NW*NB), $ FORCE(NB), RMOM(NB), CP(NW,NB), DPDS(NW,NB), EXPZ(NVM), $ EXPMZ(NVM), PSI(NW), PS(NVM) C DATA PERIOD/3./, SIG2/3./, U0/4./, MATDIM/500/, DELT/1./, $ CHORD/5./, PI/3.141592653589793/, UUPSTR/(3., 4.)/ C C STORE EXP(Z(I)) AND EXP(-Z(I)) TO REDUCE WORK IN INNER LOOP 4. C NV = 1000 PIDP = PI / PERIOD C DO 2 I = 1, NV EXPZ(I) = EXP (Z(I) * PIDP) EXPMZ(I) = 1. / EXPZ(I) 2 CONTINUE C I0 = 0 CUPST = REAL(UUPSTR) ** 2 + AIMAG(UUPSTR) ** 2 C DO 5 L = 1, NB DO 6 K = 1, NWALL(L) EXPWKL = EXP (WALL(K,L) * PIDP) EXPMWK = 1. / EXPWKL SPS = 0. DO 4 I = 1, NV DUM3 = EXPZ(I) * EXPMWK - EXPWKL * EXPMZ(I) PS(I) = GAMMA(I) * LOG (REAL(DUM3) ** 2 + & AIMAG(DUM3) ** 2 + SIG2) SPS = SPS + PS(I) 4 CONTINUE PSI(K) = AIMAG(WALL(K,L) * CONJG (UUPSTR + CMPLX (0., U0))) & - SPS * 0.25 / PI 6 CONTINUE C C COMPUTE RIGHT-HAND SIDE. C DO 8 K = 1, NWALL(L) RHS(I0+K) = PSI(K) - PSI(1) 8 CONTINUE I0 = I0 + NWALL(L) 5 CONTINUE C C SOLVE SYSTEM C DO 10 I = 1, MATDIM DO 10 J = I+1, MATDIM RHS(J) = RHS(J) - RMATRX(J,I) * RHS(I) 10 CONTINUE DO 11 I = MATDIM, 1, -1 RHS(I) = RMATRX(I,I) * RHS(I) DO 11 J = 1, I-1 RHS(J) = RHS(J) - RMATRX(J,I) * RHS(I) 11 CONTINUE C C CREATE NEW VORTICES. C NOLLD = NV I0 = 0 DO 7 L = 1, NB DO 3 K = 1, NWALL(L) C C PUT THE NEW VORTEX AT THE END OF THE ARRAY. C NV = NV+1 Z(NV) = ZCR(K,L) GAMMA(NV) = RHS(I0+K) C C RECORD THE GAIN OF LINEAR AND ANGULAR MOMENTUM C FORCE(L) = FORCE(L) + GAMMA(NV) * Z(NV) RMOM(L) = RMOM(L) + GAMMA(NV) * (REAL (Z(NV) - REFPT(L)) ** 2 & + AIMAG (Z(NV) - REFPT(L)) ** 2) DPDS(K,L) = DPDS(K,L) - GAMMA(NV) 3 CONTINUE C C FILTER AND INTEGRATE PRESSURE GRADIENT TO GET PRESSURE C CP(1,L) = 0. CPM = -1E50 DO 9 K = 2, NWALL(L) CP(K,L) = CP(K-1,L) + (3. * (DPDS(K,L) + DPDS(K-1,L)) & + DPDS(1+MOD(K,NWALL(L)), L) & + DPDS(1+MOD(K+NWALL(L)-3, NWALL(L)), L)) & / (4. * DELT * CUPST) CPM = MAX (CPM, CP(K,L)) 9 CONTINUE C C NORMALIZE PRESSURE C DO 12 K = 1, NWALL(L) CP(K,L) = CP(K,L) - CPM 12 CONTINUE C C FINISH COMPUTING FORCE AND MOMENT, AS TIME RATE OF CHANGE OF LINEAR C AND ANGULAR MOMENTUM C FORCE(L) = FORCE(L) * CMPLX (0., 2. / (DELT * CHORD * CUPST)) RMOM(L) = RMOM(L) * 2. / (DELT * CHORD ** 2 * CUPST) C I0=I0+NWALL(L) 7 CONTINUE C RETURN END SUBROUTINE VPETST (ER, FP, TM) C C VPENTA TEST PROGRAM C PARAMETER (NJA=128, NJB=128, JL=1, JU=128, KL=1, KU=128, $ F7=78125., T30=1073741824.) COMMON /ARRAYS/ A(NJA,NJB), B(NJA,NJB), C(NJA,NJB), D(NJA,NJB), $ E(NJA,NJB), F(NJA,NJB,3), X(NJA,NJB), Y(NJA,NJB), FX(NJA,NJB,3) DATA IT/400/, ANS/-0.354649411858726/ C LF = NJA * NJB * 3 T = F7 / T30 DO 100 J = KL, KU DO 100 I = JL, JU T = MOD (F7 * T, 1.) A(I,J) = T T = MOD (F7 * T, 1.) B(I,J) = T T = MOD (F7 * T, 1.) C(I,J) = T T = MOD (F7 * T, 1.) D(I,J) = T T = MOD (F7 * T, 1.) E(I,J) = T DO 100 K = 1, 3 T = MOD (F7 * T, 1.) FX(I,J,K) = T 100 CONTINUE TM = CPTIME () C C TIMING TEST C DO 110 I = 1, IT CALL COPY (LF, FX, F) CALL VPENTA 110 CONTINUE C TM = CPTIME () ER = ABS ((F(19,19,1) - ANS) / ANS) FP = IT * KU * (40. * KU - 53.) C RETURN END C SUBROUTINE VPENTA C C ROUTINE TO INVERT 3 PENTADIAGONALS SIMULTANEOUSLY C C 12/05/84 D H BAILEY MODIFIED FOR NAS KERNEL TEST C PARAMETER (NJA=128, NJB=128, JL=1, JU=128, KL=1, KU=128) COMMON /ARRAYS/ A(NJA,NJB), B(NJA,NJB), C(NJA,NJB), D(NJA,NJB), $ E(NJA,NJB), F(NJA,NJB,3), X(NJA,NJB), Y(NJA,NJB), FX(NJA,NJB,3) C C ! START FORWARD GENERATION PROCESS AND SWEEP C J = JL DO 1 K = KL,KU RLD = C(J,K) RLDI = 1./RLD F(J,K,1) = F(J,K,1)*RLDI F(J,K,2) = F(J,K,2)*RLDI F(J,K,3) = F(J,K,3)*RLDI X(J,K) = D(J,K)*RLDI Y(J,K) = E(J,K)*RLDI 1 CONTINUE C J = JL+1 DO 2 K = KL,KU RLD1 = B(J,K) RLD = C(J,K) - RLD1*X(J-1,K) RLDI = 1./RLD F(J,K,1) = (F(J,K,1) - RLD1*F(J-1,K,1))*RLDI F(J,K,2) = (F(J,K,2) - RLD1*F(J-1,K,2))*RLDI F(J,K,3) = (F(J,K,3) - RLD1*F(J-1,K,3))*RLDI X(J,K) = (D(J,K) - RLD1*Y(J-1,K))*RLDI Y(J,K) = E(J,K)*RLDI 2 CONTINUE C DO 3 J = JL+2,JU-2 DO 11 K = KL,KU RLD2 = A(J,K) RLD1 = B(J,K) - RLD2*X(J-2,K) RLD = C(J,K) - (RLD2*Y(J-2,K) + RLD1*X(J-1,K)) RLDI = 1./RLD F(J,K,1) = (F(J,K,1) - RLD2*F(J-2,K,1) - RLD1*F(J-1,K,1))*RLDI F(J,K,2) = (F(J,K,2) - RLD2*F(J-2,K,2) - RLD1*F(J-1,K,2))*RLDI F(J,K,3) = (F(J,K,3) - RLD2*F(J-2,K,3) - RLD1*F(J-1,K,3))*RLDI X(J,K) = (D(J,K) - RLD1*Y(J-1,K))*RLDI Y(J,K) = E(J,K)*RLDI 11 CONTINUE 3 CONTINUE C J = JU-1 DO 12 K = KL,KU RLD2 = A(J,K) RLD1 = B(J,K) - RLD2*X(J-2,K) RLD = C(J,K) - (RLD2*Y(J-2,K) + RLD1*X(J-1,K)) RLDI = 1./RLD F(J,K,1) = (F(J,K,1) - RLD2*F(J-2,K,1) - RLD1*F(J-1,K,1))*RLDI F(J,K,2) = (F(J,K,2) - RLD2*F(J-2,K,2) - RLD1*F(J-1,K,2))*RLDI F(J,K,3) = (F(J,K,3) - RLD2*F(J-2,K,3) - RLD1*F(J-1,K,3))*RLDI X(J,K) = (D(J,K) - RLD1*Y(J-1,K))*RLDI 12 CONTINUE C J = JU DO 13 K = KL,KU RLD2 = A(J,K) RLD1 = B(J,K) - RLD2*X(J-2,K) RLD = C(J,K) - (RLD2*Y(J-2,K) + RLD1*X(J-1,K)) RLDI = 1./RLD F(J,K,1) = (F(J,K,1) - RLD2*F(J-2,K,1) - RLD1*F(J-1,K,1))*RLDI F(J,K,2) = (F(J,K,2) - RLD2*F(J-2,K,2) - RLD1*F(J-1,K,2))*RLDI F(J,K,3) = (F(J,K,3) - RLD2*F(J-2,K,3) - RLD1*F(J-1,K,3))*RLDI 13 CONTINUE C C ! BACK SWEEP SOLUTION C DO 14 K = KL,KU F(JU,K,1) = F(JU,K,1) F(JU,K,2) = F(JU,K,2) F(JU,K,3) = F(JU,K,3) F(JU-1,K,1) = F(JU-1,K,1) - X(JU-1,K)*F(JU,K,1) F(JU-1,K,2) = F(JU-1,K,2) - X(JU-1,K)*F(JU,K,2) F(JU-1,K,3) = F(JU-1,K,3) - X(JU-1,K)*F(JU,K,3) 14 CONTINUE C DO 4 J = 2,JU-JL JX = JU-J DO 15 K = KL,KU F(JX,K,1) = F(JX,K,1) - X(JX,K)*F(JX+1,K,1) - * Y(JX,K)*F(JX+2,K,1) F(JX,K,2) = F(JX,K,2) - X(JX,K)*F(JX+1,K,2) - * Y(JX,K)*F(JX+2,K,2) F(JX,K,3) = F(JX,K,3) - X(JX,K)*F(JX+1,K,3) - * Y(JX,K)*F(JX+2,K,3) 15 CONTINUE 4 CONTINUE C RETURN END