2
SUBROUTINE HSTPLR (A, B, M, MBDCND, BDA, BDB, C, D, N, NBDCND,
3
+ BDC, BDD, ELMBDA, F, IDIMF, PERTRB, IERROR, W)
4
C***BEGIN PROLOGUE HSTPLR
5
C***PURPOSE Solve the standard five-point finite difference
6
C approximation on a staggered grid to the Helmholtz equation
7
C in polar coordinates.
8
C***LIBRARY SLATEC (FISHPACK)
10
C***TYPE SINGLE PRECISION (HSTPLR-S)
11
C***KEYWORDS ELLIPTIC, FISHPACK, HELMHOLTZ, PDE, POLAR
12
C***AUTHOR Adams, J., (NCAR)
13
C Swarztrauber, P. N., (NCAR)
17
C HSTPLR solves the standard five-point finite difference
18
C approximation on a staggered grid to the Helmholtz equation in
21
C (1/R)(d/DR)(R(dU/DR)) + (1/R**2)(d/dTHETA)(dU/dTHETA)
23
C + LAMBDA*U = F(R,THETA)
25
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
27
C * * * * * * * * Parameter Description * * * * * * * * * *
29
C * * * * * * On Input * * * * * *
32
C The range of R, i.e. A .LE. R .LE. B. A must be less than B and
33
C A must be non-negative.
36
C The number of grid points in the interval (A,B). The grid points
37
C in the R-direction are given by R(I) = A + (I-0.5)DR for
38
C I=1,2,...,M where DR =(B-A)/M. M must be greater than 2.
41
C Indicates the type of boundary conditions at R = A and R = B.
43
C = 1 If the solution is specified at R = A and R = B.
45
C = 2 If the solution is specified at R = A and the derivative
46
C of the solution with respect to R is specified at R = B.
49
C = 3 If the derivative of the solution with respect to R is
50
C specified at R = A (see note 2 below) and R = B.
52
C = 4 If the derivative of the solution with respect to R is
53
C specified at R = A (see note 2 below) and the solution is
56
C = 5 If the solution is unspecified at R = A = 0 and the solution
57
C is specified at R = B.
59
C = 6 If the solution is unspecified at R = A = 0 and the
60
C derivative of the solution with respect to R is specified at
63
C NOTE 1: If A = 0, MBDCND = 2, and NBDCND = 0 or 3, the system of
64
C equations to be solved is singular. The unique solution
65
C is determined by extrapolation to the specification of
66
C U(0,THETA(1)). But in this case the right side of the
67
C system will be perturbed by the constant PERTRB.
69
C NOTE 2: If A = 0, do not use MBDCND = 3 or 4, but instead use
70
C MBDCND = 1,2,5, or 6.
73
C A one-dimensional array of length N that specifies the boundary
74
C values (if any) of the solution at R = A. When MBDCND = 1 or 2,
76
C BDA(J) = U(A,THETA(J)) , J=1,2,...,N.
78
C When MBDCND = 3 or 4,
80
C BDA(J) = (d/dR)U(A,THETA(J)) , J=1,2,...,N.
82
C When MBDCND = 5 or 6, BDA is a dummy variable.
85
C A one-dimensional array of length N that specifies the boundary
86
C values of the solution at R = B. When MBDCND = 1,4, or 5,
88
C BDB(J) = U(B,THETA(J)) , J=1,2,...,N.
90
C When MBDCND = 2,3, or 6,
92
C BDB(J) = (d/dR)U(B,THETA(J)) , J=1,2,...,N.
95
C The range of THETA, i.e. C .LE. THETA .LE. D. C must be less
99
C The number of unknowns in the interval (C,D). The unknowns in
100
C the THETA-direction are given by THETA(J) = C + (J-0.5)DT,
101
C J=1,2,...,N, where DT = (D-C)/N. N must be greater than 2.
104
C Indicates the type of boundary conditions at THETA = C
107
C = 0 If the solution is periodic in THETA, i.e.
110
C = 1 If the solution is specified at THETA = C and THETA = D
113
C = 2 If the solution is specified at THETA = C and the derivative
114
C of the solution with respect to THETA is specified at
115
C THETA = D (see note below).
117
C = 3 If the derivative of the solution with respect to THETA is
118
C specified at THETA = C and THETA = D.
120
C = 4 If the derivative of the solution with respect to THETA is
121
C specified at THETA = C and the solution is specified at
122
C THETA = d (see note below).
124
C NOTE: When NBDCND = 1, 2, or 4, do not use MBDCND = 5 or 6 (the
125
C former indicates that the solution is specified at R = 0; the
126
C latter indicates the solution is unspecified at R = 0). Use
127
C instead MBDCND = 1 or 2.
130
C A one dimensional array of length M that specifies the boundary
131
C values of the solution at THETA = C. When NBDCND = 1 or 2,
133
C BDC(I) = U(R(I),C) , I=1,2,...,M.
135
C When NBDCND = 3 or 4,
137
C BDC(I) = (d/dTHETA)U(R(I),C), I=1,2,...,M.
139
C When NBDCND = 0, BDC is a dummy variable.
142
C A one-dimensional array of length M that specifies the boundary
143
C values of the solution at THETA = D. When NBDCND = 1 or 4,
145
C BDD(I) = U(R(I),D) , I=1,2,...,M.
147
C When NBDCND = 2 or 3,
149
C BDD(I) = (d/dTHETA)U(R(I),D) , I=1,2,...,M.
151
C When NBDCND = 0, BDD is a dummy variable.
154
C The constant LAMBDA in the Helmholtz equation. If LAMBDA is
155
C greater than 0, a solution may not exist. However, HSTPLR will
156
C attempt to find a solution.
159
C A two-dimensional array that specifies the values of the right
160
C side of the Helmholtz equation. For I=1,2,...,M and J=1,2,...,N
162
C F(I,J) = F(R(I),THETA(J)) .
164
C F must be dimensioned at least M X N.
167
C The row (or first) dimension of the array F as it appears in the
168
C program calling HSTPLR. This parameter is used to specify the
169
C variable dimension of F. IDIMF must be at least M.
172
C A one-dimensional array that must be provided by the user for
173
C work space. W may require up to 13M + 4N + M*INT(log2(N))
174
C locations. The actual number of locations used is computed by
175
C HSTPLR and is returned in the location W(1).
178
C * * * * * * On Output * * * * * *
181
C Contains the solution U(I,J) of the finite difference
182
C approximation for the grid point (R(I),THETA(J)) for
183
C I=1,2,...,M, J=1,2,...,N.
186
C If a combination of periodic, derivative, or unspecified
187
C boundary conditions is specified for a Poisson equation
188
C (LAMBDA = 0), a solution may not exist. PERTRB is a con-
189
C stant, calculated and subtracted from F, which ensures
190
C that a solution exists. HSTPLR then computes this
191
C solution, which is a least squares solution to the
192
C original approximation. This solution plus any constant is also
193
C a solution; hence, the solution is not unique. The value of
194
C PERTRB should be small compared to the right side F.
195
C Otherwise, a solution is obtained to an essentially different
196
C problem. This comparison should always be made to insure that
197
C a meaningful solution has been obtained.
200
C An error flag that indicates invalid input parameters.
201
C Except for numbers 0 and 11, a solution is not attempted.
209
C = 3 MBDCND .LT. 1 or MBDCND .GT. 6
215
C = 6 NBDCND .LT. 0 or NBDCND .GT. 4
217
C = 7 A = 0 and MBDCND = 3 or 4
219
C = 8 A .GT. 0 and MBDCND .GE. 5
221
C = 9 MBDCND .GE. 5 and NBDCND .NE. 0 or 3
229
C Since this is the only means of indicating a possibly
230
C incorrect call to HSTPLR, the user should test IERROR after
234
C W(1) contains the required length of W.
238
C * * * * * * * Program Specifications * * * * * * * * * * * *
240
C Dimension of BDA(N),BDB(N),BDC(M),BDD(M),F(IDIMF,N),
241
C Arguments W(see ARGUMENT LIST)
243
C Latest June 1, 1977
246
C Subprograms HSTPLR,POISTG,POSTG2,GENBUN,POISD2,POISN2,POISP2,
247
C Required COSGEN,MERGE,TRIX,TRI3,PIMACH
259
C Specialist Roland Sweet
263
C History Written by Roland Sweet at NCAR in February, 1977
265
C Algorithm This subroutine defines the finite-difference
266
C equations, incorporates boundary data, adjusts the
267
C right side when the system is singular and calls
268
C either POISTG or GENBUN which solves the linear
269
C system of equations.
271
C Space 8265(decimal) = 20111(octal) LOCATIONS ON THE
272
C Required NCAR Control Data 7600
274
C Timing and The execution time T on the NCAR Control Data
275
C Accuracy 7600 for subroutine HSTPLR is roughly proportional
276
C to M*N*log2(N). Some typical values are listed in
278
C The solution process employed results in a loss
279
C of no more than four significant digits for N and M
280
C as large as 64. More detailed information about
281
C accuracy can be found in the documentation for
282
C subroutine POISTG which is the routine that
283
C actually solves the finite difference equations.
286
C M(=N) MBDCND NBDCND T(MSECS)
287
C ----- ------ ------ --------
292
C Portability American National Standards Institute Fortran.
293
C The machine dependent constant PI is defined in
300
C Reference Schumann, U. and R. Sweet,'A Direct Method For
301
C The Solution Of Poisson's Equation With Neumann
302
C Boundary Conditions On A Staggered Grid of
303
C Arbitrary Size,' J. Comp. Phys. 20(1976),
306
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
308
C***REFERENCES U. Schumann and R. Sweet, A direct method for the
309
C solution of Poisson's equation with Neumann boundary
310
C conditions on a staggered grid of arbitrary size,
311
C Journal of Computational Physics 20, (1976),
313
C***ROUTINES CALLED GENBUN, POISTG
314
C***REVISION HISTORY (YYMMDD)
315
C 801001 DATE WRITTEN
316
C 890531 Changed all specific intrinsics to generic. (WRB)
317
C 890531 REVISION DATE from Version 3.2
318
C 891214 Prologue converted to Version 4.0 format. (BAB)
319
C 920501 Reformatted the REFERENCES section. (WRB)
320
C***END PROLOGUE HSTPLR
324
DIMENSION BDA(*) ,BDB(*) ,BDC(*) ,BDD(*) ,
326
C***FIRST EXECUTABLE STATEMENT HSTPLR
328
IF (A .LT. 0.) IERROR = 1
329
IF (A .GE. B) IERROR = 2
330
IF (MBDCND.LE.0 .OR. MBDCND.GE.7) IERROR = 3
331
IF (C .GE. D) IERROR = 4
332
IF (N .LE. 2) IERROR = 5
333
IF (NBDCND.LT.0 .OR. NBDCND.GE.5) IERROR = 6
334
IF (A.EQ.0. .AND. (MBDCND.EQ.3 .OR. MBDCND.EQ.4)) IERROR = 7
335
IF (A.GT.0. .AND. MBDCND.GE.5) IERROR = 8
336
IF (MBDCND.GE.5 .AND. NBDCND.NE.0 .AND. NBDCND.NE.3) IERROR = 9
337
IF (IDIMF .LT. M) IERROR = 10
338
IF (M .LE. 2) IERROR = 12
339
IF (IERROR .NE. 0) RETURN
347
IF (A.EQ.0. .AND. MBDCND.EQ.2) MB = 6
349
C DEFINE A,B,C COEFFICIENTS IN W-ARRAY.
356
W(J) = A+(I-0.5)*DELTAR
357
W(I) = (A+(I-1)*DELTAR)/DLRSQ
359
W(K) = (A+I*DELTAR)/DLRSQ
361
W(K) = (ELMBDA-2./DLRSQ)*W(J)
371
C ENTER BOUNDARY DATA FOR R-BOUNDARIES.
373
GO TO (104,104,106,106,108,108),MB
375
W(IWB+1) = W(IWB+1)-W(1)
377
F(1,J) = F(1,J)-A1*BDA(J)
381
W(IWB+1) = W(IWB+1)+W(1)
383
F(1,J) = F(1,J)+A1*BDA(J)
385
108 GO TO (109,111,111,109,109,111),MB
387
W(IWC) = W(IWC)-W(IWR)
389
F(M,J) = F(M,J)-A1*BDB(J)
392
111 A1 = DELTAR*W(IWR)
393
W(IWC) = W(IWC)+W(IWR)
395
F(M,J) = F(M,J)-A1*BDB(J)
398
C ENTER BOUNDARY DATA FOR THETA-BOUNDARIES.
401
GO TO (123,114,114,116,116),NP
404
F(I,1) = F(I,1)-A1*BDC(I)/W(J)
410
F(I,1) = F(I,1)+A1*BDC(I)/W(J)
413
GO TO (123,119,121,121,119),NP
416
F(I,N) = F(I,N)-A1*BDD(I)/W(J)
422
F(I,N) = F(I,N)-A1*BDD(I)/W(J)
426
C ADJUST RIGHT SIDE OF SINGULAR PROBLEMS TO INSURE EXISTENCE OF A
430
IF (ELMBDA) 133,125,124
433
125 GO TO (133,133,126,133,133,126),MB
434
126 GO TO (127,133,133,127,133),NP
439
PERTRB = PERTRB+F(I,J)
442
PERTRB = PERTRB/(M*N*0.5*(A+B))
457
C MULTIPLY I-TH EQUATION THROUGH BY R(I)*DELTHT**2
475
C CALL POISTG OR GENBUN TO SOLVE THE SYSTEM OF EQUATIONS.
477
IF (LP .EQ. 0) GO TO 136
478
CALL POISTG (LP,N,1,M,W,W(IWB+1),W(IWC+1),IDIMF,F,IERR1,W(IWR+1))
480
136 CALL GENBUN (LP,N,1,M,W,W(IWB+1),W(IWC+1),IDIMF,F,IERR1,W(IWR+1))
483
IF (A.NE.0. .OR. MBDCND.NE.2 .OR. ISW.NE.2) GO TO 141
488
A1 = (A1-DLRSQ*A2/16.)/N
489
IF (NBDCND .EQ. 3) A1 = A1+(BDD(1)-BDC(1))/(D-C)
added
removed
deps/openlibm/slatec/hstplr.f
