! ! file genbun.f ! ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ! . . ! . copyright (c) 2004 by UCAR . ! . . ! . UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH . ! . . ! . all rights reserved . ! . . ! . . ! . FISHPACK version 5.0 . ! . . ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ! ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! * * ! * F I S H P A C K * ! * * ! * * ! * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF * ! * * ! * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS * ! * * ! * (Version 5.0 , JUNE 2004) * ! * * ! * BY * ! * * ! * JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET * ! * * ! * OF * ! * * ! * THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH * ! * * ! * BOULDER, COLORADO (80307) U.S.A. * ! * * ! * WHICH IS SPONSORED BY * ! * * ! * THE NATIONAL SCIENCE FOUNDATION * ! * * ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ! SUBROUTINE GENBUN (NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y,IERROR) ! ! ! DIMENSION OF A(M),B(M),C(M),Y(IDIMY,N) ! ARGUMENTS ! ! LATEST REVISION JUNE 2004 ! ! PURPOSE THE NAME OF THIS PACKAGE IS A MNEMONIC FOR THE ! GENERALIZED BUNEMAN ALGORITHM. ! ! IT SOLVES THE REAL LINEAR SYSTEM OF EQUATIONS ! ! A(I)*X(I-1,J) + B(I)*X(I,J) + C(I)*X(I+1,J) ! + X(I,J-1) - 2.*X(I,J) + X(I,J+1) = Y(I,J) ! ! FOR I = 1,2,...,M AND J = 1,2,...,N. ! ! INDICES I+1 AND I-1 ARE EVALUATED MODULO M, ! I.E., X(0,J) = X(M,J) AND X(M+1,J) = X(1,J), ! AND X(I,0) MAY EQUAL 0, X(I,2), OR X(I,N), ! AND X(I,N+1) MAY EQUAL 0, X(I,N-1), OR X(I,1) ! DEPENDING ON AN INPUT PARAMETER. ! ! USAGE CALL GENBUN (NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y, ! IERROR) ! ! ARGUMENTS ! ! ON INPUT NPEROD ! ! INDICATES THE VALUES THAT X(I,0) AND ! X(I,N+1) ARE ASSUMED TO HAVE. ! ! = 0 IF X(I,0) = X(I,N) AND X(I,N+1) = ! X(I,1). ! = 1 IF X(I,0) = X(I,N+1) = 0 . ! = 2 IF X(I,0) = 0 AND X(I,N+1) = X(I,N-1). ! = 3 IF X(I,0) = X(I,2) AND X(I,N+1) = ! X(I,N-1). ! = 4 IF X(I,0) = X(I,2) AND X(I,N+1) = 0. ! ! N ! THE NUMBER OF UNKNOWNS IN THE J-DIRECTION. ! N MUST BE GREATER THAN 2. ! ! MPEROD ! = 0 IF A(1) AND C(M) ARE NOT ZERO ! = 1 IF A(1) = C(M) = 0 ! ! M ! THE NUMBER OF UNKNOWNS IN THE I-DIRECTION. ! N MUST BE GREATER THAN 2. ! ! A,B,C ! ONE-DIMENSIONAL ARRAYS OF LENGTH M THAT ! SPECIFY THE COEFFICIENTS IN THE LINEAR ! EQUATIONS GIVEN ABOVE. IF MPEROD = 0 ! THE ARRAY ELEMENTS MUST NOT DEPEND UPON ! THE INDEX I, BUT MUST BE CONSTANT. ! SPECIFICALLY, THE SUBROUTINE CHECKS THE ! FOLLOWING CONDITION . ! ! A(I) = C(1) ! C(I) = C(1) ! B(I) = B(1) ! ! FOR I=1,2,...,M. ! ! IDIMY ! THE ROW (OR FIRST) DIMENSION OF THE ! TWO-DIMENSIONAL ARRAY Y AS IT APPEARS ! IN THE PROGRAM CALLING GENBUN. ! THIS PARAMETER IS USED TO SPECIFY THE ! VARIABLE DIMENSION OF Y. ! IDIMY MUST BE AT LEAST M. ! ! Y ! A TWO-DIMENSIONAL COMPLEX ARRAY THAT ! SPECIFIES THE VALUES OF THE RIGHT SIDE ! OF THE LINEAR SYSTEM OF EQUATIONS GIVEN ! ABOVE. ! Y MUST BE DIMENSIONED AT LEAST M*N. ! ! ! ON OUTPUT Y ! ! CONTAINS THE SOLUTION X. ! ! IERROR ! AN ERROR FLAG WHICH INDICATES INVALID ! INPUT PARAMETERS EXCEPT FOR NUMBER ! ZERO, A SOLUTION IS NOT ATTEMPTED. ! ! = 0 NO ERROR. ! = 1 M .LE. 2 . ! = 2 N .LE. 2 ! = 3 IDIMY .LT. M ! = 4 NPEROD .LT. 0 OR NPEROD .GT. 4 ! = 5 MPEROD .LT. 0 OR MPEROD .GT. 1 ! = 6 A(I) .NE. C(1) OR C(I) .NE. C(1) OR ! B(I) .NE. B(1) FOR ! SOME I=1,2,...,M. ! = 7 A(1) .NE. 0 OR C(M) .NE. 0 AND ! MPEROD = 1 ! = 20 If the dynamic allocation of real and ! complex work space required for solution ! fails (for example if N,M are too large ! for your computer) ! ! ! SPECIAL CONDITONS NONE ! ! I/O NONE ! ! PRECISION SINGLE ! ! REQUIRED FILES comf.f,gnbnaux.f,fish.f ! FILES ! ! LANGUAGE FORTRAN 90 ! ! HISTORY WRITTEN IN 1979 BY ROLAND SWEET OF NCAR'S ! SCIENTIFIC COMPUTING DIVISION. MADE AVAILABLE ! ON NCAR'S PUBLIC LIBRARIES IN JANUARY, 1980. ! Revised in June 2004 by John Adams using ! Fortran 90 dynamically allocated work space. ! ! ALGORITHM THE LINEAR SYSTEM IS SOLVED BY A CYCLIC ! REDUCTION ALGORITHM DESCRIBED IN THE ! REFERENCE. ! ! PORTABILITY FORTRAN 90 -- ! THE MACHINE DEPENDENT CONSTANT PI IS ! DEFINED IN FUNCTION PIMACH. ! ! REFERENCES SWEET, R., "A CYCLIC REDUCTION ALGORITHM FOR ! SOLVING BLOCK TRIDIAGONAL SYSTEMS OF ARBITRARY ! DIMENSIONS," SIAM J. ON NUMER. ANAL., 14 (1977) ! PP. 706-720. ! ! ACCURACY THIS TEST WAS PERFORMED ON a platform with ! 64 bit floating point arithmetic. ! A UNIFORM RANDOM NUMBER GENERATOR WAS USED ! TO CREATE A SOLUTION ARRAY X FOR THE SYSTEM ! GIVEN IN THE 'PURPOSE' DESCRIPTION ABOVE ! WITH ! A(I) = C(I) = -0.5*B(I) = 1, I=1,2,...,M ! ! AND, WHEN MPEROD = 1 ! ! A(1) = C(M) = 0 ! A(M) = C(1) = 2. ! ! THE SOLUTION X WAS SUBSTITUTED INTO THE ! GIVEN SYSTEM AND, USING DOUBLE PRECISION ! A RIGHT SIDE Y WAS COMPUTED. ! USING THIS ARRAY Y, SUBROUTINE GENBUN ! WAS CALLED TO PRODUCE APPROXIMATE ! SOLUTION Z. THEN RELATIVE ERROR ! E = MAX(ABS(Z(I,J)-X(I,J)))/ ! MAX(ABS(X(I,J))) ! WAS COMPUTED, WHERE THE TWO MAXIMA ARE TAKEN ! OVER I=1,2,...,M AND J=1,...,N. ! ! THE VALUE OF E IS GIVEN IN THE TABLE ! BELOW FOR SOME TYPICAL VALUES OF M AND N. ! ! M (=N) MPEROD NPEROD E ! ------ ------ ------ ------ ! ! 31 0 0 6.E-14 ! 31 1 1 4.E-13 ! 31 1 3 3.E-13 ! 32 0 0 9.E-14 ! 32 1 1 3.E-13 ! 32 1 3 1.E-13 ! 33 0 0 9.E-14 ! 33 1 1 4.E-13 ! 33 1 3 1.E-13 ! 63 0 0 1.E-13 ! 63 1 1 1.E-12 ! 63 1 3 2.E-13 ! 64 0 0 1.E-13 ! 64 1 1 1.E-12 ! 64 1 3 6.E-13 ! 65 0 0 2.E-13 ! 65 1 1 1.E-12 ! 65 1 3 4.E-13 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE GENBUN(NPEROD, N, MPEROD, M, A, B, C, IDIMY, Y, IERROR) USE fish implicit none DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, & ONE = 1.0D0, TWO = 2.0D0, & FOUR = 4.0D0 TYPE(fishworkspace) :: w !----------------------------------------------- ! D u m m y A r g u m e n t s !----------------------------------------------- INTEGER, INTENT(IN) :: NPEROD, N, MPEROD, M, IDIMY INTEGER, INTENT(OUT) :: IERROR DOUBLE PRECISION, DIMENSION(:) :: A, B, C DOUBLE PRECISION, INTENT(INOUT) :: Y(IDIMY,*) !----------------------------------------------- ! L o c a l V a r i a b l e s !----------------------------------------------- INTEGER :: IRWK, ICWK !----------------------------------------------- write(*,*) 'fish90: genbun' IERROR = 0 ! check input arguments IF (M <= 2) then ierror = 1 return end if IF (N <= 2) then ierror = 2 return end if IF (IDIMY < M) then ierror = 3 return end if IF (NPEROD<0 .OR. NPEROD>4) then ierror = 4 return end if IF (MPEROD<0 .OR. MPEROD>1) then ierror = 5 return end if ! compute and allocate real work space for genbun CALL GEN_SPACE (N, M, IRWK) ICWK = 0 CALL ALLOCATFISH (IRWK, ICWK, W, IERROR) ! return if allocation failed (e.g., if n,m are too large) IF (IERROR == 20) THEN write(*,*) 'error call ALLOCATFISH' RETURN END IF call genbunn(NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y,IERROR,w%rew,IRWK) ! release allocated work space CALL FISHFIN (W,IERROR) IF (IERROR == 20) THEN write(*,*) 'error call FISHFIN' RETURN END IF ! END SUBROUTINE GENBUN SUBROUTINE GENBUNN(NPEROD,N,MPEROD,M,vecA,vecB,vecC,IDIMY,matY, & IERROR,vecW,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(INOUT) :: IERROR DOUBLE PRECISION, DIMENSION(M), INTENT(IN) :: vecA, vecB, vecC DOUBLE PRECISION, DIMENSION(IW), INTENT(OUT) :: vecW DOUBLE PRECISION, DIMENSION(IDIMY,N),INTENT(INOUT) :: matY !----------------------------------------------- ! L o c a l V a r i a b l e s !----------------------------------------------- INTEGER :: I, MP1, IWBA, IWBB, IWBC, IWB2, IWB3, IWW1, IWW2, IWW3, & IWD, IWTCOS, IWP, IW2, K, J, MP, NP, & IPSTOR, IREV, MH, MHM1, MODD, NBY2, MSKIP DOUBLE PRECISION :: A1 !----------------------------------------------- IF (MPEROD /= 1) THEN DO I = 2, M IF (vecA(I) /= vecC(1)) GO TO 103 IF (vecC(I) /= vecC(1)) GO TO 103 IF (vecB(I) /= vecB(1)) GO TO 103 END DO GO TO 104 ENDIF IF (vecA(1)/=ZERO .OR. vecC(M)/=ZERO) IERROR = 7 GO TO 104 103 CONTINUE IERROR = 6 104 CONTINUE IF (IERROR /= 0) RETURN MP1 = M + 1 IWBA = MP1 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 vecW(IWBA:M-1+IWBA) = -vecA(:M) vecW(IWBC:M-1+IWBC) = -vecC(:M) vecW(IWBB:M-1+IWBB) = TWO - vecB(:M) matY(:M,:N) = -matY(:M,:N) MP = MPEROD + 1 NP = NPEROD + 1 GO TO (114,107) MP 107 CONTINUE GO TO (108,109,110,111,123) NP 108 CONTINUE IW2 = IW - IWP + 1 CALL POISP2 (M, N, vecW(IWBA:IWBA+M-1), vecW(IWBB:IWBB+M-1), & vecW(IWBC:IWBC+M-1), matY, IDIMY, & vecW(1:M), vecW(IWB2:IWB2+M-1), vecW(IWB3:IWB3+M-1), & vecW(IWW1:IWW1+M-1), vecW(IWW2:IWW2+M-1), & vecW(IWW3:IWW3+M-1), vecW(IWD:IWD+M-1), & vecW(IWTCOS:IWTCOS+4*N-1), vecW(IWP:),IW2) GO TO 112 109 CONTINUE IW2 = IW - IWP + 1 CALL POISD2 (M, N, 1, vecW(IWBA:IWBA+M-1), vecW(IWBB:IWBB+M-1), & vecW(IWBC::IWBC+M-1), matY, IDIMY, & vecW(1:M), vecW(IWW1:IWW1+M-1), vecW(IWD:IWD+M-1), & vecW(IWTCOS:IWTCOS+4*N-1), vecW(IWP:), IW2) GO TO 112 110 CONTINUE IW2 = IW - IWP + 1 CALL POISN2 (M, N, 1, 2, vecW(IWBA:IWBA+M-1), vecW(IWBB:IWBB+M-1), & vecW(IWBC:IWBC+M-1), matY, IDIMY, vecW(1:M), & vecW(IWB2:IWB2+M-1), vecW(IWB3:IWB3+M-1), & vecW(IWW1:IWW1+M-1), vecW(IWW2:IWW2+M-1), & vecW(IWW3::IWW3+M-1), vecW(IWD:IWD+M-1), & vecW(IWTCOS:IWTCOS+4*N-1),vecW(IWP:),IW2) GO TO 112 111 CONTINUE IW2 = IW - IWP + 1 CALL POISN2 (M, N, 1, 1, vecW(IWBA:IWBA+M-1), vecW(IWBB:IWBB+M-1), & vecW(IWBC::IWBC+M-1), matY, IDIMY, vecW(1:M), & vecW(IWB2:IWW1+M-1), vecW(IWB3:IWB3+M-1), & vecW(IWW1:IWW1+M-1), vecW(IWW2:IWW2+M-1), & vecW(IWW3::IWW3+M-1), vecW(IWD:IWD+M-1), & vecW(IWTCOS:IWTCOS+4*N-1),vecW(IWP:),IW2) 112 CONTINUE IPSTOR = vecW(IWW1) IREV = 2 IF (NPEROD == 4) GO TO 124 113 CONTINUE GO TO (127,133) MP 114 CONTINUE MH = (M + 1)/2 MHM1 = MH - 1 MODD = 1 IF (MH*2 == M) MODD = 2 DO J = 1, N vecW(:MHM1) = matY(MH-1:MH-MHM1:(-1),J) - matY(MH+1:MHM1+MH,J) vecW(MH+1:MHM1+MH) = matY(MH-1:MH-MHM1:(-1),J) + matY(MH+1:MHM1+MH,J) vecW(MH) = TWO*matY(MH,J) GO TO (117,116) MODD 116 CONTINUE vecW(M) = TWO*matY(M,J) 117 CONTINUE matY(:M,J) = vecW(:M) END DO K = IWBC + MHM1 - 1 I = IWBA + MHM1 vecW(K) = ZERO vecW(I) = ZERO vecW(K+1) = TWO*vecW(K+1) SELECT CASE (MODD) CASE DEFAULT K = IWBB + MHM1 - 1 vecW(K) = vecW(K) - vecW(I-1) vecW(IWBC-1) = vecW(IWBC-1) + vecW(IWBB-1) CASE (2) vecW(IWBB-1) = vecW(K+1) END SELECT GO TO 107 ! ! REVERSE COLUMNS WHEN NPEROD = 4. ! 123 CONTINUE IREV = 1 NBY2 = N/2 124 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 (110,113) IREV 127 CONTINUE DO J = 1, N vecW(MH-1:MH-MHM1:(-1)) = HALF*(matY(MH+1:MHM1+MH,J)+matY(:MHM1,J)) vecW(MH+1:MHM1+MH) = HALF*(matY(MH+1:MHM1+MH,J)-matY(:MHM1,J)) vecW(MH) = HALF*matY(MH,J) GO TO (130,129) MODD 129 CONTINUE vecW(M) = HALF*matY(M,J) 130 CONTINUE matY(:M,J) = vecW(:M) END DO 133 CONTINUE vecW(1) = IPSTOR + IWP - 1 RETURN END SUBROUTINE GENBUNN SUBROUTINE POISD2(MR,NR,ISTAG,vecBA,vecBB,vecBC,matY,IDIMY, & vecB,vecW,vecD,TCOS,vecP,IDIMP) 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) :: MR, NR, ISTAG, IDIMY,IDIMP DOUBLE PRECISION, DIMENSION(MR), INTENT(IN) :: vecBA, vecBB, & vecBC DOUBLE PRECISION, DIMENSION(IDIMY,NR),INTENT(INOUT) :: matY DOUBLE PRECISION, DIMENSION(MR), INTENT(INOUT) :: vecB, vecD, vecW DOUBLE PRECISION, DIMENSION(4*NR), INTENT(INOUT) :: TCOS DOUBLE PRECISION, DIMENSION(IDIMP), INTENT(INOUT) :: vecP !----------------------------------------------- ! L o c a l V a r i a b l e s !----------------------------------------------- INTEGER :: M, N, JSH, IP, IPSTOR, KR, IRREG, JSTSAV, I, LR, NUN, & JST, JSP, L, NODD, J, JM1, JP1, JM2, JP2, JM3, JP3, & NODDPR, KRPI, IDEG, JDEG DOUBLE PRECISION :: FI, T !----------------------------------------------- ! ! SUBROUTINE TO SOLVE POISSON'S EQUATION FOR DIRICHLET BOUNDARY ! CONDITIONS. ! ! ISTAG = 1 IF THE LAST DIAGONAL BLOCK IS THE MATRIX A. ! ISTAG = 2 IF THE LAST DIAGONAL BLOCK IS THE MATRIX A+I. ! M = MR N = NR JSH = 0 FI = ONE/FLOAT(ISTAG) IP = -M IPSTOR = 0 SELECT CASE (ISTAG) CASE DEFAULT KR = 0 IRREG = 1 IF (N > 1) GO TO 106 TCOS(1) = ZERO CASE (2) KR = 1 JSTSAV = 1 IRREG = 2 IF (N > 1) GO TO 106 TCOS(1) = -ONE END SELECT vecB(:M) = matY(:M,1) CALL TRIX (1, 0, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW) matY(:M,1) = vecB(:M) GO TO 183 106 CONTINUE LR = 0 vecP(:M) = ZERO NUN = N JST = 1 JSP = N ! ! IRREG = 1 WHEN NO IRREGULARITIES HAVE OCCURRED, OTHERWISE IT IS 2. ! 108 CONTINUE L = 2*JST NODD = 2 - 2*((NUN + 1)/2) + NUN ! ! NODD = 1 WHEN NUN IS ODD, OTHERWISE IT IS 2. ! SELECT CASE (NODD) CASE DEFAULT JSP = JSP - L CASE (1) JSP = JSP - JST IF (IRREG /= 1) JSP = JSP - L END SELECT CALL COSGEN (JST, 1, HALF, ZERO, TCOS, 4*NR) IF (L <= JSP) THEN DO J = L, JSP, L JM1 = J - JSH JP1 = J + JSH JM2 = J - JST JP2 = J + JST JM3 = JM2 - JSH JP3 = JP2 + JSH IF (JST == 1) THEN vecB(:M) = TWO*matY(:M,J) matY(:M,J) = matY(:M,JM2) + matY(:M,JP2) ELSE DO I = 1, M T = matY(I,J) - matY(I,JM1) - matY(I,JP1) + matY(I,JM2) + & matY(I,JP2) vecB(I) = T + matY(I,J) - matY(I,JM3) - matY(I,JP3) matY(I,J) = T END DO ENDIF CALL TRIX (JST, 0, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW) matY(:M,J) = matY(:M,J) + vecB(:M) END DO ENDIF ! ! REDUCTION FOR LAST UNKNOWN ! SELECT CASE (NODD) CASE DEFAULT GO TO (152,120) IRREG ! ! ODD NUMBER OF UNKNOWNS ! 120 CONTINUE JSP = JSP + L J = JSP JM1 = J - JSH JP1 = J + JSH JM2 = J - JST JP2 = J + JST JM3 = JM2 - JSH GO TO (123,121) ISTAG 121 CONTINUE IF (JST /= 1) GO TO 123 vecB(:M) = matY(:M,J) matY(:M,J) = ZERO GO TO 130 123 CONTINUE SELECT CASE (NODDPR) CASE DEFAULT vecB(:M) = HALF*(matY(:M,JM2)-matY(:M,JM1)-matY(:M,JM3)) + & vecP(IP+1:M+IP)+ matY(:M,J) CASE (2) vecB(:M) = HALF*(matY(:M,JM2)-matY(:M,JM1)-matY(:M,JM3)) + & matY(:M,JP2) - matY(:M,JP1) + matY(:M,J) END SELECT matY(:M,J) = HALF*(matY(:M,J)-matY(:M,JM1)-matY(:M,JP1)) 130 CONTINUE CALL TRIX (JST, 0, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW) IP = IP + M IPSTOR = MAX0(IPSTOR,IP + M) vecP(IP+1:M+IP) = matY(:M,J) + vecB(:M) vecB(:M) = matY(:M,JP2) + vecP(IP+1:M+IP) IF (LR == 0) THEN DO I = 1, JST KRPI = KR + I TCOS(KRPI) = TCOS(I) END DO ELSE CALL COSGEN (LR, JSTSAV, ZERO, FI, TCOS(JST+1), 4*NR-JST) CALL MERGE (TCOS, 0, JST, JST, LR, KR, 4*NR) ENDIF CALL COSGEN (KR, JSTSAV, ZERO, FI, TCOS, 4*NR) CALL TRIX (KR, KR, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW) matY(:M,J) = matY(:M,JM2) + vecB(:M) + vecP(IP+1:M+IP) LR = KR KR = KR + L ! ! EVEN NUMBER OF UNKNOWNS ! CASE (2) JSP = JSP + L J = JSP JM1 = J - JSH JP1 = J + JSH JM2 = J - JST JP2 = J + JST JM3 = JM2 - JSH SELECT CASE (IRREG) CASE DEFAULT JSTSAV = JST IDEG = JST KR = L CASE (2) CALL COSGEN (KR, JSTSAV, ZERO, FI, TCOS, 4*NR) CALL COSGEN (LR, JSTSAV, ZERO, FI, TCOS(KR+1), 4*NR-KR) IDEG = KR KR = KR + JST END SELECT IF (JST == 1) THEN IRREG = 2 vecB(:M) = matY(:M,J) matY(:M,J) = matY(:M,JM2) ELSE vecB(:M) = matY(:M,J) + HALF*(matY(:M,JM2)-matY(:M,JM1)- & matY(:M,JM3)) SELECT CASE (IRREG) CASE DEFAULT matY(:M,J) = matY(:M,JM2) + HALF*(matY(:M,J)- & matY(:M,JM1)-matY(:M,JP1)) IRREG = 2 CASE (2) SELECT CASE (NODDPR) CASE DEFAULT matY(:M,J) = matY(:M,JM2) + vecP(IP+1:M+IP) IP = IP - M CASE (2) matY(:M,J) = matY(:M,JM2) + matY(:M,J) - matY(:M,JM1) END SELECT END SELECT ENDIF CALL TRIX (IDEG, LR, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW) matY(:M,J) = matY(:M,J) + vecB(:M) END SELECT 152 CONTINUE NUN = NUN/2 NODDPR = NODD JSH = JST JST = 2*JST IF (NUN >= 2) GO TO 108 ! ! START SOLUTION. ! J = JSP vecB(:M) = matY(:M,J) SELECT CASE (IRREG) CASE DEFAULT CALL COSGEN (JST, 1, HALF, ZERO, TCOS, 4*NR) IDEG = JST CASE (2) KR = LR + JST CALL COSGEN (KR, JSTSAV, ZERO, FI, TCOS, 4*NR) CALL COSGEN (LR, JSTSAV, ZERO, FI, TCOS(KR+1), 4*NR-KR) IDEG = KR END SELECT CALL TRIX (IDEG, LR, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW) JM1 = J - JSH JP1 = J + JSH SELECT CASE (IRREG) CASE DEFAULT matY(:M,J) = HALF*(matY(:M,J)-matY(:M,JM1)-matY(:M,JP1)) + vecB(:M) CASE (2) SELECT CASE (NODDPR) CASE DEFAULT matY(:M,J) = vecP(IP+1:M+IP) + vecB(:M) IP = IP - M CASE (2) matY(:M,J) = matY(:M,J) - matY(:M,JM1) + vecB(:M) END SELECT END SELECT 164 CONTINUE JST = JST/2 JSH = JST/2 NUN = 2*NUN IF (NUN > N) GO TO 183 DO J = JST, N, L JM1 = J - JSH JP1 = J + JSH JM2 = J - JST JP2 = J + JST IF (J <= JST) THEN vecB(:M) = matY(:M,J) + matY(:M,JP2) ELSE IF (JP2 <= N) GO TO 168 vecB(:M) = matY(:M,J) + matY(:M,JM2) IF (JST < JSTSAV) IRREG = 1 GO TO (170,171) IRREG 168 CONTINUE vecB(:M) = matY(:M,J) + matY(:M,JM2) + matY(:M,JP2) ENDIF 170 CONTINUE CALL COSGEN (JST, 1, HALF, ZERO, TCOS, 4*NR) IDEG = JST JDEG = 0 GO TO 172 171 CONTINUE IF (J + L > N) LR = LR - JST KR = JST + LR CALL COSGEN (KR, JSTSAV, ZERO, FI, TCOS, 4*NR) CALL COSGEN (LR, JSTSAV, ZERO, FI, TCOS(KR+1), 4*NR-KR) IDEG = KR JDEG = LR 172 CONTINUE CALL TRIX (IDEG, JDEG, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW) IF (JST <= 1) THEN matY(:M,J) = vecB(:M) ELSE IF (JP2 > N) GO TO 177 175 CONTINUE matY(:M,J) = HALF*(matY(:M,J)-matY(:M,JM1)-matY(:M,JP1)) + vecB(:M) CYCLE 177 CONTINUE GO TO (175,178) IRREG 178 CONTINUE IF (J + JSH <= N) THEN matY(:M,J) = vecB(:M) + vecP(IP+1:M+IP) IP = IP - M ELSE matY(:M,J) = vecB(:M) + matY(:M,J) - matY(:M,JM1) ENDIF ENDIF END DO L = L/2 GO TO 164 183 CONTINUE vecW(1) = IPSTOR RETURN END SUBROUTINE POISD2 SUBROUTINE POISN2(M, N, ISTAG, MIXBND, vecA, vecBB, vecC, & matQ, IDIMQ, vecB, vecB2,vecB3, & vecW, vecW2, vecW3, vecD, TCOS, vecP, IDIMP) 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) :: M, N, ISTAG, MIXBND, IDIMQ,IDIMP DOUBLE PRECISION, DIMENSION(M), INTENT(IN) :: vecA, vecBB, vecC DOUBLE PRECISION, DIMENSION(IDIMQ,N), INTENT(INOUT) :: matQ DOUBLE PRECISION, DIMENSION(M),INTENT(INOUT) :: vecB, vecB2, vecB3, & vecD, vecW, vecW2, vecW3 DOUBLE PRECISION, DIMENSION(IDIMP),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, MR, IP, IPSTOR, I2R, JR, NR, NLAST, & KR, LR, I, NROD, JSTART, JSTOP, I2RBY2, & J, JP1, JP2, JP3, JM1,JM2, JM3, NRODPR, II, I1, I2, & JR2, NLASTP, JSTEP DOUBLE PRECISION :: FISTAG, FNUM, FDEN, FI, T !----------------------------------------------- ! ! SUBROUTINE TO SOLVE POISSON'S EQUATION WITH NEUMANN BOUNDARY ! CONDITIONS. ! ! ISTAG = 1 IF THE LAST DIAGONAL BLOCK IS A. ! ISTAG = 2 IF THE LAST DIAGONAL BLOCK IS A-I. ! MIXBND = 1 IF HAVE NEUMANN BOUNDARY CONDITIONS AT BOTH BOUNDARIES. ! MIXBND = 2 IF HAVE NEUMANN BOUNDARY CONDITIONS AT BOTTOM AND ! DIRICHLET CONDITION AT TOP. (FOR THIS CASE, MUST HAVE ISTAG = 1.) ! EQUIVALENCE (K(1), K1), (K(2), K2), (K(3), K3), (K(4), K4) FISTAG = 3 - ISTAG FNUM = ONE/FLOAT(ISTAG) FDEN = HALF*FLOAT(ISTAG - 1) MR = M IP = -MR IPSTOR = 0 I2R = 1 JR = 2 NR = N NLAST = N KR = 1 LR = 0 GO TO (101,103) ISTAG 101 CONTINUE matQ(:MR,N) = HALF*matQ(:MR,N) GO TO (103,104) MIXBND 103 CONTINUE IF (N <= 3) GO TO 155 104 CONTINUE JR = 2*I2R NROD = 1 IF ((NR/2)*2 == NR) NROD = 0 SELECT CASE (MIXBND) CASE DEFAULT JSTART = 1 CASE (2) JSTART = JR NROD = 1 - NROD END SELECT JSTOP = NLAST - JR IF (NROD == 0) JSTOP = JSTOP - I2R CALL COSGEN (I2R, 1, HALF, ZERO, TCOS, 4*NR) I2RBY2 = I2R/2 IF (JSTOP < JSTART) THEN J = JR ELSE 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 JM1 = JP1 JM2 = JP2 JM3 = JP3 ENDIF IF (I2R == 1) THEN IF (J == 1) JM2 = JP2 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 CALL TRIX (I2R, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, 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) = FISTAG*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) THEN CALL COSGEN (LR, 1, HALF, FDEN, TCOS(KR+1), 4*NR-KR) ELSE vecB(:MR) = FISTAG*vecB(:MR) ENDIF ENDIF CALL COSGEN (KR, 1, HALF, FDEN, TCOS, 4*NR) CALL TRIX (KR, LR, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, 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) CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) IP = 0 IPSTOR = MR SELECT CASE (ISTAG) CASE DEFAULT vecP(:MR) = vecB(:MR) vecB(:MR) = vecB(:MR) + matQ(:MR,N) TCOS(1) = ONE TCOS(2) = ZERO CALL TRIX (1, 1, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,J) = matQ(:MR,JM2) + vecP(:MR) + vecB(:MR) GO TO 150 CASE (1) vecP(:MR) = vecB(:MR) matQ(:MR,J) = matQ(:MR,JM2) + TWO*matQ(:MR,JP2) + 3.*vecB(:MR) GO TO 150 END SELECT ENDIF 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 TRIX (I2R, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, 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, HALF, FDEN, TCOS(I2R+1), 4*NR-I2R) CALL MERGE (TCOS, 0, I2R, I2R, LR, KR, 4*NR) ELSE DO I = 1, I2R II = KR + I TCOS(II) = TCOS(I) END DO ENDIF CALL COSGEN (KR, 1, HALF, FDEN, TCOS, 4*NR) IF (LR == 0) THEN GO TO (146,145) ISTAG ENDIF 145 CONTINUE CALL TRIX (KR, KR, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) GO TO 148 146 CONTINUE vecB(:MR) = FISTAG*vecB(:MR) 148 CONTINUE matQ(:MR,J) = matQ(:MR,JM2) + vecP(IP+1:MR+IP) + vecB(:MR) 150 CONTINUE LR = KR KR = KR + JR ENDIF SELECT CASE (MIXBND) CASE DEFAULT NR = (NLAST - 1)/JR + 1 IF (NR <= 3) GO TO 155 CASE (2) NR = NLAST/JR IF (NR <= 1) GO TO 192 END SELECT I2R = JR NRODPR = NROD GO TO 104 155 CONTINUE J = 1 + JR JM1 = J - I2R JP1 = J + I2R JM2 = NLAST - I2R IF (NR /= 2) THEN IF (LR /= 0) GO TO 170 IF (N == 3) THEN ! ! CASE N = 3. ! GO TO (156,168) ISTAG 156 CONTINUE vecB(:MR) = matQ(:MR,2) TCOS(1) = ZERO CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,2) = vecB(:MR) vecB(:MR) = 4.*vecB(:MR) + matQ(:MR,1) + TWO*matQ(:MR,3) TCOS(1) = -TWO TCOS(2) = TWO I1 = 2 I2 = 0 CALL TRIX (I1, I2, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,2) = matQ(:MR,2) + vecB(:MR) vecB(:MR) = matQ(:MR,1) + TWO*matQ(:MR,2) TCOS(1) = ZERO CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,1) = vecB(:MR) JR = 1 I2R = 0 GO TO 194 ENDIF ! ! CASE N = 2**P+1 ! GO TO (162,170) ISTAG 162 CONTINUE vecB(:MR) = matQ(:MR,J) + HALF*matQ(:MR,1) - & matQ(:MR,JM1) + matQ(:MR,NLAST) - matQ(:MR,JM2) CALL COSGEN (JR, 1, HALF, ZERO, TCOS, 4*NR) CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,J) = HALF*(matQ(:MR,J)-matQ(:MR,JM1)-matQ(:MR,JP1)) + vecB(:MR) vecB(:MR) = matQ(:MR,1) + TWO*matQ(:MR,NLAST) + 4.*matQ(:MR,J) JR2 = 2*JR CALL COSGEN (JR, 1, ZERO, ZERO, TCOS, 4*NR) TCOS(JR+1:JR*2) = -TCOS(JR:1:(-1)) CALL TRIX (JR2, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,J) = matQ(:MR,J) + vecB(:MR) vecB(:MR) = matQ(:MR,1) + TWO*matQ(:MR,J) CALL COSGEN (JR, 1, HALF, ZERO, TCOS, 4*NR) CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,1) = HALF*matQ(:MR,1) - matQ(:MR,JM1) + vecB(:MR) GO TO 194 ! ! CASE OF GENERAL N WITH NR = 3 . ! 168 CONTINUE vecB(:MR) = matQ(:MR,2) matQ(:MR,2) = ZERO vecB2(:MR) = matQ(:MR,3) vecB3(:MR) = matQ(:MR,1) JR = 1 I2R = 0 J = 2 GO TO 177 170 CONTINUE vecB(:MR) = HALF*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) + TWO*T END DO 177 CONTINUE K1 = KR + 2*JR - 1 K2 = KR + JR TCOS(K1+1) = -TWO K4 = K1 + 3 - ISTAG CALL COSGEN (K2 + ISTAG - 2, 1, ZERO, FNUM, TCOS(K4), 4*NR-K4+1) K4 = K1 + K2 + 1 CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(K4), 4*NR-K4+1) CALL MERGE (TCOS, K1, K2, K1 + K2, JR - 1, 0, 4*NR) K3 = K1 + K2 + LR CALL COSGEN (JR, 1, HALF, ZERO, TCOS(K3+1), 4*NR-K3) K4 = K3 + JR + 1 CALL COSGEN (KR, 1, HALF, FDEN, TCOS(K4), 4*NR-K4+1) CALL MERGE (TCOS, K3, JR, K3 + JR, KR, K1, 4*NR) IF (LR /= 0) THEN CALL COSGEN (LR, 1, HALF, FDEN, TCOS(K4), 4*NR-K4+1) CALL MERGE (TCOS, K3, JR, K3 + JR, LR, K3 - LR, 4*NR) CALL COSGEN (KR, 1, HALF, FDEN, TCOS(K4), 4*NR-K4+1) ENDIF K3 = KR K4 = KR CALL TRI3 (MR, vecA, vecBB, vecC, K, vecB, vecB2, vecB3, & TCOS, 4*NR, vecD, vecW, vecW2, vecW3) vecB(:MR) = vecB(:MR) + vecB2(:MR) + vecB3(:MR) TCOS(1) = TWO CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,J) = matQ(:MR,J) + vecB(:MR) vecB(:MR) = matQ(:MR,1) + TWO*matQ(:MR,J) CALL COSGEN (JR, 1, HALF, ZERO, TCOS, 4*NR) CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) IF (JR == 1) THEN matQ(:MR,1) = vecB(:MR) GO TO 194 ENDIF matQ(:MR,1) = HALF*matQ(:MR,1) - matQ(:MR,JM1) + vecB(:MR) GO TO 194 ENDIF IF (N == 2) THEN ! ! CASE N = 2 ! vecB(:MR) = matQ(:MR,1) TCOS(1) = ZERO CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,1) = vecB(:MR) vecB(:MR) = TWO*(matQ(:MR,2)+vecB(:MR))*FISTAG TCOS(1) = -FISTAG TCOS(2) = TWO CALL TRIX (2, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,1) = matQ(:MR,1) + vecB(:MR) JR = 1 I2R = 0 GO TO 194 ENDIF vecB3(:MR) = ZERO vecB(:MR) = matQ(:MR,1) + TWO*vecP(IP+1:MR+IP) matQ(:MR,1) = HALF*matQ(:MR,1) - matQ(:MR,JM1) vecB2(:MR) = TWO*(matQ(:MR,1)+matQ(:MR,NLAST)) K1 = KR + JR - 1 TCOS(K1+1) = -TWO K4 = K1 + 3 - ISTAG CALL COSGEN (KR + ISTAG - 2, 1, ZERO, FNUM, TCOS(K4), 4*NR-K4+1) K4 = K1 + KR + 1 CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(K4), 4*NR-K4+1) CALL MERGE (TCOS, K1, KR, K1 + KR, JR - 1, 0, 4*NR) CALL COSGEN (KR, 1, HALF, FDEN, TCOS(K1+1), 4*NR-K1) K2 = KR K4 = K1 + K2 + 1 CALL COSGEN (LR, 1, HALF, FDEN, TCOS(K4), 4*NR-K4+1) K3 = LR K4 = 0 CALL TRI3 (MR, vecA, vecBB, vecC, K, vecB, vecB2, vecB3, TCOS, & 4*NR, vecD, vecW, vecW2, vecW3) vecB(:MR) = vecB(:MR) + vecB2(:MR) TCOS(1) = TWO CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,1) = matQ(:MR,1) + vecB(:MR) GO TO 194 192 CONTINUE vecB(:MR) = matQ(:MR,NLAST) GO TO 196 194 CONTINUE J = NLAST - JR vecB(:MR) = matQ(:MR,NLAST) + matQ(:MR,J) 196 CONTINUE JM2 = NLAST - I2R IF (JR == 1) THEN matQ(:MR,NLAST) = ZERO 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, HALF, FDEN, TCOS, 4*NR) CALL COSGEN (LR, 1, HALF, FDEN, TCOS(KR+1), 4*NR-KR) IF (LR == 0) THEN vecB(:MR) = FISTAG*vecB(:MR) ENDIF CALL TRIX (KR, LR, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW) matQ(:MR,NLAST) = matQ(:MR,NLAST) + vecB(:MR) NLASTP = NLAST 206 CONTINUE JSTEP = JR JR = I2R I2R = I2R/2 IF (JR == 0) GO TO 222 SELECT CASE (MIXBND) CASE DEFAULT JSTART = 1 + JR CASE (2) JSTART = JR END SELECT 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, 4*NR) 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) = ZERO 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, 4*NR, 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 194 GO TO 206 222 CONTINUE vecW(1) = IPSTOR RETURN END SUBROUTINE POISN2 SUBROUTINE POISP2(M, N, vecA, vecBB, vecC, matQ, IDIMQ, & vecB, vecB2, vecB3, vecW, vecW2, vecW3, & vecD, TCOS, vecP, IP) 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) :: M, N, IDIMQ, IP DOUBLE PRECISION, DIMENSION(M), INTENT(IN) :: vecA, vecBB, vecC DOUBLE PRECISION, DIMENSION(4*N), INTENT(INOUT) :: TCOS DOUBLE PRECISION, DIMENSION(IDIMQ,N), INTENT(INOUT) :: matQ DOUBLE PRECISION, DIMENSION(IP),INTENT(INOUT) :: vecP DOUBLE PRECISION, DIMENSION(M),INTENT(INOUT) :: vecB, vecB2, vecB3, & vecD, vecW, vecW2, vecW3 !----------------------------------------------- ! L o c a l V a r i a b l e s !----------------------------------------------- INTEGER :: MR, NR, NRM1, J, NRMJ, NRPJ, I, IPSTOR, LH DOUBLE PRECISION :: S, T !----------------------------------------------- ! ! SUBROUTINE TO SOLVE POISSON EQUATION WITH PERIODIC BOUNDARY ! CONDITIONS. ! MR = M NR = (N + 1)/2 NRM1 = NR - 1 IF (2*NR == N) THEN ! ! EVEN NUMBER OF UNKNOWNS ! DO J = 1, NRM1 NRMJ = NR - J NRPJ = NR + J DO I = 1, MR S = matQ(I,NRMJ) - matQ(I,NRPJ) T = matQ(I,NRMJ) + matQ(I,NRPJ) matQ(I,NRMJ) = S matQ(I,NRPJ) = T END DO END DO matQ(:MR,NR) = TWO*matQ(:MR,NR) matQ(:MR,N) = TWO*matQ(:MR,N) CALL POISD2 (MR, NRM1, 1, vecA, vecBB, vecC, matQ, IDIMQ, & vecB, vecW, vecD, TCOS, vecP, IP) IPSTOR = vecW(1) CALL POISN2 (MR, NR + 1, 1, 1, vecA, vecBB, vecC, matQ(:,NR:), & IDIMQ, vecB, vecB2, vecB3, vecW, vecW2, vecW3, vecD, & TCOS, vecP, IP) IPSTOR = MAX0(IPSTOR,INT(vecW(1))) DO J = 1, NRM1 NRMJ = NR - J NRPJ = NR + J DO I = 1, MR S = HALF*(matQ(I,NRPJ)+matQ(I,NRMJ)) T = HALF*(matQ(I,NRPJ)-matQ(I,NRMJ)) matQ(I,NRMJ) = S matQ(I,NRPJ) = T END DO END DO matQ(:MR,NR) = HALF*matQ(:MR,NR) matQ(:MR,N) = HALF*matQ(:MR,N) ELSE DO J = 1, NRM1 NRPJ = N + 1 - J DO I = 1, MR S = matQ(I,J) - matQ(I,NRPJ) T = matQ(I,J) + matQ(I,NRPJ) matQ(I,J) = S matQ(I,NRPJ) = T END DO END DO matQ(:MR,NR) = TWO*matQ(:MR,NR) LH = NRM1/2 DO J = 1, LH NRMJ = NR - J DO I = 1, MR S = matQ(I,J) matQ(I,J) = matQ(I,NRMJ) matQ(I,NRMJ) = S END DO END DO CALL POISD2 (MR, NRM1, 2, vecA, vecBB, vecC, matQ, IDIMQ, & vecB, vecW, vecD, TCOS, vecP, IP) IPSTOR = vecW(1) CALL POISN2 (MR, NR, 2, 1, vecA, vecBB, vecC, matQ(:,NR:), IDIMQ, & vecB, vecB2, vecB3, vecW, vecW2, vecW3, vecD, TCOS, & vecP, IP) IPSTOR = MAX0(IPSTOR,INT(vecW(1))) DO J = 1, NRM1 NRPJ = NR + J DO I = 1, MR S = HALF*(matQ(I,NRPJ)+matQ(I,J)) T = HALF*(matQ(I,NRPJ)-matQ(I,J)) matQ(I,NRPJ) = T matQ(I,J) = S END DO END DO matQ(:MR,NR) = HALF*matQ(:MR,NR) DO J = 1, LH NRMJ = NR - J DO I = 1, MR S = matQ(I,J) matQ(I,J) = matQ(I,NRMJ) matQ(I,NRMJ) = S END DO END DO ENDIF 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 POISP2