SLATEC Routines --- HW3CRT ---
*DECK HW3CRT
SUBROUTINE HW3CRT (XS, XF, L, LBDCND, BDXS, BDXF, YS, YF, M,
+ MBDCND, BDYS, BDYF, ZS, ZF, N, NBDCND, BDZS, BDZF, ELMBDA,
+ LDIMF, MDIMF, F, PERTRB, IERROR, W)
C***BEGIN PROLOGUE HW3CRT
C***PURPOSE Solve the standard seven-point finite difference
C approximation to the Helmholtz equation in Cartesian
C coordinates.
C***LIBRARY SLATEC (FISHPACK)
C***CATEGORY I2B1A1A
C***TYPE SINGLE PRECISION (HW3CRT-S)
C***KEYWORDS CARTESIAN, ELLIPTIC, FISHPACK, HELMHOLTZ, PDE
C***AUTHOR Adams, J., (NCAR)
C Swarztrauber, P. N., (NCAR)
C Sweet, R., (NCAR)
C***DESCRIPTION
C
C Subroutine HW3CRT solves the standard seven-point finite
C difference approximation to the Helmholtz equation in Cartesian
C coordinates:
C
C (d/dX)(dU/dX) + (d/dY)(dU/dY) + (d/dZ)(dU/dZ)
C
C + LAMBDA*U = F(X,Y,Z) .
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
C * * * * * * * * Parameter Description * * * * * * * * * *
C
C
C * * * * * * On Input * * * * * *
C
C XS,XF
C The range of X, i.e. XS .LE. X .LE. XF .
C XS must be less than XF.
C
C L
C The number of panels into which the interval (XS,XF) is
C subdivided. Hence, there will be L+1 grid points in the
C X-direction given by X(I) = XS+(I-1)DX for I=1,2,...,L+1,
C where DX = (XF-XS)/L is the panel width. L must be at
C least 5 .
C
C LBDCND
C Indicates the type of boundary conditions at X = XS and X = XF.
C
C = 0 If the solution is periodic in X, i.e.
C U(L+I,J,K) = U(I,J,K).
C = 1 If the solution is specified at X = XS and X = XF.
C = 2 If the solution is specified at X = XS and the derivative
C of the solution with respect to X is specified at X = XF.
C = 3 If the derivative of the solution with respect to X is
C specified at X = XS and X = XF.
C = 4 If the derivative of the solution with respect to X is
C specified at X = XS and the solution is specified at X=XF.
C
C BDXS
C A two-dimensional array that specifies the values of the
C derivative of the solution with respect to X at X = XS.
C when LBDCND = 3 or 4,
C
C BDXS(J,K) = (d/dX)U(XS,Y(J),Z(K)), J=1,2,...,M+1,
C K=1,2,...,N+1.
C
C When LBDCND has any other value, BDXS is a dummy variable.
C BDXS must be dimensioned at least (M+1)*(N+1).
C
C BDXF
C A two-dimensional array that specifies the values of the
C derivative of the solution with respect to X at X = XF.
C When LBDCND = 2 or 3,
C
C BDXF(J,K) = (d/dX)U(XF,Y(J),Z(K)), J=1,2,...,M+1,
C K=1,2,...,N+1.
C
C When LBDCND has any other value, BDXF is a dummy variable.
C BDXF must be dimensioned at least (M+1)*(N+1).
C
C YS,YF
C The range of Y, i.e. YS .LE. Y .LE. YF.
C YS must be less than YF.
C
C M
C The number of panels into which the interval (YS,YF) is
C subdivided. Hence, there will be M+1 grid points in the
C Y-direction given by Y(J) = YS+(J-1)DY for J=1,2,...,M+1,
C where DY = (YF-YS)/M is the panel width. M must be at
C least 5 .
C
C MBDCND
C Indicates the type of boundary conditions at Y = YS and Y = YF.
C
C = 0 If the solution is periodic in Y, i.e.
C U(I,M+J,K) = U(I,J,K).
C = 1 If the solution is specified at Y = YS and Y = YF.
C = 2 If the solution is specified at Y = YS and the derivative
C of the solution with respect to Y is specified at Y = YF.
C = 3 If the derivative of the solution with respect to Y is
C specified at Y = YS and Y = YF.
C = 4 If the derivative of the solution with respect to Y is
C specified at Y = YS and the solution is specified at Y=YF.
C
C BDYS
C A two-dimensional array that specifies the values of the
C derivative of the solution with respect to Y at Y = YS.
C When MBDCND = 3 or 4,
C
C BDYS(I,K) = (d/dY)U(X(I),YS,Z(K)), I=1,2,...,L+1,
C K=1,2,...,N+1.
C
C When MBDCND has any other value, BDYS is a dummy variable.
C BDYS must be dimensioned at least (L+1)*(N+1).
C
C BDYF
C A two-dimensional array that specifies the values of the
C derivative of the solution with respect to Y at Y = YF.
C When MBDCND = 2 or 3,
C
C BDYF(I,K) = (d/dY)U(X(I),YF,Z(K)), I=1,2,...,L+1,
C K=1,2,...,N+1.
C
C When MBDCND has any other value, BDYF is a dummy variable.
C BDYF must be dimensioned at least (L+1)*(N+1).
C
C ZS,ZF
C The range of Z, i.e. ZS .LE. Z .LE. ZF.
C ZS must be less than ZF.
C
C N
C The number of panels into which the interval (ZS,ZF) is
C subdivided. Hence, there will be N+1 grid points in the
C Z-direction given by Z(K) = ZS+(K-1)DZ for K=1,2,...,N+1,
C where DZ = (ZF-ZS)/N is the panel width. N must be at least 5.
C
C NBDCND
C Indicates the type of boundary conditions at Z = ZS and Z = ZF.
C
C = 0 If the solution is periodic in Z, i.e.
C U(I,J,N+K) = U(I,J,K).
C = 1 If the solution is specified at Z = ZS and Z = ZF.
C = 2 If the solution is specified at Z = ZS and the derivative
C of the solution with respect to Z is specified at Z = ZF.
C = 3 If the derivative of the solution with respect to Z is
C specified at Z = ZS and Z = ZF.
C = 4 If the derivative of the solution with respect to Z is
C specified at Z = ZS and the solution is specified at Z=ZF.
C
C BDZS
C A two-dimensional array that specifies the values of the
C derivative of the solution with respect to Z at Z = ZS.
C When NBDCND = 3 or 4,
C
C BDZS(I,J) = (d/dZ)U(X(I),Y(J),ZS), I=1,2,...,L+1,
C J=1,2,...,M+1.
C
C When NBDCND has any other value, BDZS is a dummy variable.
C BDZS must be dimensioned at least (L+1)*(M+1).
C
C BDZF
C A two-dimensional array that specifies the values of the
C derivative of the solution with respect to Z at Z = ZF.
C When NBDCND = 2 or 3,
C
C BDZF(I,J) = (d/dZ)U(X(I),Y(J),ZF), I=1,2,...,L+1,
C J=1,2,...,M+1.
C
C When NBDCND has any other value, BDZF is a dummy variable.
C BDZF must be dimensioned at least (L+1)*(M+1).
C
C ELMBDA
C The constant LAMBDA in the Helmholtz equation. If
C LAMBDA .GT. 0, a solution may not exist. However, HW3CRT will
C attempt to find a solution.
C
C F
C A three-dimensional array that specifies the values of the
C right side of the Helmholtz equation and boundary values (if
C any). For I=2,3,...,L, J=2,3,...,M, and K=2,3,...,N
C
C F(I,J,K) = F(X(I),Y(J),Z(K)).
C
C On the boundaries F is defined by
C
C LBDCND F(1,J,K) F(L+1,J,K)
C ------ --------------- ---------------
C
C 0 F(XS,Y(J),Z(K)) F(XS,Y(J),Z(K))
C 1 U(XS,Y(J),Z(K)) U(XF,Y(J),Z(K))
C 2 U(XS,Y(J),Z(K)) F(XF,Y(J),Z(K)) J=1,2,...,M+1
C 3 F(XS,Y(J),Z(K)) F(XF,Y(J),Z(K)) K=1,2,...,N+1
C 4 F(XS,Y(J),Z(K)) U(XF,Y(J),Z(K))
C
C MBDCND F(I,1,K) F(I,M+1,K)
C ------ --------------- ---------------
C
C 0 F(X(I),YS,Z(K)) F(X(I),YS,Z(K))
C 1 U(X(I),YS,Z(K)) U(X(I),YF,Z(K))
C 2 U(X(I),YS,Z(K)) F(X(I),YF,Z(K)) I=1,2,...,L+1
C 3 F(X(I),YS,Z(K)) F(X(I),YF,Z(K)) K=1,2,...,N+1
C 4 F(X(I),YS,Z(K)) U(X(I),YF,Z(K))
C
C NBDCND F(I,J,1) F(I,J,N+1)
C ------ --------------- ---------------
C
C 0 F(X(I),Y(J),ZS) F(X(I),Y(J),ZS)
C 1 U(X(I),Y(J),ZS) U(X(I),Y(J),ZF)
C 2 U(X(I),Y(J),ZS) F(X(I),Y(J),ZF) I=1,2,...,L+1
C 3 F(X(I),Y(J),ZS) F(X(I),Y(J),ZF) J=1,2,...,M+1
C 4 F(X(I),Y(J),ZS) U(X(I),Y(J),ZF)
C
C F must be dimensioned at least (L+1)*(M+1)*(N+1).
C
C NOTE:
C
C If the table calls for both the solution U and the right side F
C on a boundary, then the solution must be specified.
C
C LDIMF
C The row (or first) dimension of the arrays F,BDYS,BDYF,BDZS,
C and BDZF as it appears in the program calling HW3CRT. this
C parameter is used to specify the variable dimension of these
C arrays. LDIMF must be at least L+1.
C
C MDIMF
C The column (or second) dimension of the array F and the row (or
C first) dimension of the arrays BDXS and BDXF as it appears in
C the program calling HW3CRT. This parameter is used to specify
C the variable dimension of these arrays.
C MDIMF 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. The length of W must be at least 30 + L + M + 5*N
C + MAX(L,M,N) + 7*(INT((L+1)/2) + INT((M+1)/2))
C
C
C * * * * * * On Output * * * * * *
C
C F
C Contains the solution U(I,J,K) of the finite difference
C approximation for the grid point (X(I),Y(J),Z(K)) for
C I=1,2,...,L+1, J=1,2,...,M+1, and K=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
C may not exist. PERTRB is a constant, calculated and subtracted
C from F, which ensures that a solution exists. PWSCRT then
C computes this solution, which is a least squares solution to
C the original approximation. This solution is not unique and is
C unnormalized. The value of PERTRB should be small compared to
C the right side F. Otherwise, a solution is obtained to an
C essentially different problem. This comparison should always
C be made to insure that a meaningful solution has been obtained.
C
C IERROR
C An error flag that indicates invalid input parameters. Except
C for numbers 0 and 12, a solution is not attempted.
C
C = 0 No error
C = 1 XS .GE. XF
C = 2 L .LT. 5
C = 3 LBDCND .LT. 0 .OR. LBDCND .GT. 4
C = 4 YS .GE. YF
C = 5 M .LT. 5
C = 6 MBDCND .LT. 0 .OR. MBDCND .GT. 4
C = 7 ZS .GE. ZF
C = 8 N .LT. 5
C = 9 NBDCND .LT. 0 .OR. NBDCND .GT. 4
C = 10 LDIMF .LT. L+1
C = 11 MDIMF .LT. M+1
C = 12 LAMBDA .GT. 0
C
C Since this is the only means of indicating a possibly incorrect
C call to HW3CRT, the user should test IERROR after the call.
C
C *Long Description:
C
C * * * * * * * Program Specifications * * * * * * * * * * * *
C
C Dimension of BDXS(MDIMF,N+1),BDXF(MDIMF,N+1),BDYS(LDIMF,N+1),
C Arguments BDYF(LDIMF,N+1),BDZS(LDIMF,M+1),BDZF(LDIMF,M+1),
C F(LDIMF,MDIMF,N+1),W(see argument list)
C
C Latest December 1, 1978
C Revision
C
C Subprograms HW3CRT,POIS3D,POS3D1,TRIDQ,RFFTI,RFFTF,RFFTF1,
C Required RFFTB,RFFTB1,COSTI,COST,SINTI,SINT,COSQI,COSQF,
C COSQF1,COSQB,COSQB1,SINQI,SINQF,SINQB,CFFTI,
C CFFTI1,CFFTB,CFFTB1,PASSB2,PASSB3,PASSB4,PASSB,
C CFFTF,CFFTF1,PASSF1,PASSF2,PASSF3,PASSF4,PASSF,
C 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 Written by Roland Sweet at NCAR in July 1977
C
C Algorithm This subroutine defines the finite difference
C equations, incorporates boundary data, and
C adjusts the right side of singular systems and
C then calls POIS3D to solve the system.
C
C Space 7862(decimal) = 17300(octal) locations on the
C Required NCAR Control Data 7600
C
C Timing and The execution time T on the NCAR Control Data
C Accuracy 7600 for subroutine HW3CRT is roughly proportional
C to L*M*N*(log2(L)+log2(M)+5), but also depends on
C input parameters LBDCND and MBDCND. Some typical
C values are listed in the table below.
C The solution process employed results in a loss
C of no more than three significant digits for L,M
C and N as large as 32. More detailed information
C about accuracy can be found in the documentation
C for subroutine POIS3D which is the routine that
C actually solves the finite difference equations.
C
C
C L(=M=N) LBDCND(=MBDCND=NBDCND) T(MSECS)
C ------- ---------------------- --------
C
C 16 0 300
C 16 1 302
C 16 3 348
C 32 0 1925
C 32 1 1929
C 32 3 2109
C
C Portability American National Standards Institute FORTRAN.
C The machine dependent constant PI is defined in
C function PIMACH.
C
C Required COS,SIN,ATAN
C Resident
C Routines
C
C Reference NONE
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C***REFERENCES (NONE)
C***ROUTINES CALLED POIS3D
C***REVISION HISTORY (YYMMDD)
C 801001 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890531 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C***END PROLOGUE HW3CRT