diff --git a/fish90/src/poisson.f90 b/fish90/src/poisson.f90 new file mode 100644 index 0000000000000000000000000000000000000000..57f8b42e6808d4819133a77e11a3325dec07701c --- /dev/null +++ b/fish90/src/poisson.f90 @@ -0,0 +1,608 @@ + 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