c Fishpack calls routines from fftpack. If you don't already have c those routines, submit your request in the form c send blktri from fishpack fftpack Careful! Anything free comes with no guarantee. C FISHPAK2 FROM PORTLIB 12/30/83 SUBROUTINE HWSCRT (A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC,BDD, 1 ELMBDA,F,IDIMF,PERTRB,IERROR,W) c c This is the fast Poisson solver used with the nitsol test codes. c c This code was originally single precision. Quick and dirty conversion c to double precision done on 1/30/90. --HFW c Further mods to change all declarations of arrays of length 1 to c assumed-size arrays. 3/96 --MAP C C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * F I S H P A K * C * * C * * C * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF * C * * C * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS * C * * C * (VERSION 3.1 , OCTOBER 1980) * C * * C * BY * C * * C * JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET * C * * C * OF * C * * C * THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH * C * * C * BOULDER, COLORADO (80307) U.S.A. * C * * C * WHICH IS SPONSORED BY * C * * C * THE NATIONAL SCIENCE FOUNDATION * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C C * * * * * * * * * PURPOSE * * * * * * * * * * * * * * * * * * C C SUBROUTINE HWSCRT SOLVES THE STANDARD FIVE-POINT FINITE C DIFFERENCE APPROXIMATION TO THE HELMHOLTZ EQUATION IN CARTESIAN C COORDINATES: C C (D/DX)(DU/DX) + (D/DY)(DU/DY) + LAMBDA*U = F(X,Y). C C C C * * * * * * * * PARAMETER DESCRIPTION * * * * * * * * * * C C * * * * * * ON INPUT * * * * * * C C A,B C THE RANGE OF X, I.E., A .LE. X .LE. B. A MUST BE LESS THAN B. C C M C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (A,B) IS C SUBDIVIDED. HENCE, THERE WILL BE M+1 GRID POINTS IN THE C X-DIRECTION GIVEN BY X(I) = A+(I-1)DX FOR I = 1,2,...,M+1, C WHERE DX = (B-A)/M IS THE PANEL WIDTH. M MUST BE GREATER THAN 3. C C MBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS AT X = A AND X = B. C C = 0 IF THE SOLUTION IS PERIODIC IN X, I.E., U(I,J) = U(M+I,J). C = 1 IF THE SOLUTION IS SPECIFIED AT X = A AND X = B. C = 2 IF THE SOLUTION IS SPECIFIED AT X = A AND THE DERIVATIVE OF C THE SOLUTION WITH RESPECT TO X IS SPECIFIED AT X = B. C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO X IS C SPECIFIED AT X = A AND X = B. C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO X IS C SPECIFIED AT X = A AND THE SOLUTION IS SPECIFIED AT X = B. C C BDA C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO X AT X = A. C WHEN MBDCND = 3 OR 4, C C BDA(J) = (D/DX)U(A,Y(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDA IS A DUMMY VARIABLE. C C BDB C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO X AT X = B. C WHEN MBDCND = 2 OR 3, C C BDB(J) = (D/DX)U(B,Y(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE BDB IS A DUMMY VARIABLE. C C C,D C THE RANGE OF Y, I.E., C .LE. Y .LE. D. C MUST BE LESS THAN D. C C N C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (C,D) IS C SUBDIVIDED. HENCE, THERE WILL BE N+1 GRID POINTS IN THE C Y-DIRECTION GIVEN BY Y(J) = C+(J-1)DY FOR J = 1,2,...,N+1, WHERE C DY = (D-C)/N IS THE PANEL WIDTH. N MUST BE GREATER THAN 3. C C NBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS AT Y = C AND Y = D. C C = 0 IF THE SOLUTION IS PERIODIC IN Y, I.E., U(I,J) = U(I,N+J). C = 1 IF THE SOLUTION IS SPECIFIED AT Y = C AND Y = D. C = 2 IF THE SOLUTION IS SPECIFIED AT Y = C AND THE DERIVATIVE OF C THE SOLUTION WITH RESPECT TO Y IS SPECIFIED AT Y = D. C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Y IS C SPECIFIED AT Y = C AND Y = D. C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Y IS C SPECIFIED AT Y = C AND THE SOLUTION IS SPECIFIED AT Y = D. C C BDC C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Y AT Y = C. C WHEN NBDCND = 3 OR 4, C C BDC(I) = (D/DY)U(X(I),C), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDC IS A DUMMY VARIABLE. C C BDD C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Y AT Y = D. C WHEN NBDCND = 2 OR 3, C C BDD(I) = (D/DY)U(X(I),D), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDD IS A DUMMY VARIABLE. C C ELMBDA C THE CONSTANT LAMBDA IN THE HELMHOLTZ EQUATION. IF C LAMBDA .GT. 0, A SOLUTION MAY NOT EXIST. HOWEVER, HWSCRT WILL C ATTEMPT TO FIND A SOLUTION. C C F C A TWO-DIMENSIONAL ARRAY WHICH SPECIFIES THE VALUES OF THE RIGHT C SIDE OF THE HELMHOLTZ EQUATION AND BOUNDARY VALUES (IF ANY). C FOR I = 2,3,...,M AND J = 2,3,...,N C C F(I,J) = F(X(I),Y(J)). C C ON THE BOUNDARIES F IS DEFINED BY C C MBDCND F(1,J) F(M+1,J) C ------ --------- -------- C C 0 F(A,Y(J)) F(A,Y(J)) C 1 U(A,Y(J)) U(B,Y(J)) C 2 U(A,Y(J)) F(B,Y(J)) J = 1,2,...,N+1 C 3 F(A,Y(J)) F(B,Y(J)) C 4 F(A,Y(J)) U(B,Y(J)) C C C NBDCND F(I,1) F(I,N+1) C ------ --------- -------- C C 0 F(X(I),C) F(X(I),C) C 1 U(X(I),C) U(X(I),D) C 2 U(X(I),C) F(X(I),D) I = 1,2,...,M+1 C 3 F(X(I),C) F(X(I),D) C 4 F(X(I),C) U(X(I),D) C C F MUST BE DIMENSIONED AT LEAST (M+1)*(N+1). C C NOTE C C IF THE TABLE CALLS FOR BOTH THE SOLUTION U AND THE RIGHT SIDE F C AT A CORNER THEN THE SOLUTION MUST BE SPECIFIED. C C IDIMF C THE ROW (OR FIRST) DIMENSION OF THE ARRAY F AS IT APPEARS IN THE C PROGRAM CALLING HWSCRT. THIS PARAMETER IS USED TO SPECIFY THE C VARIABLE DIMENSION OF F. IDIMF MUST BE AT LEAST M+1 . C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE USER FOR C WORK SPACE. W MAY REQUIRE UP TO 4*(N+1) + C (13 + INT(LOG2(N+1)))*(M+1) LOCATIONS. THE ACTUAL NUMBER OF C LOCATIONS USED IS COMPUTED BY HWSCRT AND IS RETURNED IN LOCATION C W(1). C C C * * * * * * ON OUTPUT * * * * * * C C F C CONTAINS THE SOLUTION U(I,J) OF THE FINITE DIFFERENCE C APPROXIMATION FOR THE GRID POINT (X(I),Y(J)), I = 1,2,...,M+1, C J = 1,2,...,N+1 . C C PERTRB C IF A COMBINATION OF PERIODIC OR DERIVATIVE BOUNDARY CONDITIONS C IS SPECIFIED FOR A POISSON EQUATION (LAMBDA = 0), A SOLUTION MAY C NOT EXIST. PERTRB IS A CONSTANT, CALCULATED AND SUBTRACTED FROM C F, WHICH ENSURES THAT A SOLUTION EXISTS. HWSCRT THEN COMPUTES C THIS SOLUTION, WHICH IS A LEAST SQUARES SOLUTION TO THE ORIGINAL C APPROXIMATION. THIS SOLUTION PLUS ANY CONSTANT IS ALSO A C SOLUTION. HENCE, THE SOLUTION IS NOT UNIQUE. THE VALUE OF C PERTRB SHOULD BE SMALL COMPARED TO THE RIGHT SIDE F. OTHERWISE, C A SOLUTION IS OBTAINED TO AN ESSENTIALLY DIFFERENT PROBLEM. C THIS COMPARISON SHOULD ALWAYS BE MADE TO INSURE THAT A C MEANINGFUL SOLUTION HAS BEEN OBTAINED. C C IERROR C AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS. EXCEPT C FOR NUMBERS 0 AND 6, A SOLUTION IS NOT ATTEMPTED. C C = 0 NO ERROR. C = 1 A .GE. B. C = 2 MBDCND .LT. 0 OR MBDCND .GT. 4 . C = 3 C .GE. D. C = 4 N .LE. 3 C = 5 NBDCND .LT. 0 OR NBDCND .GT. 4 . C = 6 LAMBDA .GT. 0 . C = 7 IDIMF .LT. M+1 . C = 8 M .LE. 3 C C SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSIBLY INCORRECT C CALL TO HWSCRT, THE USER SHOULD TEST IERROR AFTER THE CALL. C C W C W(1) CONTAINS THE REQUIRED LENGTH OF W. C C C * * * * * * * PROGRAM SPECIFICATIONS * * * * * * * * * * * * C C C DIMENSION OF BDA(N+1),BDB(N+1),BDC(M+1),BDD(M+1),F(IDIMF,N+1), C ARGUMENTS W(SEE ARGUMENT LIST) C C LATEST JUNE 1, 1976 C REVISION C C SUBPROGRAMS HWSCRT,GENBUN,POISD2,POISN2,POISP2,COSGEN,MERGE, C REQUIRED TRIX,TRI3,PIMACH C C SPECIAL NONE C CONDITIONS C C COMMON NONE C BLOCKS C C I/O NONE C C PRECISION SINGLE --- Converted to double ... HFW, 1/30/90. C C SPECIALIST ROLAND SWEET C C LANGUAGE FORTRAN C C HISTORY STANDARDIZED SEPTEMBER 1, 1973 C REVISED APRIL 1, 1976 C C ALGORITHM THE ROUTINE DEFINES THE FINITE DIFFERENCE C EQUATIONS, INCORPORATES BOUNDARY DATA, AND ADJUSTS C THE RIGHT SIDE OF SINGULAR SYSTEMS AND THEN CALLS C GENBUN TO SOLVE THE SYSTEM. C C SPACE 13110(OCTAL) = 5704(DECIMAL) LOCATIONS ON THE NCAR C REQUIRED CONTROL DATA 7600 C C TIMING AND THE EXECUTION TIME T ON THE NCAR CONTROL DATA C ACCURACY 7600 FOR SUBROUTINE HWSCRT IS ROUGHLY PROPORTIONAL C TO M*N*LOG2(N), BUT ALSO DEPENDS ON THE INPUT C PARAMETERS NBDCND AND MBDCND. SOME TYPICAL VALUES C ARE LISTED IN THE TABLE BELOW. C THE SOLUTION PROCESS EMPLOYED RESULTS IN A LOSS C OF NO MORE THAN THREE SIGNIFICANT DIGITS FOR N AND C M AS LARGE AS 64. MORE DETAILED INFORMATION ABOUT C ACCURACY CAN BE FOUND IN THE DOCUMENTATION FOR C SUBROUTINE GENBUN WHICH IS THE ROUTINE THAT C SOLVES THE FINITE DIFFERENCE EQUATIONS. C C C M(=N) MBDCND NBDCND T(MSECS) C ----- ------ ------ -------- C C 32 0 0 31 C 32 1 1 23 C 32 3 3 36 C 64 0 0 128 C 64 1 1 96 C 64 3 3 142 C C PORTABILITY AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN. C ALL MACHINE DEPENDENT CONSTANTS ARE LOCATED IN THE C FUNCTION PIMACH. C C REFERENCE SWARZTRAUBER,P. AND R. SWEET, 'EFFICIENT FORTRAN C SUBPROGRAMS FOR THE SOLUTION OF ELLIPTIC EQUATIONS' C NCAR TN/IA-109, JULY, 1975, 138 PP. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C implicit double precision (a-h,o-z) DIMENSION F(IDIMF,*) DIMENSION BDA(*) ,BDB(*) ,BDC(*) ,BDD(*) , 1 W(*) C C CHECK FOR INVALID PARAMETERS. C IERROR = 0 IF (A .GE. B) IERROR = 1 IF (MBDCND.LT.0 .OR. MBDCND.GT.4) IERROR = 2 IF (C .GE. D) IERROR = 3 IF (N .LE. 3) IERROR = 4 IF (NBDCND.LT.0 .OR. NBDCND.GT.4) IERROR = 5 IF (IDIMF .LT. M+1) IERROR = 7 IF (M .LE. 3) IERROR = 8 IF (IERROR .NE. 0) RETURN NPEROD = NBDCND MPEROD = 0 IF (MBDCND .GT. 0) MPEROD = 1 DELTAX = (B-A)/dfloat(M) TWDELX = 2.d0/DELTAX DELXSQ = 1.d0/DELTAX**2 DELTAY = (D-C)/dfloat(N) TWDELY = 2.d0/DELTAY DELYSQ = 1.d0/DELTAY**2 NP = NBDCND+1 NP1 = N+1 MP = MBDCND+1 MP1 = M+1 NSTART = 1 NSTOP = N NSKIP = 1 GO TO (104,101,102,103,104),NP 101 NSTART = 2 GO TO 104 102 NSTART = 2 103 NSTOP = NP1 NSKIP = 2 104 NUNK = NSTOP-NSTART+1 C C ENTER BOUNDARY DATA FOR X-BOUNDARIES. C MSTART = 1 MSTOP = M MSKIP = 1 GO TO (117,105,106,109,110),MP 105 MSTART = 2 GO TO 107 106 MSTART = 2 MSTOP = MP1 MSKIP = 2 107 DO 108 J=NSTART,NSTOP F(2,J) = F(2,J)-F(1,J)*DELXSQ 108 CONTINUE GO TO 112 109 MSTOP = MP1 MSKIP = 2 110 DO 111 J=NSTART,NSTOP F(1,J) = F(1,J)+BDA(J)*TWDELX 111 CONTINUE 112 GO TO (113,115),MSKIP 113 DO 114 J=NSTART,NSTOP F(M,J) = F(M,J)-F(MP1,J)*DELXSQ 114 CONTINUE GO TO 117 115 DO 116 J=NSTART,NSTOP F(MP1,J) = F(MP1,J)-BDB(J)*TWDELX 116 CONTINUE 117 MUNK = MSTOP-MSTART+1 C C ENTER BOUNDARY DATA FOR Y-BOUNDARIES. C GO TO (127,118,118,120,120),NP 118 DO 119 I=MSTART,MSTOP F(I,2) = F(I,2)-F(I,1)*DELYSQ 119 CONTINUE GO TO 122 120 DO 121 I=MSTART,MSTOP F(I,1) = F(I,1)+BDC(I)*TWDELY 121 CONTINUE 122 GO TO (123,125),NSKIP 123 DO 124 I=MSTART,MSTOP F(I,N) = F(I,N)-F(I,NP1)*DELYSQ 124 CONTINUE GO TO 127 125 DO 126 I=MSTART,MSTOP F(I,NP1) = F(I,NP1)-BDD(I)*TWDELY 126 CONTINUE C C MULTIPLY RIGHT SIDE BY DELTAY**2. C 127 DELYSQ = DELTAY*DELTAY DO 129 I=MSTART,MSTOP DO 128 J=NSTART,NSTOP F(I,J) = F(I,J)*DELYSQ 128 CONTINUE 129 CONTINUE C C DEFINE THE A,B,C COEFFICIENTS IN W-ARRAY. C ID2 = MUNK ID3 = ID2+MUNK ID4 = ID3+MUNK S = DELYSQ*DELXSQ ST2 = 2.d0*S DO 130 I=1,MUNK W(I) = S J = ID2+I W(J) = -ST2+ELMBDA*DELYSQ J = ID3+I W(J) = S 130 CONTINUE IF (MP .EQ. 1) GO TO 131 W(1) = 0.d0 W(ID4) = 0.d0 131 CONTINUE GO TO (135,135,132,133,134),MP 132 W(ID2) = ST2 GO TO 135 133 W(ID2) = ST2 134 W(ID3+1) = ST2 135 CONTINUE PERTRB = 0.d0 IF (ELMBDA) 144,137,136 136 IERROR = 6 GO TO 144 137 IF ((NBDCND.EQ.0 .OR. NBDCND.EQ.3) .AND. 1 (MBDCND.EQ.0 .OR. MBDCND.EQ.3)) GO TO 138 GO TO 144 C C FOR SINGULAR PROBLEMS MUST ADJUST DATA TO INSURE THAT A SOLUTION C WILL EXIST. C 138 A1 = 1.d0 A2 = 1.d0 IF (NBDCND .EQ. 3) A2 = 2.d0 IF (MBDCND .EQ. 3) A1 = 2.d0 S1 = 0.d0 MSP1 = MSTART+1 MSTM1 = MSTOP-1 NSP1 = NSTART+1 NSTM1 = NSTOP-1 DO 140 J=NSP1,NSTM1 S = 0.d0 DO 139 I=MSP1,MSTM1 S = S+F(I,J) 139 CONTINUE S1 = S1+S*A1+F(MSTART,J)+F(MSTOP,J) 140 CONTINUE S1 = A2*S1 S = 0.d0 DO 141 I=MSP1,MSTM1 S = S+F(I,NSTART)+F(I,NSTOP) 141 CONTINUE S1 = S1+S*A1+F(MSTART,NSTART)+F(MSTART,NSTOP)+F(MSTOP,NSTART)+ 1 F(MSTOP,NSTOP) S = (2.d0+dfloat(NUNK-2)*A2)*(2.d0+dfloat(MUNK-2)*A1) PERTRB = S1/S DO 143 J=NSTART,NSTOP DO 142 I=MSTART,MSTOP F(I,J) = F(I,J)-PERTRB 142 CONTINUE 143 CONTINUE PERTRB = PERTRB/DELYSQ C C SOLVE THE EQUATION. C 144 CALL GENBUN (NPEROD,NUNK,MPEROD,MUNK,W(1),W(ID2+1),W(ID3+1), 1 IDIMF,F(MSTART,NSTART),IERR1,W(ID4+1)) W(1) = W(ID4+1)+3.d0*dfloat(MUNK) C C FILL IN IDENTICAL VALUES WHEN HAVE PERIODIC BOUNDARY CONDITIONS. C IF (NBDCND .NE. 0) GO TO 146 DO 145 I=MSTART,MSTOP F(I,NP1) = F(I,1) 145 CONTINUE 146 IF (MBDCND .NE. 0) GO TO 148 DO 147 J=NSTART,NSTOP F(MP1,J) = F(1,J) 147 CONTINUE IF (NBDCND .EQ. 0) F(MP1,NP1) = F(1,NP1) 148 CONTINUE RETURN END SUBROUTINE GENBUN (NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y,IERROR,W) C C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * F I S H P A K * C * * C * * C * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF * C * * C * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS * C * * C * (VERSION 3.1 , OCTOBER 1980) * C * * C * BY * C * * C * JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET * C * * C * OF * C * * C * THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH * C * * C * BOULDER, COLORADO (80307) U.S.A. * C * * C * WHICH IS SPONSORED BY * C * * C * THE NATIONAL SCIENCE FOUNDATION * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C C * * * * * * * * * PURPOSE * * * * * * * * * * * * * * * * * * C C C SUBROUTINE GENBUN SOLVES THE LINEAR SYSTEM OF EQUATIONS C C A(I)*X(I-1,J) + B(I)*X(I,J) + C(I)*X(I+1,J) C C + X(I,J-1) - 2.*X(I,J) + X(I,J+1) = Y(I,J) C C FOR I = 1,2,...,M AND J = 1,2,...,N. C C THE INDICES I+1 AND I-1 ARE EVALUATED MODULO M, I.E., C X(0,J) = X(M,J) AND X(M+1,J) = X(1,J), AND X(I,0) MAY BE EQUAL TO C 0, X(I,2), OR X(I,N) AND X(I,N+1) MAY BE EQUAL TO 0, X(I,N-1), OR C X(I,1) DEPENDING ON AN INPUT PARAMETER. C C C * * * * * * * * PARAMETER DESCRIPTION * * * * * * * * * * C C * * * * * * ON INPUT * * * * * * C C NPEROD C INDICATES THE VALUES THAT X(I,0) AND X(I,N+1) ARE ASSUMED TO C HAVE. C C = 0 IF X(I,0) = X(I,N) AND X(I,N+1) = X(I,1). C = 1 IF X(I,0) = X(I,N+1) = 0 . C = 2 IF X(I,0) = 0 AND X(I,N+1) = X(I,N-1). C = 3 IF X(I,0) = X(I,2) AND X(I,N+1) = X(I,N-1). C = 4 IF X(I,0) = X(I,2) AND X(I,N+1) = 0. C C N C THE NUMBER OF UNKNOWNS IN THE J-DIRECTION. N MUST BE GREATER C THAN 2. C C MPEROD C = 0 IF A(1) AND C(M) ARE NOT ZERO C = 1 IF A(1) = C(M) = 0 C C M C THE NUMBER OF UNKNOWNS IN THE I-DIRECTION. M MUST BE GREATER C THAN 2. C C A,B,C C ONE-DIMENSIONAL ARRAYS OF LENGTH M THAT SPECIFY THE C COEFFICIENTS IN THE LINEAR EQUATIONS GIVEN ABOVE. IF MPEROD = 0 C THE ARRAY ELEMENTS MUST NOT DEPEND UPON THE INDEX I, BUT MUST BE C CONSTANT. SPECIFICALLY, THE SUBROUTINE CHECKS THE FOLLOWING C CONDITION C C A(I) = C(1) C C(I) = C(1) C B(I) = B(1) C C FOR I=1,2,...,M. C C IDIMY C THE ROW (OR FIRST) DIMENSION OF THE TWO-DIMENSIONAL ARRAY Y AS C IT APPEARS IN THE PROGRAM CALLING GENBUN. THIS PARAMETER IS C USED TO SPECIFY THE VARIABLE DIMENSION OF Y. IDIMY MUST BE AT C LEAST M. C C Y C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE VALUES OF THE RIGHT C SIDE OF THE LINEAR SYSTEM OF EQUATIONS GIVEN ABOVE. Y MUST BE C DIMENSIONED AT LEAST M*N. C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE USER FOR C WORK SPACE. W MAY REQUIRE UP TO 4*N + (10 + INT(LOG2(N)))*M C LOCATIONS. THE ACTUAL NUMBER OF LOCATIONS USED IS COMPUTED BY C GENBUN AND IS RETURNED IN LOCATION W(1). C C C * * * * * * ON OUTPUT * * * * * * C C Y C CONTAINS THE SOLUTION X. C C IERROR C AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS EXCEPT C FOR NUMBER ZERO, A SOLUTION IS NOT ATTEMPTED. C C = 0 NO ERROR. C = 1 M .LE. 2 . C = 2 N .LE. 2 C = 3 IDIMY .LT. M C = 4 NPEROD .LT. 0 OR NPEROD .GT. 4 C = 5 MPEROD .LT. 0 OR MPEROD .GT. 1 C = 6 A(I) .NE. C(1) OR C(I) .NE. C(1) OR B(I) .NE. B(1) FOR C SOME I=1,2,...,M. C = 7 A(1) .NE. 0 OR C(M) .NE. 0 AND MPEROD = 1 C C W C W(1) CONTAINS THE REQUIRED LENGTH OF W. C C * * * * * * * PROGRAM SPECIFICATIONS * * * * * * * * * * * * C C DIMENSION OF A(M),B(M),C(M),Y(IDIMY,N),W(SEE PARAMETER LIST) C ARGUMENTS C C LATEST JUNE 1, 1976 C REVISION C C SUBPROGRAMS GENBUN,POISD2,POISN2,POISP2,COSGEN,MERGE,TRIX,TRI3, C REQUIRED PIMACH C C SPECIAL NONE C CONDITIONS C C COMMON NONE C BLOCKS C C I/O NONE C C PRECISION SINGLE C C SPECIALIST ROLAND SWEET C C LANGUAGE FORTRAN C C HISTORY STANDARDIZED APRIL 1, 1973 C REVISED AUGUST 20,1973 C REVISED JANUARY 1, 1976 C C ALGORITHM THE LINEAR SYSTEM IS SOLVED BY A CYCLIC REDUCTION C ALGORITHM DESCRIBED IN THE REFERENCE. C C SPACE 4944(DECIMAL) = 11520(OCTAL) LOCATIONS ON THE NCAR C REQUIRED CONTROL DATA 7600 C C TIMING AND THE EXECUTION TIME T ON THE NCAR CONTROL DATA C ACCURACY 7600 FOR SUBROUTINE GENBUN IS ROUGHLY PROPORTIONAL C TO M*N*LOG2(N), BUT ALSO DEPENDS ON THE INPUT C PARAMETER NPEROD. SOME TYPICAL VALUES ARE LISTED C IN THE TABLE BELOW. MORE COMPREHENSIVE TIMING C CHARTS MAY BE FOUND IN THE REFERENCE. C TO MEASURE THE ACCURACY OF THE ALGORITHM A C UNIFORM RANDOM NUMBER GENERATOR WAS USED TO CREATE C A SOLUTION ARRAY X FOR THE SYSTEM GIVEN IN THE C 'PURPOSE' WITH C C A(I) = C(I) = -0.5*B(I) = 1, I=1,2,...,M C C AND, WHEN MPEROD = 1 C C A(1) = C(M) = 0 C A(M) = C(1) = 2. C C THE SOLUTION X WAS SUBSTITUTED INTO THE GIVEN SYS- C TEM AND, USING DOUBLE PRECISION, A RIGHT SIDE Y WAS C COMPUTED. USING THIS ARRAY Y SUBROUTINE GENBUN WAS C CALLED TO PRODUCE AN APPROXIMATE SOLUTION Z. THEN C THE RELATIVE ERROR, DEFINED AS C C E = MAX(ABS(Z(I,J)-X(I,J)))/MAX(ABS(X(I,J))) C C WHERE THE TWO MAXIMA ARE TAKEN OVER ALL I=1,2,...,M C AND J=1,2,...,N, WAS COMPUTED. THE VALUE OF E IS C GIVEN IN THE TABLE BELOW FOR SOME TYPICAL VALUES OF C M AND N. C C C M (=N) MPEROD NPEROD T(MSECS) E C ------ ------ ------ -------- ------ C C 31 0 0 36 6.E-14 C 31 1 1 21 4.E-13 C 31 1 3 41 3.E-13 C 32 0 0 29 9.E-14 C 32 1 1 32 3.E-13 C 32 1 3 48 1.E-13 C 33 0 0 36 9.E-14 C 33 1 1 30 4.E-13 C 33 1 3 34 1.E-13 C 63 0 0 150 1.E-13 C 63 1 1 91 1.E-12 C 63 1 3 173 2.E-13 C 64 0 0 122 1.E-13 C 64 1 1 128 1.E-12 C 64 1 3 199 6.E-13 C 65 0 0 143 2.E-13 C 65 1 1 120 1.E-12 C 65 1 3 138 4.E-13 C C PORTABILITY AMERICAN NATIONAL STANDARDS INSTITUE FORTRAN. C ALL MACHINE DEPENDENT CONSTANTS ARE LOCATED IN THE C FUNCTION PIMACH. C C REQUIRED COS C RESIDENT C ROUTINES C C REFERENCE SWEET, R., 'A CYCLIC REDUCTION ALGORITHM FOR C SOLVING BLOCK TRIDIAGONAL SYSTEMS OF ARBITRARY C DIMENSIONS,' SIAM J. ON NUMER. ANAL., C 14(SEPT., 1977), PP. 706-720. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C implicit double precision (a-h,o-z) DIMENSION Y(IDIMY,*) DIMENSION W(*) ,B(*) ,A(*) ,C(*) IERROR = 0 IF (M .LE. 2) IERROR = 1 IF (N .LE. 2) IERROR = 2 IF (IDIMY .LT. M) IERROR = 3 IF (NPEROD.LT.0 .OR. NPEROD.GT.4) IERROR = 4 IF (MPEROD.LT.0 .OR. MPEROD.GT.1) IERROR = 5 IF (MPEROD .EQ. 1) GO TO 102 DO 101 I=2,M IF (A(I) .NE. C(1)) GO TO 103 IF (C(I) .NE. C(1)) GO TO 103 IF (B(I) .NE. B(1)) GO TO 103 101 CONTINUE GO TO 104 102 IF (A(1).NE.0. .OR. C(M).NE.0.) IERROR = 7 GO TO 104 103 IERROR = 6 104 IF (IERROR .NE. 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 DO 106 I=1,M K = IWBA+I-1 W(K) = -A(I) K = IWBC+I-1 W(K) = -C(I) K = IWBB+I-1 W(K) = 2.d0-B(I) DO 105 J=1,N Y(I,J) = -Y(I,J) 105 CONTINUE 106 CONTINUE MP = MPEROD+1 NP = NPEROD+1 GO TO (114,107),MP 107 GO TO (108,109,110,111,123),NP 108 CALL POISP2 (M,N,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2), 1 W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS), 2 W(IWP)) GO TO 112 109 CALL POISD2 (M,N,1,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWW1), 1 W(IWD),W(IWTCOS),W(IWP)) GO TO 112 110 CALL POISN2 (M,N,1,2,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2), 1 W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS), 2 W(IWP)) GO TO 112 111 CALL POISN2 (M,N,1,1,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2), 1 W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS), 2 W(IWP)) 112 IPSTOR = W(IWW1) IREV = 2 IF (NPEROD .EQ. 4) GO TO 124 113 GO TO (127,133),MP 114 CONTINUE C C REORDER UNKNOWNS WHEN MP =0 C MH = (M+1)/2 MHM1 = MH-1 MODD = 1 IF (MH*2 .EQ. M) MODD = 2 DO 119 J=1,N DO 115 I=1,MHM1 MHPI = MH+I MHMI = MH-I W(I) = Y(MHMI,J)-Y(MHPI,J) W(MHPI) = Y(MHMI,J)+Y(MHPI,J) 115 CONTINUE W(MH) = 2.d0*Y(MH,J) GO TO (117,116),MODD 116 W(M) = 2.d0*Y(M,J) 117 CONTINUE DO 118 I=1,M Y(I,J) = W(I) 118 CONTINUE 119 CONTINUE K = IWBC+MHM1-1 I = IWBA+MHM1 W(K) = 0.d0 W(I) = 0.d0 W(K+1) = 2.d0*W(K+1) GO TO (120,121),MODD 120 CONTINUE K = IWBB+MHM1-1 W(K) = W(K)-W(I-1) W(IWBC-1) = W(IWBC-1)+W(IWBB-1) GO TO 122 121 W(IWBB-1) = W(K+1) 122 CONTINUE GO TO 107 C C REVERSE COLUMNS WHEN NPEROD = 4. C 123 IREV = 1 NBY2 = N/2 124 DO 126 J=1,NBY2 MSKIP = N+1-J DO 125 I=1,M A1 = Y(I,J) Y(I,J) = Y(I,MSKIP) Y(I,MSKIP) = A1 125 CONTINUE 126 CONTINUE GO TO (110,113),IREV 127 CONTINUE DO 132 J=1,N DO 128 I=1,MHM1 MHMI = MH-I MHPI = MH+I W(MHMI) = .5d0*(Y(MHPI,J)+Y(I,J)) W(MHPI) = .5d0*(Y(MHPI,J)-Y(I,J)) 128 CONTINUE W(MH) = .5d0*Y(MH,J) GO TO (130,129),MODD 129 W(M) = .5d0*Y(M,J) 130 CONTINUE DO 131 I=1,M Y(I,J) = W(I) 131 CONTINUE 132 CONTINUE 133 CONTINUE C C RETURN STORAGE REQUIREMENTS FOR W ARRAY. C W(1) = IPSTOR+IWP-1 RETURN END SUBROUTINE POISP2 (M,N,A,BB,C,Q,IDIMQ,B,B2,B3,W,W2,W3,D,TCOS,P) C C SUBROUTINE TO SOLVE POISSON EQUATION WITH PERIODIC BOUNDARY C CONDITIONS. C implicit double precision (a-h,o-z) DIMENSION A(*) ,BB(*) ,C(*) ,Q(IDIMQ,*) , 1 B(*) ,B2(*) ,B3(*) ,W(*) , 2 W2(*) ,W3(*) ,D(*) ,TCOS(*) , 3 P(*) MR = M NR = (N+1)/2 NRM1 = NR-1 IF (2*NR .NE. N) GO TO 107 C C EVEN NUMBER OF UNKNOWNS C DO 102 J=1,NRM1 NRMJ = NR-J NRPJ = NR+J DO 101 I=1,MR S = Q(I,NRMJ)-Q(I,NRPJ) T = Q(I,NRMJ)+Q(I,NRPJ) Q(I,NRMJ) = S Q(I,NRPJ) = T 101 CONTINUE 102 CONTINUE DO 103 I=1,MR Q(I,NR) = 2.d0*Q(I,NR) Q(I,N) = 2.d0*Q(I,N) 103 CONTINUE CALL POISD2 (MR,NRM1,1,A,BB,C,Q,IDIMQ,B,W,D,TCOS,P) IPSTOR = W(1) CALL POISN2 (MR,NR+1,1,1,A,BB,C,Q(1,NR),IDIMQ,B,B2,B3,W,W2,W3,D, 1 TCOS,P) IPSTOR = MAX0(IPSTOR,INT(W(1))) DO 105 J=1,NRM1 NRMJ = NR-J NRPJ = NR+J DO 104 I=1,MR S = .5d0*(Q(I,NRPJ)+Q(I,NRMJ)) T = .5d0*(Q(I,NRPJ)-Q(I,NRMJ)) Q(I,NRMJ) = S Q(I,NRPJ) = T 104 CONTINUE 105 CONTINUE DO 106 I=1,MR Q(I,NR) = .5d0*Q(I,NR) Q(I,N) = .5d0*Q(I,N) 106 CONTINUE GO TO 118 107 CONTINUE C C ODD NUMBER OF UNKNOWNS C DO 109 J=1,NRM1 NRPJ = N+1-J DO 108 I=1,MR S = Q(I,J)-Q(I,NRPJ) T = Q(I,J)+Q(I,NRPJ) Q(I,J) = S Q(I,NRPJ) = T 108 CONTINUE 109 CONTINUE DO 110 I=1,MR Q(I,NR) = 2.d0*Q(I,NR) 110 CONTINUE LH = NRM1/2 DO 112 J=1,LH NRMJ = NR-J DO 111 I=1,MR S = Q(I,J) Q(I,J) = Q(I,NRMJ) Q(I,NRMJ) = S 111 CONTINUE 112 CONTINUE CALL POISD2 (MR,NRM1,2,A,BB,C,Q,IDIMQ,B,W,D,TCOS,P) IPSTOR = W(1) CALL POISN2 (MR,NR,2,1,A,BB,C,Q(1,NR),IDIMQ,B,B2,B3,W,W2,W3,D, 1 TCOS,P) IPSTOR = MAX0(IPSTOR,INT(W(1))) DO 114 J=1,NRM1 NRPJ = NR+J DO 113 I=1,MR S = .5d0*(Q(I,NRPJ)+Q(I,J)) T = .5d0*(Q(I,NRPJ)-Q(I,J)) Q(I,NRPJ) = T Q(I,J) = S 113 CONTINUE 114 CONTINUE DO 115 I=1,MR Q(I,NR) = .5d0*Q(I,NR) 115 CONTINUE DO 117 J=1,LH NRMJ = NR-J DO 116 I=1,MR S = Q(I,J) Q(I,J) = Q(I,NRMJ) Q(I,NRMJ) = S 116 CONTINUE 117 CONTINUE 118 CONTINUE C C RETURN STORAGE REQUIREMENTS FOR P VECTORS. C W(1) = IPSTOR RETURN END SUBROUTINE POISD2 (MR,NR,ISTAG,BA,BB,BC,Q,IDIMQ,B,W,D,TCOS,P) C C SUBROUTINE TO SOLVE POISSON'S EQUATION FOR DIRICHLET BOUNDARY C CONDITIONS. C C ISTAG = 1 IF THE LAST DIAGONAL BLOCK IS THE MATRIX A. C ISTAG = 2 IF THE LAST DIAGONAL BLOCK IS THE MATRIX A+I. C implicit double precision (a-h,o-z) DIMENSION Q(IDIMQ,*) ,BA(*) ,BB(*) ,BC(*) , 1 TCOS(*) ,B(*) ,D(*) ,W(*) , 2 P(*) M = MR N = NR JSH = 0 FI = 1.d0/dfloat(ISTAG) IP = -M IPSTOR = 0 GO TO (101,102),ISTAG 101 KR = 0 IRREG = 1 IF (N .GT. 1) GO TO 106 TCOS(1) = 0.d0 GO TO 103 102 KR = 1 JSTSAV = 1 IRREG = 2 IF (N .GT. 1) GO TO 106 TCOS(1) = -1.d0 103 DO 104 I=1,M B(I) = Q(I,1) 104 CONTINUE CALL TRIX (1,0,M,BA,BB,BC,B,TCOS,D,W) DO 105 I=1,M Q(I,1) = B(I) 105 CONTINUE GO TO 183 106 LR = 0 DO 107 I=1,M P(I) = 0.d0 107 CONTINUE NUN = N JST = 1 JSP = N C C IRREG = 1 WHEN NO IRREGULARITIES HAVE OCCURRED, OTHERWISE IT IS 2. C 108 L = 2*JST NODD = 2-2*((NUN+1)/2)+NUN C C NODD = 1 WHEN NUN IS ODD, OTHERWISE IT IS 2. C GO TO (110,109),NODD 109 JSP = JSP-L GO TO 111 110 JSP = JSP-JST IF (IRREG .NE. 1) JSP = JSP-L 111 CONTINUE C C REGULAR REDUCTION C CALL COSGEN (JST,1,0.5d0,0.0d0,TCOS) IF (L .GT. JSP) GO TO 118 DO 117 J=L,JSP,L JM1 = J-JSH JP1 = J+JSH JM2 = J-JST JP2 = J+JST JM3 = JM2-JSH JP3 = JP2+JSH IF (JST .NE. 1) GO TO 113 DO 112 I=1,M B(I) = 2.d0*Q(I,J) Q(I,J) = Q(I,JM2)+Q(I,JP2) 112 CONTINUE GO TO 115 113 DO 114 I=1,M T = Q(I,J)-Q(I,JM1)-Q(I,JP1)+Q(I,JM2)+Q(I,JP2) B(I) = T+Q(I,J)-Q(I,JM3)-Q(I,JP3) Q(I,J) = T 114 CONTINUE 115 CONTINUE CALL TRIX (JST,0,M,BA,BB,BC,B,TCOS,D,W) DO 116 I=1,M Q(I,J) = Q(I,J)+B(I) 116 CONTINUE 117 CONTINUE C C REDUCTION FOR LAST UNKNOWN C 118 GO TO (119,136),NODD 119 GO TO (152,120),IRREG C C ODD NUMBER OF UNKNOWNS C 120 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 .NE. 1) GO TO 123 DO 122 I=1,M B(I) = Q(I,J) Q(I,J) = 0.d0 122 CONTINUE GO TO 130 123 GO TO (124,126),NODDPR 124 DO 125 I=1,M IP1 = IP+I B(I) = .5d0*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))+P(IP1)+Q(I,J) 125 CONTINUE GO TO 128 126 DO 127 I=1,M B(I) = .5d0*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))+Q(I,JP2)-Q(I,JP1)+Q(I,J) 127 CONTINUE 128 DO 129 I=1,M Q(I,J) = .5d0*(Q(I,J)-Q(I,JM1)-Q(I,JP1)) 129 CONTINUE 130 CALL TRIX (JST,0,M,BA,BB,BC,B,TCOS,D,W) IP = IP+M IPSTOR = MAX0(IPSTOR,IP+M) DO 131 I=1,M IP1 = IP+I P(IP1) = Q(I,J)+B(I) B(I) = Q(I,JP2)+P(IP1) 131 CONTINUE IF (LR .NE. 0) GO TO 133 DO 132 I=1,JST KRPI = KR+I TCOS(KRPI) = TCOS(I) 132 CONTINUE GO TO 134 133 CONTINUE CALL COSGEN (LR,JSTSAV,0.d0,FI,TCOS(JST+1)) CALL MERGE (TCOS,0,JST,JST,LR,KR) 134 CONTINUE CALL COSGEN (KR,JSTSAV,0.0d0,FI,TCOS) CALL TRIX (KR,KR,M,BA,BB,BC,B,TCOS,D,W) DO 135 I=1,M IP1 = IP+I Q(I,J) = Q(I,JM2)+B(I)+P(IP1) 135 CONTINUE LR = KR KR = KR+L GO TO 152 C C EVEN NUMBER OF UNKNOWNS C 136 JSP = JSP+L J = JSP JM1 = J-JSH JP1 = J+JSH JM2 = J-JST JP2 = J+JST JM3 = JM2-JSH GO TO (137,138),IRREG 137 CONTINUE JSTSAV = JST IDEG = JST KR = L GO TO 139 138 CALL COSGEN (KR,JSTSAV,0.0d0,FI,TCOS) CALL COSGEN (LR,JSTSAV,0.0d0,FI,TCOS(KR+1)) IDEG = KR KR = KR+JST 139 IF (JST .NE. 1) GO TO 141 IRREG = 2 DO 140 I=1,M B(I) = Q(I,J) Q(I,J) = Q(I,JM2) 140 CONTINUE GO TO 150 141 DO 142 I=1,M B(I) = Q(I,J)+.5d0*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3)) 142 CONTINUE GO TO (143,145),IRREG 143 DO 144 I=1,M Q(I,J) = Q(I,JM2)+.5d0*(Q(I,J)-Q(I,JM1)-Q(I,JP1)) 144 CONTINUE IRREG = 2 GO TO 150 145 CONTINUE GO TO (146,148),NODDPR 146 DO 147 I=1,M IP1 = IP+I Q(I,J) = Q(I,JM2)+P(IP1) 147 CONTINUE IP = IP-M GO TO 150 148 DO 149 I=1,M Q(I,J) = Q(I,JM2)+Q(I,J)-Q(I,JM1) 149 CONTINUE 150 CALL TRIX (IDEG,LR,M,BA,BB,BC,B,TCOS,D,W) DO 151 I=1,M Q(I,J) = Q(I,J)+B(I) 151 CONTINUE 152 NUN = NUN/2 NODDPR = NODD JSH = JST JST = 2*JST IF (NUN .GE. 2) GO TO 108 C C START SOLUTION. C J = JSP DO 153 I=1,M B(I) = Q(I,J) 153 CONTINUE GO TO (154,155),IRREG 154 CONTINUE CALL COSGEN (JST,1,0.5d0,0.0d0,TCOS) IDEG = JST GO TO 156 155 KR = LR+JST CALL COSGEN (KR,JSTSAV,0.0d0,FI,TCOS) CALL COSGEN (LR,JSTSAV,0.0d0,FI,TCOS(KR+1)) IDEG = KR 156 CONTINUE CALL TRIX (IDEG,LR,M,BA,BB,BC,B,TCOS,D,W) JM1 = J-JSH JP1 = J+JSH GO TO (157,159),IRREG 157 DO 158 I=1,M Q(I,J) = .5d0*(Q(I,J)-Q(I,JM1)-Q(I,JP1))+B(I) 158 CONTINUE GO TO 164 159 GO TO (160,162),NODDPR 160 DO 161 I=1,M IP1 = IP+I Q(I,J) = P(IP1)+B(I) 161 CONTINUE IP = IP-M GO TO 164 162 DO 163 I=1,M Q(I,J) = Q(I,J)-Q(I,JM1)+B(I) 163 CONTINUE 164 CONTINUE C C START BACK SUBSTITUTION. C JST = JST/2 JSH = JST/2 NUN = 2*NUN IF (NUN .GT. N) GO TO 183 DO 182 J=JST,N,L JM1 = J-JSH JP1 = J+JSH JM2 = J-JST JP2 = J+JST IF (J .GT. JST) GO TO 166 DO 165 I=1,M B(I) = Q(I,J)+Q(I,JP2) 165 CONTINUE GO TO 170 166 IF (JP2 .LE. N) GO TO 168 DO 167 I=1,M B(I) = Q(I,J)+Q(I,JM2) 167 CONTINUE IF (JST .LT. JSTSAV) IRREG = 1 GO TO (170,171),IRREG 168 DO 169 I=1,M B(I) = Q(I,J)+Q(I,JM2)+Q(I,JP2) 169 CONTINUE 170 CONTINUE CALL COSGEN (JST,1,0.5d0,0.0d0,TCOS) IDEG = JST JDEG = 0 GO TO 172 171 IF (J+L .GT. N) LR = LR-JST KR = JST+LR CALL COSGEN (KR,JSTSAV,0.0d0,FI,TCOS) CALL COSGEN (LR,JSTSAV,0.0d0,FI,TCOS(KR+1)) IDEG = KR JDEG = LR 172 CONTINUE CALL TRIX (IDEG,JDEG,M,BA,BB,BC,B,TCOS,D,W) IF (JST .GT. 1) GO TO 174 DO 173 I=1,M Q(I,J) = B(I) 173 CONTINUE GO TO 182 174 IF (JP2 .GT. N) GO TO 177 175 DO 176 I=1,M Q(I,J) = .5d0*(Q(I,J)-Q(I,JM1)-Q(I,JP1))+B(I) 176 CONTINUE GO TO 182 177 GO TO (175,178),IRREG 178 IF (J+JSH .GT. N) GO TO 180 DO 179 I=1,M IP1 = IP+I Q(I,J) = B(I)+P(IP1) 179 CONTINUE IP = IP-M GO TO 182 180 DO 181 I=1,M Q(I,J) = B(I)+Q(I,J)-Q(I,JM1) 181 CONTINUE 182 CONTINUE L = L/2 GO TO 164 183 CONTINUE C C RETURN STORAGE REQUIREMENTS FOR P VECTORS. C W(1) = IPSTOR RETURN END SUBROUTINE POISN2 (M,N,ISTAG,MIXBND,A,BB,C,Q,IDIMQ,B,B2,B3,W,W2, 1 W3,D,TCOS,P) C C SUBROUTINE TO SOLVE POISSON'S EQUATION WITH NEUMANN BOUNDARY C CONDITIONS. C C ISTAG = 1 IF THE LAST DIAGONAL BLOCK IS A. C ISTAG = 2 IF THE LAST DIAGONAL BLOCK IS A-I. C MIXBND = 1 IF HAVE NEUMANN BOUNDARY CONDITIONS AT BOTH BOUNDARIES. C MIXBND = 2 IF HAVE NEUMANN BOUNDARY CONDITIONS AT BOTTOM AND C DIRICHLET CONDITION AT TOP. (FOR THIS CASE, MUST HAVE ISTAG = 1.) C implicit double precision (a-h,o-z) DIMENSION A(*) ,BB(*) ,C(*) ,Q(IDIMQ,*) , 1 B(*) ,B2(*) ,B3(*) ,W(*) , 2 W2(*) ,W3(*) ,D(*) ,TCOS(*) , 3 K(4) ,P(*) EQUIVALENCE (K(1),K1) ,(K(2),K2) ,(K(3),K3) ,(K(4),K4) FISTAG = 3-ISTAG FNUM = 1.d0/dfloat(ISTAG) FDEN = 0.5d0*dfloat(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 DO 102 I=1,MR Q(I,N) = .5d0*Q(I,N) 102 CONTINUE GO TO (103,104),MIXBND 103 IF (N .LE. 3) GO TO 155 104 CONTINUE JR = 2*I2R NROD = 1 IF ((NR/2)*2 .EQ. NR) NROD = 0 GO TO (105,106),MIXBND 105 JSTART = 1 GO TO 107 106 JSTART = JR NROD = 1-NROD 107 CONTINUE JSTOP = NLAST-JR IF (NROD .EQ. 0) JSTOP = JSTOP-I2R CALL COSGEN (I2R,1,0.5d0,0.0d0,TCOS) I2RBY2 = I2R/2 IF (JSTOP .GE. JSTART) GO TO 108 J = JR GO TO 116 108 CONTINUE C C REGULAR REDUCTION. C DO 115 J=JSTART,JSTOP,JR JP1 = J+I2RBY2 JP2 = J+I2R JP3 = JP2+I2RBY2 JM1 = J-I2RBY2 JM2 = J-I2R JM3 = JM2-I2RBY2 IF (J .NE. 1) GO TO 109 JM1 = JP1 JM2 = JP2 JM3 = JP3 109 CONTINUE IF (I2R .NE. 1) GO TO 111 IF (J .EQ. 1) JM2 = JP2 DO 110 I=1,MR B(I) = 2.d0*Q(I,J) Q(I,J) = Q(I,JM2)+Q(I,JP2) 110 CONTINUE GO TO 113 111 CONTINUE DO 112 I=1,MR FI = Q(I,J) Q(I,J) = Q(I,J)-Q(I,JM1)-Q(I,JP1)+Q(I,JM2)+Q(I,JP2) B(I) = FI+Q(I,J)-Q(I,JM3)-Q(I,JP3) 112 CONTINUE 113 CONTINUE CALL TRIX (I2R,0,MR,A,BB,C,B,TCOS,D,W) DO 114 I=1,MR Q(I,J) = Q(I,J)+B(I) 114 CONTINUE C C END OF REDUCTION FOR REGULAR UNKNOWNS. C 115 CONTINUE C C BEGIN SPECIAL REDUCTION FOR LAST UNKNOWN. C J = JSTOP+JR 116 NLAST = J JM1 = J-I2RBY2 JM2 = J-I2R JM3 = JM2-I2RBY2 IF (NROD .EQ. 0) GO TO 128 C C ODD NUMBER OF UNKNOWNS C IF (I2R .NE. 1) GO TO 118 DO 117 I=1,MR B(I) = FISTAG*Q(I,J) Q(I,J) = Q(I,JM2) 117 CONTINUE GO TO 126 118 DO 119 I=1,MR B(I) = Q(I,J)+.5d0*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3)) 119 CONTINUE IF (NRODPR .NE. 0) GO TO 121 DO 120 I=1,MR II = IP+I Q(I,J) = Q(I,JM2)+P(II) 120 CONTINUE IP = IP-MR GO TO 123 121 CONTINUE DO 122 I=1,MR Q(I,J) = Q(I,J)-Q(I,JM1)+Q(I,JM2) 122 CONTINUE 123 IF (LR .EQ. 0) GO TO 124 CALL COSGEN (LR,1,0.5d0,FDEN,TCOS(KR+1)) GO TO 126 124 CONTINUE DO 125 I=1,MR B(I) = FISTAG*B(I) 125 CONTINUE 126 CONTINUE CALL COSGEN (KR,1,0.5d0,FDEN,TCOS) CALL TRIX (KR,LR,MR,A,BB,C,B,TCOS,D,W) DO 127 I=1,MR Q(I,J) = Q(I,J)+B(I) 127 CONTINUE KR = KR+I2R GO TO 151 128 CONTINUE C C EVEN NUMBER OF UNKNOWNS C JP1 = J+I2RBY2 JP2 = J+I2R IF (I2R .NE. 1) GO TO 135 DO 129 I=1,MR B(I) = Q(I,J) 129 CONTINUE CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W) IP = 0 IPSTOR = MR GO TO (133,130),ISTAG 130 DO 131 I=1,MR P(I) = B(I) B(I) = B(I)+Q(I,N) 131 CONTINUE TCOS(1) = 1.d0 TCOS(2) = 0.d0 CALL TRIX (1,1,MR,A,BB,C,B,TCOS,D,W) DO 132 I=1,MR Q(I,J) = Q(I,JM2)+P(I)+B(I) 132 CONTINUE GO TO 150 133 CONTINUE DO 134 I=1,MR P(I) = B(I) Q(I,J) = Q(I,JM2)+2.d0*Q(I,JP2)+3.d0*B(I) 134 CONTINUE GO TO 150 135 CONTINUE DO 136 I=1,MR B(I) = Q(I,J)+.5d0*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3)) 136 CONTINUE IF (NRODPR .NE. 0) GO TO 138 DO 137 I=1,MR II = IP+I B(I) = B(I)+P(II) 137 CONTINUE GO TO 140 138 CONTINUE DO 139 I=1,MR B(I) = B(I)+Q(I,JP2)-Q(I,JP1) 139 CONTINUE 140 CONTINUE CALL TRIX (I2R,0,MR,A,BB,C,B,TCOS,D,W) IP = IP+MR IPSTOR = MAX0(IPSTOR,IP+MR) DO 141 I=1,MR II = IP+I P(II) = B(I)+.5d0*(Q(I,J)-Q(I,JM1)-Q(I,JP1)) B(I) = P(II)+Q(I,JP2) 141 CONTINUE IF (LR .EQ. 0) GO TO 142 CALL COSGEN (LR,1,0.5d0,FDEN,TCOS(I2R+1)) CALL MERGE (TCOS,0,I2R,I2R,LR,KR) GO TO 144 142 DO 143 I=1,I2R II = KR+I TCOS(II) = TCOS(I) 143 CONTINUE 144 CALL COSGEN (KR,1,0.5d0,FDEN,TCOS) IF (LR .NE. 0) GO TO 145 GO TO (146,145),ISTAG 145 CONTINUE CALL TRIX (KR,KR,MR,A,BB,C,B,TCOS,D,W) GO TO 148 146 CONTINUE DO 147 I=1,MR B(I) = FISTAG*B(I) 147 CONTINUE 148 CONTINUE DO 149 I=1,MR II = IP+I Q(I,J) = Q(I,JM2)+P(II)+B(I) 149 CONTINUE 150 CONTINUE LR = KR KR = KR+JR 151 CONTINUE GO TO (152,153),MIXBND 152 NR = (NLAST-1)/JR+1 IF (NR .LE. 3) GO TO 155 GO TO 154 153 NR = NLAST/JR IF (NR .LE. 1) GO TO 192 154 I2R = JR NRODPR = NROD GO TO 104 155 CONTINUE C C BEGIN SOLUTION C J = 1+JR JM1 = J-I2R JP1 = J+I2R JM2 = NLAST-I2R IF (NR .EQ. 2) GO TO 184 IF (LR .NE. 0) GO TO 170 IF (N .NE. 3) GO TO 161 C C CASE N = 3. C GO TO (156,168),ISTAG 156 CONTINUE DO 157 I=1,MR B(I) = Q(I,2) 157 CONTINUE TCOS(1) = 0.d0 CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W) DO 158 I=1,MR Q(I,2) = B(I) B(I) = 4.d0*B(I)+Q(I,1)+2.d0*Q(I,3) 158 CONTINUE TCOS(1) = -2.d0 TCOS(2) = 2.d0 I1 = 2 I2 = 0 CALL TRIX (I1,I2,MR,A,BB,C,B,TCOS,D,W) DO 159 I=1,MR Q(I,2) = Q(I,2)+B(I) B(I) = Q(I,1)+2.d0*Q(I,2) 159 CONTINUE TCOS(1) = 0.d0 CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W) DO 160 I=1,MR Q(I,1) = B(I) 160 CONTINUE JR = 1 I2R = 0 GO TO 194 C C CASE N = 2**P+1 C 161 CONTINUE GO TO (162,170),ISTAG 162 CONTINUE DO 163 I=1,MR B(I) = Q(I,J)+.5d0*Q(I,1)-Q(I,JM1)+Q(I,NLAST)-Q(I,JM2) 163 CONTINUE CALL COSGEN (JR,1,0.5d0,0.0d0,TCOS) CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W) DO 164 I=1,MR Q(I,J) = .5d0*(Q(I,J)-Q(I,JM1)-Q(I,JP1))+B(I) B(I) = Q(I,1)+2.d0*Q(I,NLAST)+4.d0*Q(I,J) 164 CONTINUE JR2 = 2*JR CALL COSGEN (JR,1,0.0d0,0.0d0,TCOS) DO 165 I=1,JR I1 = JR+I I2 = JR+1-I TCOS(I1) = -TCOS(I2) 165 CONTINUE CALL TRIX (JR2,0,MR,A,BB,C,B,TCOS,D,W) DO 166 I=1,MR Q(I,J) = Q(I,J)+B(I) B(I) = Q(I,1)+2.d0*Q(I,J) 166 CONTINUE CALL COSGEN (JR,1,0.5d0,0.0d0,TCOS) CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W) DO 167 I=1,MR Q(I,1) = .5d0*Q(I,1)-Q(I,JM1)+B(I) 167 CONTINUE GO TO 194 C C CASE OF GENERAL N WITH NR = 3 . C 168 DO 169 I=1,MR B(I) = Q(I,2) Q(I,2) = 0.d0 B2(I) = Q(I,3) B3(I) = Q(I,1) 169 CONTINUE JR = 1 I2R = 0 J = 2 GO TO 177 170 CONTINUE DO 171 I=1,MR B(I) = .5d0*Q(I,1)-Q(I,JM1)+Q(I,J) 171 CONTINUE IF (NROD .NE. 0) GO TO 173 DO 172 I=1,MR II = IP+I B(I) = B(I)+P(II) 172 CONTINUE GO TO 175 173 DO 174 I=1,MR B(I) = B(I)+Q(I,NLAST)-Q(I,JM2) 174 CONTINUE 175 CONTINUE DO 176 I=1,MR T = .5d0*(Q(I,J)-Q(I,JM1)-Q(I,JP1)) Q(I,J) = T B2(I) = Q(I,NLAST)+T B3(I) = Q(I,1)+2.d0*T 176 CONTINUE 177 CONTINUE K1 = KR+2*JR-1 K2 = KR+JR TCOS(K1+1) = -2.d0 K4 = K1+3-ISTAG CALL COSGEN (K2+ISTAG-2,1,0.0d0,FNUM,TCOS(K4)) K4 = K1+K2+1 CALL COSGEN (JR-1,1,0.0d0,1.0d0,TCOS(K4)) CALL MERGE (TCOS,K1,K2,K1+K2,JR-1,0) K3 = K1+K2+LR CALL COSGEN (JR,1,0.5d0,0.0d0,TCOS(K3+1)) K4 = K3+JR+1 CALL COSGEN (KR,1,0.5d0,FDEN,TCOS(K4)) CALL MERGE (TCOS,K3,JR,K3+JR,KR,K1) IF (LR .EQ. 0) GO TO 178 CALL COSGEN (LR,1,0.5d0,FDEN,TCOS(K4)) CALL MERGE (TCOS,K3,JR,K3+JR,LR,K3-LR) CALL COSGEN (KR,1,0.5d0,FDEN,TCOS(K4)) 178 K3 = KR K4 = KR CALL TRI3 (MR,A,BB,C,K,B,B2,B3,TCOS,D,W,W2,W3) DO 179 I=1,MR B(I) = B(I)+B2(I)+B3(I) 179 CONTINUE TCOS(1) = 2.d0 CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W) DO 180 I=1,MR Q(I,J) = Q(I,J)+B(I) B(I) = Q(I,1)+2.d0*Q(I,J) 180 CONTINUE CALL COSGEN (JR,1,0.5d0,0.0d0,TCOS) CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W) IF (JR .NE. 1) GO TO 182 DO 181 I=1,MR Q(I,1) = B(I) 181 CONTINUE GO TO 194 182 CONTINUE DO 183 I=1,MR Q(I,1) = .5d0*Q(I,1)-Q(I,JM1)+B(I) 183 CONTINUE GO TO 194 184 CONTINUE IF (N .NE. 2) GO TO 188 C C CASE N = 2 C DO 185 I=1,MR B(I) = Q(I,1) 185 CONTINUE TCOS(1) = 0.d0 CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W) DO 186 I=1,MR Q(I,1) = B(I) B(I) = 2.d0*(Q(I,2)+B(I))*FISTAG 186 CONTINUE TCOS(1) = -FISTAG TCOS(2) = 2.d0 CALL TRIX (2,0,MR,A,BB,C,B,TCOS,D,W) DO 187 I=1,MR Q(I,1) = Q(I,1)+B(I) 187 CONTINUE JR = 1 I2R = 0 GO TO 194 188 CONTINUE C C CASE OF GENERAL N AND NR = 2 . C DO 189 I=1,MR II = IP+I B3(I) = 0.d0 B(I) = Q(I,1)+2.d0*P(II) Q(I,1) = .5d0*Q(I,1)-Q(I,JM1) B2(I) = 2.d0*(Q(I,1)+Q(I,NLAST)) 189 CONTINUE K1 = KR+JR-1 TCOS(K1+1) = -2.d0 K4 = K1+3-ISTAG CALL COSGEN (KR+ISTAG-2,1,0.0d0,FNUM,TCOS(K4)) K4 = K1+KR+1 CALL COSGEN (JR-1,1,0.0d0,1.0d0,TCOS(K4)) CALL MERGE (TCOS,K1,KR,K1+KR,JR-1,0) CALL COSGEN (KR,1,0.5d0,FDEN,TCOS(K1+1)) K2 = KR K4 = K1+K2+1 CALL COSGEN (LR,1,0.5d0,FDEN,TCOS(K4)) K3 = LR K4 = 0 CALL TRI3 (MR,A,BB,C,K,B,B2,B3,TCOS,D,W,W2,W3) DO 190 I=1,MR B(I) = B(I)+B2(I) 190 CONTINUE TCOS(1) = 2.d0 CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W) DO 191 I=1,MR Q(I,1) = Q(I,1)+B(I) 191 CONTINUE GO TO 194 192 DO 193 I=1,MR B(I) = Q(I,NLAST) 193 CONTINUE GO TO 196 194 CONTINUE C C START BACK SUBSTITUTION. C J = NLAST-JR DO 195 I=1,MR B(I) = Q(I,NLAST)+Q(I,J) 195 CONTINUE 196 JM2 = NLAST-I2R IF (JR .NE. 1) GO TO 198 DO 197 I=1,MR Q(I,NLAST) = 0.d0 197 CONTINUE GO TO 202 198 CONTINUE IF (NROD .NE. 0) GO TO 200 DO 199 I=1,MR II = IP+I Q(I,NLAST) = P(II) 199 CONTINUE IP = IP-MR GO TO 202 200 DO 201 I=1,MR Q(I,NLAST) = Q(I,NLAST)-Q(I,JM2) 201 CONTINUE 202 CONTINUE CALL COSGEN (KR,1,0.5d0,FDEN,TCOS) CALL COSGEN (LR,1,0.5d0,FDEN,TCOS(KR+1)) IF (LR .NE. 0) GO TO 204 DO 203 I=1,MR B(I) = FISTAG*B(I) 203 CONTINUE 204 CONTINUE CALL TRIX (KR,LR,MR,A,BB,C,B,TCOS,D,W) DO 205 I=1,MR Q(I,NLAST) = Q(I,NLAST)+B(I) 205 CONTINUE NLASTP = NLAST 206 CONTINUE JSTEP = JR JR = I2R I2R = I2R/2 IF (JR .EQ. 0) GO TO 222 GO TO (207,208),MIXBND 207 JSTART = 1+JR GO TO 209 208 JSTART = JR 209 CONTINUE KR = KR-JR IF (NLAST+JR .GT. N) GO TO 210 KR = KR-JR NLAST = NLAST+JR JSTOP = NLAST-JSTEP GO TO 211 210 CONTINUE JSTOP = NLAST-JR 211 CONTINUE LR = KR-JR CALL COSGEN (JR,1,0.5d0,0.0d0,TCOS) DO 221 J=JSTART,JSTOP,JSTEP JM2 = J-JR JP2 = J+JR IF (J .NE. JR) GO TO 213 DO 212 I=1,MR B(I) = Q(I,J)+Q(I,JP2) 212 CONTINUE GO TO 215 213 CONTINUE DO 214 I=1,MR B(I) = Q(I,J)+Q(I,JM2)+Q(I,JP2) 214 CONTINUE 215 CONTINUE IF (JR .NE. 1) GO TO 217 DO 216 I=1,MR Q(I,J) = 0.d0 216 CONTINUE GO TO 219 217 CONTINUE JM1 = J-I2R JP1 = J+I2R DO 218 I=1,MR Q(I,J) = .5d0*(Q(I,J)-Q(I,JM1)-Q(I,JP1)) 218 CONTINUE 219 CONTINUE CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W) DO 220 I=1,MR Q(I,J) = Q(I,J)+B(I) 220 CONTINUE 221 CONTINUE NROD = 1 IF (NLAST+I2R .LE. N) NROD = 0 IF (NLASTP .NE. NLAST) GO TO 194 GO TO 206 222 CONTINUE C C RETURN STORAGE REQUIREMENTS FOR P VECTORS. C W(1) = IPSTOR RETURN END SUBROUTINE TRIX (IDEGBR,IDEGCR,M,A,B,C,Y,TCOS,D,W) C C SUBROUTINE TO SOLVE A SYSTEM OF LINEAR EQUATIONS WHERE THE C COEFFICIENT MATRIX IS A RATIONAL FUNCTION IN THE MATRIX GIVEN BY C TRIDIAGONAL ( . . . , A(I), B(I), C(I), . . . ). C implicit double precision (a-h,o-z) DIMENSION A(*) ,B(*) ,C(*) ,Y(*) , 1 TCOS(*) ,D(*) ,W(*) MM1 = M-1 FB = IDEGBR+1 FC = IDEGCR+1 L = FB/FC LINT = 1 DO 108 K=1,IDEGBR X = TCOS(K) IF (K .NE. L) GO TO 102 I = IDEGBR+LINT XX = X-TCOS(I) DO 101 I=1,M W(I) = Y(I) Y(I) = XX*Y(I) 101 CONTINUE 102 CONTINUE Z = 1.d0/(B(1)-X) D(1) = C(1)*Z Y(1) = Y(1)*Z DO 103 I=2,MM1 Z = 1.d0/(B(I)-X-A(I)*D(I-1)) D(I) = C(I)*Z Y(I) = (Y(I)-A(I)*Y(I-1))*Z 103 CONTINUE Z = B(M)-X-A(M)*D(MM1) IF (Z .NE. 0.d0) GO TO 104 Y(M) = 0.d0 GO TO 105 104 Y(M) = (Y(M)-A(M)*Y(MM1))/Z 105 CONTINUE DO 106 IP=1,MM1 I = M-IP Y(I) = Y(I)-D(I)*Y(I+1) 106 CONTINUE IF (K .NE. L) GO TO 108 DO 107 I=1,M Y(I) = Y(I)+W(I) 107 CONTINUE LINT = LINT+1 L = (dfloat(LINT)*FB)/FC 108 CONTINUE RETURN END SUBROUTINE COSGEN (N,IJUMP,FNUM,FDEN,A) implicit double precision (a-h,o-z) DIMENSION A(*) C C C THIS SUBROUTINE COMPUTES REQUIRED COSINE VALUES IN ASCENDING C ORDER. WHEN IJUMP .GT. 1 THE ROUTINE COMPUTES VALUES C C 2*COS(J*PI/L) , J=1,2,...,L AND J .NE. 0(MOD N/IJUMP+1) C C WHERE L = IJUMP*(N/IJUMP+1). C C C WHEN IJUMP = 1 IT COMPUTES C C 2*COS((J-FNUM)*PI/(N+FDEN)) , J=1, 2, ... ,N C C WHERE C FNUM = 0.5, FDEN = 0.0, FOR REGULAR REDUCTION VALUES C FNUM = 0.0, FDEN = 1.0, FOR B-R AND C-R WHEN ISTAG = 1 C FNUM = 0.0, FDEN = 0.5, FOR B-R AND C-R WHEN ISTAG = 2 C FNUM = 0.5, FDEN = 0.5, FOR B-R AND C-R WHEN ISTAG = 2 C IN POISN2 ONLY. C C PI = PIMACH(DUM) IF (N .EQ. 0) GO TO 105 IF (IJUMP .EQ. 1) GO TO 103 K3 = N/IJUMP+1 K4 = K3-1 PIBYN = PI/dfloat(N+IJUMP) DO 102 K=1,IJUMP K1 = (K-1)*K3 K5 = (K-1)*K4 DO 101 I=1,K4 X = K1+I K2 = K5+I A(K2) = -2.d0*dcos(X*PIBYN) 101 CONTINUE 102 CONTINUE GO TO 105 103 CONTINUE NP1 = N+1 Y = PI/(dfloat(N)+FDEN) DO 104 I=1,N X = dfloat(NP1-I)-FNUM A(I) = 2.d0*dcos(X*Y) 104 CONTINUE 105 CONTINUE RETURN END SUBROUTINE MERGE (TCOS,I1,M1,I2,M2,I3) implicit double precision (a-h,o-z) DIMENSION TCOS(*) C C C THIS SUBROUTINE MERGES TWO ASCENDING STRINGS OF NUMBERS IN THE C ARRAY TCOS. THE FIRST STRING IS OF LENGTH M1 AND STARTS AT C TCOS(I1+1). THE SECOND STRING IS OF LENGTH M2 AND STARTS AT C TCOS(I2+1). THE MERGED STRING GOES INTO TCOS(I3+1). C C J1 = 1 J2 = 1 J = I3 IF (M1 .EQ. 0) GO TO 107 IF (M2 .EQ. 0) GO TO 104 101 J = J+1 L = J1+I1 X = TCOS(L) L = J2+I2 Y = TCOS(L) IF (X-Y) 102,102,103 102 TCOS(J) = X J1 = J1+1 IF (J1 .GT. M1) GO TO 106 GO TO 101 103 TCOS(J) = Y J2 = J2+1 IF (J2 .LE. M2) GO TO 101 IF (J1 .GT. M1) GO TO 109 104 K = J-J1+1 DO 105 J=J1,M1 M = K+J L = J+I1 TCOS(M) = TCOS(L) 105 CONTINUE GO TO 109 106 CONTINUE IF (J2 .GT. M2) GO TO 109 107 K = J-J2+1 DO 108 J=J2,M2 M = K+J L = J+I2 TCOS(M) = TCOS(L) 108 CONTINUE 109 CONTINUE RETURN END SUBROUTINE TRI3 (M,A,B,C,K,Y1,Y2,Y3,TCOS,D,W1,W2,W3) implicit double precision (a-h,o-z) DIMENSION A(*) ,B(*) ,C(*) ,K(4) , 1 TCOS(*) ,Y1(*) ,Y2(*) ,Y3(*) , 2 D(*) ,W1(*) ,W2(*) ,W3(*) C C SUBROUTINE TO SOLVE THREE LINEAR SYSTEMS WHOSE COMMON COEFFICIENT C MATRIX IS A RATIONAL FUNCTION IN THE MATRIX GIVEN BY C C TRIDIAGONAL (...,A(I),B(I),C(I),...) C MM1 = M-1 K1 = K(1) K2 = K(2) K3 = K(3) K4 = K(4) F1 = K1+1 F2 = K2+1 F3 = K3+1 F4 = K4+1 K2K3K4 = K2+K3+K4 IF (K2K3K4 .EQ. 0) GO TO 101 L1 = F1/F2 L2 = F1/F3 L3 = F1/F4 LINT1 = 1 LINT2 = 1 LINT3 = 1 KINT1 = K1 KINT2 = KINT1+K2 KINT3 = KINT2+K3 101 CONTINUE DO 115 N=1,K1 X = TCOS(N) IF (K2K3K4 .EQ. 0) GO TO 107 IF (N .NE. L1) GO TO 103 DO 102 I=1,M W1(I) = Y1(I) 102 CONTINUE 103 IF (N .NE. L2) GO TO 105 DO 104 I=1,M W2(I) = Y2(I) 104 CONTINUE 105 IF (N .NE. L3) GO TO 107 DO 106 I=1,M W3(I) = Y3(I) 106 CONTINUE 107 CONTINUE Z = 1.d0/(B(1)-X) D(1) = C(1)*Z Y1(1) = Y1(1)*Z Y2(1) = Y2(1)*Z Y3(1) = Y3(1)*Z DO 108 I=2,M Z = 1.d0/(B(I)-X-A(I)*D(I-1)) D(I) = C(I)*Z Y1(I) = (Y1(I)-A(I)*Y1(I-1))*Z Y2(I) = (Y2(I)-A(I)*Y2(I-1))*Z Y3(I) = (Y3(I)-A(I)*Y3(I-1))*Z 108 CONTINUE DO 109 IP=1,MM1 I = M-IP Y1(I) = Y1(I)-D(I)*Y1(I+1) Y2(I) = Y2(I)-D(I)*Y2(I+1) Y3(I) = Y3(I)-D(I)*Y3(I+1) 109 CONTINUE IF (K2K3K4 .EQ. 0) GO TO 115 IF (N .NE. L1) GO TO 111 I = LINT1+KINT1 XX = X-TCOS(I) DO 110 I=1,M Y1(I) = XX*Y1(I)+W1(I) 110 CONTINUE LINT1 = LINT1+1 L1 = (dfloat(LINT1)*F1)/F2 111 IF (N .NE. L2) GO TO 113 I = LINT2+KINT2 XX = X-TCOS(I) DO 112 I=1,M Y2(I) = XX*Y2(I)+W2(I) 112 CONTINUE LINT2 = LINT2+1 L2 = (dfloat(LINT2)*F1)/F3 113 IF (N .NE. L3) GO TO 115 I = LINT3+KINT3 XX = X-TCOS(I) DO 114 I=1,M Y3(I) = XX*Y3(I)+W3(I) 114 CONTINUE LINT3 = LINT3+1 L3 = (dfloat(LINT3)*F1)/F4 115 CONTINUE RETURN END double precision FUNCTION PIMACH (DUM) C C THIS SUBPROGRAM SUPPLIES THE VALUE OF THE CONSTANT PI CORRECT TO C MACHINE PRECISION WHERE C C PI=3.1415926535897932384626433832795028841971693993751058209749446 C implicit double precision (a-h,o-z) PIMACH = 3.141592653589793238d0 RETURN END