SUBROUTINE POISTGG(NPEROD,N,MPEROD,M,vecA,vecB,vecC,IDIMY,matY, & IERROR,W,IW) implicit none DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, & ONE = 1.0D0, TWO = 2.0D0, & FOUR = 4.0D0 !----------------------------------------------- ! D u m m y A r g u m e n t s !----------------------------------------------- INTEGER, INTENT(IN) :: NPEROD, N, MPEROD, M, IDIMY,IW INTEGER, INTENT(OUT) :: IERROR DOUBLE PRECISION, DIMENSION(M), INTENT(IN) :: vecA,vecB,vecC DOUBLE PRECISION, DIMENSION(IDIMY,N),INTENT(INOUT) :: matY DOUBLE PRECISION, DIMENSION(IW) :: W !----------------------------------------------- ! L o c a l V a r i a b l e s !----------------------------------------------- INTEGER :: IWBA, IWBB, IWBC, IWB2, IWB3, IWW1, IWW2, IWW3, & IWD,IWTCOS, IWP, I, K, J, NP, MP, IPSTOR, & IREV, MH, MHM1, MODD, NBY2, MSKIP DOUBLE PRECISION :: A1 !----------------------------------------------- IERROR = 0 IWBA = M + 1 IWBB = IWBA + M IWBC = IWBB + M IWB2 = IWBC + M IWB3 = IWB2 + M IWW1 = IWB3 + M IWW2 = IWW1 + M IWW3 = IWW2 + M IWD = IWW3 + M IWTCOS = IWD + M IWP = IWTCOS + 4*N DO I = 1, M K = IWBA + I - 1 ! write(*,*) '0: i,k:',i,k ! write(*,*) '0: vecA(I):',vecA(I) W(K) = -vecA(I) ! write(*,*) 'A: i,k,w(k):',i,k,w(k) K = IWBC + I - 1 W(K) = -vecC(I) ! write(*,*) 'B: i,k,w(k):',i,k,w(k) K = IWBB + I - 1 W(K) = TWO - vecB(I) ! write(*,*) 'C: i,k,w(k):',i,k,w(k) matY(I,:N) = -matY(I,:N) ! write(*,*) 'D: Y(i,:N):',Y(I,:N) END DO NP = NPEROD MP = MPEROD + 1 GO TO (110,107) MP 107 CONTINUE GO TO (108,108,108,119) NPEROD 108 CONTINUE CALL POSTG2 (NP, N, M, W(IWBA:IWBA+M-1), W(IWBB:IWBB+M-1), & W(IWBC:IWBC+M-1), IDIMY, matY, W(1:M), W(IWB2:IWB2+M-1), & W(IWB3:IWB3+M-1), W(IWW1:IWW1+M-1), W(IWW2:IWW2+M-1), & W(IWW3:IWW3+M-1), W(IWD:IWD+M-1), W(IWTCOS:IWTCOS+4*N-1), & W(IWP:),IW-IWP) IPSTOR = W(IWW1) IREV = 2 IF (NPEROD == 4) GO TO 120 109 CONTINUE GO TO (123,129) MP 110 CONTINUE MH = (M + 1)/2 MHM1 = MH - 1 MODD = 1 IF (MH*2 == M) MODD = 2 DO J = 1, N DO I = 1, MHM1 W(I) = matY(MH-I,J) - matY(I+MH,J) W(I+MH) = matY(MH-I,J) + matY(I+MH,J) END DO W(MH) = TWO*matY(MH,J) GO TO (113,112) MODD 112 CONTINUE W(M) = TWO*matY(M,J) 113 CONTINUE matY(:M,J) = W(:M) END DO K = IWBC + MHM1 - 1 I = IWBA + MHM1 W(K) = 0. W(I) = 0. W(K+1) = TWO*W(K+1) SELECT CASE (MODD) CASE DEFAULT K = IWBB + MHM1 - 1 W(K) = W(K) - W(I-1) W(IWBC-1) = W(IWBC-1) + W(IWBB-1) CASE (2) W(IWBB-1) = W(K+1) END SELECT GO TO 107 119 CONTINUE IREV = 1 NBY2 = N/2 NP = 2 120 CONTINUE DO J = 1, NBY2 MSKIP = N + 1 - J DO I = 1, M A1 = matY(I,J) matY(I,J) = matY(I,MSKIP) matY(I,MSKIP) = A1 END DO END DO GO TO (108,109) IREV 123 CONTINUE DO J = 1, N W(MH-1:MH-MHM1:(-1)) = HALF*(matY(MH+1:MHM1+MH,J)+matY(:MHM1,J)) W(MH+1:MHM1+MH) = HALF*(matY(MH+1:MHM1+MH,J)-matY(:MHM1,J)) W(MH) = HALF*matY(MH,J) GO TO (126,125) MODD 125 CONTINUE W(M) = HALF*matY(M,J) 126 CONTINUE matY(:M,J) = W(:M) END DO 129 CONTINUE W(1) = IPSTOR + IWP - 1 ! END SUBROUTINE POISTGG SUBROUTINE POSTG2(NPEROD, N, M, vecA, vecBB, vecC, IDIMQ, matQ, & vecB, vecB2, vecB3, vecW,vecW2, vecW3, vecD, & TCOS, vecP,IvecP) implicit none DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, & ONE = 1.0D0, TWO = 2.0D0, & FOUR = 4.0D0 !----------------------------------------------- ! D u m m y A r g u m e n t s !----------------------------------------------- INTEGER , INTENT(IN) :: NPEROD,N,M,IDIMQ,IvecP DOUBLE PRECISION, DIMENSION(M),INTENT(IN) :: vecA,vecBB,vecC DOUBLE PRECISION, DIMENSION(M),INTENT(OUT) :: vecB DOUBLE PRECISION, DIMENSION(IDIMQ,N),INTENT(INOUT) :: matQ DOUBLE PRECISION, DIMENSION(M),INTENT(INOUT) :: vecB2,vecB3,vecD, & vecW2,vecW3, vecW DOUBLE PRECISION, DIMENSION(IvecP),INTENT(INOUT) :: vecP DOUBLE PRECISION, DIMENSION(4*N),INTENT(INOUT) :: TCOS !----------------------------------------------- ! L o c a l V a r i a b l e s !----------------------------------------------- INTEGER, DIMENSION(4) :: K INTEGER :: K1, K2, K3, K4, NP, MR, IP, IPSTOR, I2R, JR, NR, NLAST, & KR, LR, NROD, JSTART, JSTOP, I2RBY2, & J, IJUMP, JP1, JP2, JP3, JM1, JM2, JM3, I, NRODPR, II, & NLASTP, JSTEP DOUBLE PRECISION :: FNUM, FNUM2, FI, T !----------------------------------------------- ! ! SUBROUTINE TO SOLVE POISSON'S EQUATION ON A STAGGERED GRID. ! EQUIVALENCE (K(1), K1), (K(2), K2), (K(3), K3), (K(4), K4) NP = NPEROD FNUM = HALF*DBLE(NP/3) FNUM2 = HALF*DBLE(NP/2) MR = M IP = -MR IPSTOR = 0 I2R = 1 JR = 2 NR = N NLAST = N KR = 1 LR = 0 IF (NR > 3) THEN 101 CONTINUE JR = 2*I2R NROD = 1 IF ((NR/2)*2 == NR) NROD = 0 JSTART = 1 JSTOP = NLAST - JR IF (NROD == 0) JSTOP = JSTOP - I2R I2RBY2 = I2R/2 IF (JSTOP < JSTART) THEN J = JR ELSE IJUMP = 1 DO J = JSTART, JSTOP, JR JP1 = J + I2RBY2 JP2 = J + I2R JP3 = JP2 + I2RBY2 JM1 = J - I2RBY2 JM2 = J - I2R JM3 = JM2 - I2RBY2 IF (J == 1) THEN CALL COSGEN (I2R, 1, FNUM, HALF, TCOS, size(TCOS) ) IF (I2R == 1) THEN vecB(:MR) = matQ(:MR,1) matQ(:MR,1) = matQ(:MR,2) GO TO 112 ENDIF vecB(:MR) = matQ(:MR,1) + & HALF*(matQ(:MR,JP2)-matQ(:MR,JP1)-matQ(:MR,JP3)) matQ(:MR,1) = matQ(:MR,JP2) + matQ(:MR,1) - matQ(:MR,JP1) GO TO 112 ENDIF GO TO (107,108) IJUMP 107 CONTINUE IJUMP = 2 CALL COSGEN (I2R, 1, HALF, ZERO, TCOS, size(TCOS) ) 108 CONTINUE IF (I2R == 1) THEN vecB(:MR) = TWO*matQ(:MR,J) matQ(:MR,J) = matQ(:MR,JM2) + matQ(:MR,JP2) ELSE DO I = 1, MR FI = matQ(I,J) matQ(I,J)=matQ(I,J)-matQ(I,JM1)-matQ(I,JP1)+ & matQ(I,JM2)+matQ(I,JP2) vecB(I) =FI + matQ(I,J) - matQ(I,JM3) - matQ(I,JP3) END DO ENDIF 112 CONTINUE CALL TRIX (I2R, 0, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) matQ(:MR,J) = matQ(:MR,J) + vecB(:MR) ! ! END OF REDUCTION FOR REGULAR UNKNOWNS. ! END DO ! ! BEGIN SPECIAL REDUCTION FOR LAST UNKNOWN. ! J = JSTOP + JR ENDIF NLAST = J JM1 = J - I2RBY2 JM2 = J - I2R JM3 = JM2 - I2RBY2 IF (NROD /= 0) THEN ! ! ODD NUMBER OF UNKNOWNS ! IF (I2R == 1) THEN vecB(:MR) = matQ(:MR,J) matQ(:MR,J) = matQ(:MR,JM2) ELSE vecB(:MR)=matQ(:MR,J)+ & HALF*(matQ(:MR,JM2)-matQ(:MR,JM1)-matQ(:MR,JM3)) IF (NRODPR == 0) THEN matQ(:MR,J) = matQ(:MR,JM2) + vecP(IP+1:MR+IP) IP = IP - MR ELSE matQ(:MR,J) = matQ(:MR,J) - matQ(:MR,JM1) + matQ(:MR,JM2) ENDIF IF (LR /= 0) CALL COSGEN (LR, 1, FNUM2, HALF, TCOS(KR+1), & size(tcos)-KR) ENDIF CALL COSGEN (KR, 1, FNUM2, HALF, TCOS, size(tcos)) CALL TRIX (KR, LR, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) matQ(:MR,J) = matQ(:MR,J) + vecB(:MR) KR = KR + I2R ELSE JP1 = J + I2RBY2 JP2 = J + I2R IF (I2R == 1) THEN vecB(:MR) = matQ(:MR,J) TCOS(1) = 0. CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) IP = 0 IPSTOR = MR vecP(:MR) = vecB(:MR) vecB(:MR) = vecB(:MR) + matQ(:MR,N) TCOS(1) = -ONE + TWO*DBLE(NP/2) TCOS(2) = 0. CALL TRIX (1, 1, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) matQ(:MR,J) = matQ(:MR,JM2) + vecP(:MR) + vecB(:MR) ELSE vecB(:MR)=matQ(:MR,J)+ & HALF*(matQ(:MR,JM2)-matQ(:MR,JM1)-matQ(:MR,JM3)) IF (NRODPR == 0) THEN vecB(:MR) = vecB(:MR) + vecP(IP+1:MR+IP) ELSE vecB(:MR) = vecB(:MR) + matQ(:MR,JP2) - matQ(:MR,JP1) ENDIF CALL COSGEN (I2R, 1, HALF, ZERO, TCOS, size(tcos)) CALL TRIX (I2R, 0, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) IP = IP + MR IPSTOR = MAX0(IPSTOR,IP + MR) vecP(IP+1:MR+IP) = vecB(:MR) + & HALF*(matQ(:MR,J)-matQ(:MR,JM1)-matQ(:MR,JP1)) vecB(:MR) = vecP(IP+1:MR+IP) + matQ(:MR,JP2) IF (LR /= 0) THEN CALL COSGEN (LR, 1, FNUM2, HALF, TCOS(I2R+1), size(tcos)-I2R) CALL MERGE (TCOS, 0, I2R, I2R, LR, KR, size(TCOS)) ELSE DO I = 1, I2R II = KR + I TCOS(II) = TCOS(I) END DO ENDIF CALL COSGEN (KR, 1, FNUM2, HALF, TCOS, size(TCOS)) CALL TRIX (KR, KR, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) matQ(:MR,J) = matQ(:MR,JM2) + vecP(IP+1:MR+IP) + vecB(:MR) ENDIF LR = KR KR = KR + JR ENDIF NR = (NLAST - 1)/JR + 1 IF (NR <= 3) GO TO 142 I2R = JR NRODPR = NROD GO TO 101 ENDIF ! IF (NR > 3) 142 CONTINUE J = 1 + JR JM1 = J - I2R JP1 = J + I2R JM2 = NLAST - I2R IF (NR /= 2) THEN IF (LR == 0) THEN IF (N == 3) THEN ! ! CASE N = 3. ! GO TO (143,148,143) NP 143 CONTINUE vecB(:MR) = matQ(:MR,2) vecB2(:MR) = matQ(:MR,1) + matQ(:MR,3) vecB3(:MR) = 0. SELECT CASE (NP) CASE DEFAULT TCOS(1) = -ONE TCOS(2) = ONE K1 = 1 CASE (1:2) TCOS(1) = -TWO TCOS(2) = ONE TCOS(3) = -ONE K1 = 2 END SELECT K2 = 1 K3 = 0 K4 = 0 GO TO 150 148 CONTINUE vecB(:MR) = matQ(:MR,2) vecB2(:MR) = matQ(:MR,3) vecB3(:MR) = matQ(:MR,1) CALL COSGEN (3, 1, HALF, ZERO, TCOS, size(TCOS)) TCOS(4) = -ONE TCOS(5) = ONE TCOS(6) = -ONE TCOS(7) = ONE K1 = 3 K2 = 2 K3 = 1 K4 = 1 150 CONTINUE CALL TRI3(MR,vecA,vecBB,vecC,K,vecB,vecB2,vecB3,TCOS, & size(TCOS),vecD,vecW,vecW2,vecW3) vecB(:MR) = vecB(:MR) + vecB2(:MR) + vecB3(:MR) GO TO (153,153,152) NP 152 CONTINUE TCOS(1) = TWO CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) 153 CONTINUE matQ(:MR,2) = vecB(:MR) vecB(:MR) = matQ(:MR,1) + vecB(:MR) TCOS(1) = -ONE + FOUR*FNUM CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) matQ(:MR,1) = vecB(:MR) JR = 1 I2R = 0 GO TO 188 ENDIF ! ! CASE N = 2**P+1 ! vecB(:MR)=matQ(:MR,J)+matQ(:MR,1)-matQ(:MR,JM1)+ & matQ(:MR,NLAST)-matQ(:MR,JM2) GO TO (158,160,158) NP 158 CONTINUE vecB2(:MR) = matQ(:MR,1) + matQ(:MR,NLAST) + & matQ(:MR,J) - matQ(:MR,JM1) -matQ(:MR,JP1) vecB3(:MR) = 0. K1 = NLAST - 1 K2 = NLAST + JR - 1 CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(NLAST), size(TCOS)-NLAST+1) TCOS(K2) = TWO*DBLE(NP - 2) CALL COSGEN (JR, 1, HALF - FNUM, HALF, TCOS(K2+1),size(TCOS)-K2) K3 = (3 - NP)/2 CALL MERGE (TCOS, K1, JR - K3, K2 - K3, JR + K3, 0, size(TCOS)) K1 = K1 - 1 + K3 CALL COSGEN (JR, 1, FNUM, HALF, TCOS(K1+1),size(TCOS)-K1) K2 = JR K3 = 0 K4 = 0 GO TO 162 160 CONTINUE DO I = 1, MR FI = HALF*(matQ(I,J)-matQ(I,JM1)-matQ(I,JP1)) vecB2(I) = matQ(I,1) + FI vecB3(I) = matQ(I,NLAST) + FI END DO K1 = NLAST + JR - 1 K2 = K1 + JR - 1 CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(K1+1), size(TCOS)-K1) CALL COSGEN (NLAST, 1, HALF, ZERO, TCOS(K2+1), size(TCOS)-K2) CALL MERGE (TCOS, K1, JR - 1, K2, NLAST, 0, size(TCOS)) K3 = K1 + NLAST - 1 K4 = K3 + JR CALL COSGEN (JR, 1, HALF, HALF, TCOS(K3+1),SIZE(TCOS)-K3) CALL COSGEN (JR, 1, ZERO, HALF, TCOS(K4+1),SIZE(TCOS)-K4) CALL MERGE (TCOS, K3, JR, K4, JR, K1, size(TCOS) ) K2 = NLAST - 1 K3 = JR K4 = JR 162 CONTINUE CALL TRI3 (MR, vecA, vecBB, vecC, K, vecB, vecB2, vecB3, & TCOS, size(TCOS), vecD, vecW, vecW2, vecW3) vecB(:MR) = vecB(:MR) + vecB2(:MR) + vecB3(:MR) IF (NP == 3) THEN TCOS(1) = TWO CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) ENDIF matQ(:MR,J) = vecB(:MR) + & HALF*(matQ(:MR,J)-matQ(:MR,JM1)-matQ(:MR,JP1)) vecB(:MR) = matQ(:MR,J) + matQ(:MR,1) CALL COSGEN (JR, 1, FNUM, HALF, TCOS,size(TCOS)) CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) matQ(:MR,1) = matQ(:MR,1) - matQ(:MR,JM1) + vecB(:MR) GO TO 188 ENDIF ! ! CASE OF GENERAL N WITH NR = 3 . ! vecB(:MR) = matQ(:MR,1) - matQ(:MR,JM1) + matQ(:MR,J) IF (NROD == 0) THEN vecB(:MR) = vecB(:MR) + vecP(IP+1:MR+IP) ELSE vecB(:MR) = vecB(:MR) + matQ(:MR,NLAST) - matQ(:MR,JM2) ENDIF DO I = 1, MR T = HALF*(matQ(I,J)-matQ(I,JM1)-matQ(I,JP1)) matQ(I,J) = T vecB2(I) = matQ(I,NLAST) + T vecB3(I) = matQ(I,1) + T END DO K1 = KR + 2*JR CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(K1+1),size(TCOS)-K1) K2 = K1 + JR TCOS(K2) = TWO*DBLE(NP - 2) K4 = (NP - 1)*(3 - NP) K3 = K2 + 1 - K4 CALL COSGEN(KR+JR+K4,1,HALF*DBLE(K4),ONE-DBLE(K4),TCOS(K3), & size(TCOS)-K3+1 ) K4 = 1 - NP/3 CALL MERGE (TCOS, K1, JR - K4, K2 - K4, KR + JR + K4, 0, size(TCOS) ) IF (NP == 3) K1 = K1 - 1 K2 = KR + JR K4 = K1 + K2 CALL COSGEN (KR, 1, FNUM2, HALF, TCOS(K4+1),size(TCOS)-K4) K3 = K4 + KR CALL COSGEN (JR, 1, FNUM, HALF, TCOS(K3+1),size(TCOS)-K3) CALL MERGE (TCOS, K4, KR, K3, JR, K1, size(TCOS) ) K4 = K3 + JR CALL COSGEN (LR, 1, FNUM2, HALF, TCOS(K4+1),size(TCOS)-K4) CALL MERGE (TCOS, K3, JR, K4, LR, K1 + K2, size(TCOS) ) CALL COSGEN (KR, 1, FNUM2, HALF, TCOS(K3+1),size(TCOS)-K3) K3 = KR K4 = KR CALL TRI3 (MR, vecA, vecBB, vecC, K, vecB, vecB2, vecB3, & TCOS, SIZE(TCOS), vecD, vecW, vecW2, vecW3) vecB(:MR) = vecB(:MR) + vecB2(:MR) + vecB3(:MR) IF (NP == 3) THEN TCOS(1) = TWO CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) ENDIF matQ(:MR,J) = matQ(:MR,J) + vecB(:MR) vecB(:MR) = matQ(:MR,1) + matQ(:MR,J) CALL COSGEN (JR, 1, FNUM, HALF, TCOS,size(TCOS)) CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) IF (JR == 1) THEN matQ(:MR,1) = vecB(:MR) GO TO 188 ENDIF matQ(:MR,1) = matQ(:MR,1) - matQ(:MR,JM1) + vecB(:MR) GO TO 188 ENDIF vecB3(:MR) = 0. vecB (:MR) = matQ(:MR,1) + vecP(IP+1:MR+IP) matQ (:MR,1) = matQ(:MR,1) - matQ(:MR,JM1) vecB2(:MR) = matQ(:MR,1) + matQ(:MR,NLAST) K1 = KR + JR K2 = K1 + JR CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(K1+1),SIZE(TCOS)-K1) GO TO (182,183,182) NP 182 CONTINUE TCOS(K2) = TWO*DBLE(NP - 2) CALL COSGEN (KR, 1, ZERO, ONE, TCOS(K2+1),SIZE(TCOS)-K2) GO TO 184 183 CONTINUE CALL COSGEN (KR + 1, 1, HALF, ZERO, TCOS(K2),SIZE(TCOS)-K2) 184 CONTINUE K4 = 1 - NP/3 CALL MERGE (TCOS, K1, JR - K4, K2 - K4, KR + K4, 0, size(TCOS) ) IF (NP == 3) K1 = K1 - 1 K2 = KR CALL COSGEN (KR, 1, FNUM2, HALF, TCOS(K1+1),SIZE(TCOS)-K1) K4 = K1 + KR CALL COSGEN (LR, 1, FNUM2, HALF, TCOS(K4+1),SIZE(TCOS)-K4) K3 = LR K4 = 0 CALL TRI3 (MR, vecA, vecBB, vecC, K, vecB, vecB2, vecB3, & TCOS, SIZE(TCOS), vecD, vecW, vecW2, vecW3) vecB(:MR) = vecB(:MR) + vecB2(:MR) IF (NP == 3) THEN TCOS(1) = TWO CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) ENDIF matQ(:MR,1) = matQ(:MR,1) + vecB(:MR) 188 CONTINUE J = NLAST - JR vecB(:MR) = matQ(:MR,NLAST) + matQ(:MR,J) JM2 = NLAST - I2R IF (JR == 1) THEN matQ(:MR,NLAST) = 0. ELSE IF (NROD == 0) THEN matQ(:MR,NLAST) = vecP(IP+1:MR+IP) IP = IP - MR ELSE matQ(:MR,NLAST) = matQ(:MR,NLAST) - matQ(:MR,JM2) ENDIF ENDIF CALL COSGEN (KR, 1, FNUM2, HALF, TCOS,SIZE(TCOS)) CALL COSGEN (LR, 1, FNUM2, HALF, TCOS(KR+1),SIZE(TCOS)-KR) CALL TRIX (KR, LR, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) matQ(:MR,NLAST) = matQ(:MR,NLAST) + vecB(:MR) NLASTP = NLAST 197 CONTINUE JSTEP = JR JR = I2R I2R = I2R/2 IF (JR == 0) GO TO 210 JSTART = 1 + JR KR = KR - JR IF (NLAST + JR <= N) THEN KR = KR - JR NLAST = NLAST + JR JSTOP = NLAST - JSTEP ELSE JSTOP = NLAST - JR ENDIF LR = KR - JR CALL COSGEN (JR, 1, HALF, ZERO, TCOS,SIZE(TCOS)) DO J = JSTART, JSTOP, JSTEP JM2 = J - JR JP2 = J + JR IF (J == JR) THEN vecB(:MR) = matQ(:MR,J) + matQ(:MR,JP2) ELSE vecB(:MR) = matQ(:MR,J) + matQ(:MR,JM2) + matQ(:MR,JP2) ENDIF IF (JR == 1) THEN matQ(:MR,J) = 0. ELSE JM1 = J - I2R JP1 = J + I2R matQ(:MR,J) = HALF*(matQ(:MR,J)-matQ(:MR,JM1)-matQ(:MR,JP1)) ENDIF CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS, & size(TCOS), vecD, vecW) matQ(:MR,J) = matQ(:MR,J) + vecB(:MR) END DO NROD = 1 IF (NLAST + I2R <= N) NROD = 0 IF (NLASTP /= NLAST) GO TO 188 GO TO 197 210 CONTINUE vecW(1) = IPSTOR ! ! REVISION HISTORY--- ! ! SEPTEMBER 1973 VERSION 1 ! APRIL 1976 VERSION 2 ! JANUARY 1978 VERSION 3 ! DECEMBER 1979 VERSION 3.1 ! FEBRUARY 1985 DOCUMENTATION UPGRADE ! NOVEMBER 1988 VERSION 3.2, FORTRAN 77 CHANGES ! June 2004 Version 5.0, Fortran 90 Changes !----------------------------------------------------------------------- END SUBROUTINE POSTG2