2
SUBROUTINE HSTCSP (INTL, A, B, M, MBDCND, BDA, BDB, C, D, N,
3
+ NBDCND, BDC, BDD, ELMBDA, F, IDIMF, PERTRB, IERROR, W)
4
C***BEGIN PROLOGUE HSTCSP
5
C***PURPOSE Solve the standard five-point finite difference
6
C approximation on a staggered grid to the modified Helmholtz
7
C equation in spherical coordinates assuming axisymmetry
8
C (no dependence on longitude).
9
C***LIBRARY SLATEC (FISHPACK)
11
C***TYPE SINGLE PRECISION (HSTCSP-S)
12
C***KEYWORDS ELLIPTIC, FISHPACK, HELMHOLTZ, PDE, SPHERICAL
13
C***AUTHOR Adams, J., (NCAR)
14
C Swarztrauber, P. N., (NCAR)
18
C HSTCSP solves the standard five-point finite difference
19
C approximation on a staggered grid to the modified Helmholtz
20
C equation spherical coordinates assuming axisymmetry (no dependence
23
C (1/R**2)(d/dR)(R**2(dU/dR)) +
25
C 1/(R**2*SIN(THETA))(d/dTHETA)(SIN(THETA)(dU/dTHETA)) +
27
C (LAMBDA/(R*SIN(THETA))**2)U = F(THETA,R)
29
C where THETA is colatitude and R is the radial coordinate.
30
C This two-dimensional modified Helmholtz equation results from
31
C the Fourier transform of the three-dimensional Poisson equation.
33
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
36
C * * * * * * * * Parameter Description * * * * * * * * * *
39
C * * * * * * On Input * * * * * *
42
C = 0 On initial entry to HSTCSP or if any of the arguments
43
C C, D, N, or NBDCND are changed from a previous call.
45
C = 1 If C, D, N, and NBDCND are all unchanged from previous
48
C NOTE: A call with INTL = 0 takes approximately 1.5 times as much
49
C time as a call with INTL = 1. Once a call with INTL = 0
50
C has been made then subsequent solutions corresponding to
51
C different F, BDA, BDB, BDC, and BDD can be obtained
52
C faster with INTL = 1 since initialization is not repeated.
55
C The range of THETA (colatitude), i.e. A .LE. THETA .LE. B. A
56
C must be less than B and A must be non-negative. A and B are in
57
C radians. A = 0 corresponds to the north pole and B = PI
58
C corresponds to the south pole.
60
C * * * IMPORTANT * * *
62
C If B is equal to PI, then B must be computed using the statement
66
C This insures that B in the user's program is equal to PI in this
67
C program which permits several tests of the input parameters that
68
C otherwise would not be possible.
70
C * * * * * * * * * * * *
73
C The number of grid points in the interval (A,B). The grid points
74
C in the THETA-direction are given by THETA(I) = A + (I-0.5)DTHETA
75
C for I=1,2,...,M where DTHETA =(B-A)/M. M must be greater than 4.
78
C Indicates the type of boundary conditions at THETA = A and
81
C = 1 If the solution is specified at THETA = A and THETA = B.
82
C (See notes 1, 2 below)
84
C = 2 If the solution is specified at THETA = A and the derivative
85
C of the solution with respect to THETA is specified at
86
C THETA = B (See notes 1, 2 below).
88
C = 3 If the derivative of the solution with respect to THETA is
89
C specified at THETA = A (See notes 1, 2 below) and THETA = B.
91
C = 4 If the derivative of the solution with respect to THETA is
92
C specified at THETA = A (See notes 1, 2 below) and the
93
C solution is specified at THETA = B.
95
C = 5 If the solution is unspecified at THETA = A = 0 and the
96
C solution is specified at THETA = B. (See note 2 below)
98
C = 6 If the solution is unspecified at THETA = A = 0 and the
99
C derivative of the solution with respect to THETA is
100
C specified at THETA = B (See note 2 below).
102
C = 7 If the solution is specified at THETA = A and the
103
C solution is unspecified at THETA = B = PI.
105
C = 8 If the derivative of the solution with respect to
106
C THETA is specified at THETA = A (See note 1 below)
107
C and the solution is unspecified at THETA = B = PI.
109
C = 9 If the solution is unspecified at THETA = A = 0 and
112
C NOTES: 1. If A = 0, do not use MBDCND = 1,2,3,4,7 or 8,
113
C but instead use MBDCND = 5, 6, or 9.
115
C 2. if B = PI, do not use MBDCND = 1,2,3,4,5 or 6,
116
C but instead use MBDCND = 7, 8, or 9.
118
C When A = 0 and/or B = PI the only meaningful boundary
119
C condition is dU/dTHETA = 0. (See D. Greenspan, 'Numerical
120
C Analysis of Elliptic Boundary Value Problems,' Harper and
121
C Row, 1965, Chapter 5.)
124
C A one-dimensional array of length N that specifies the boundary
125
C values (if any) of the solution at THETA = A. When
126
C MBDCND = 1, 2, or 7,
128
C BDA(J) = U(A,R(J)) , J=1,2,...,N.
130
C When MBDCND = 3, 4, or 8,
132
C BDA(J) = (d/dTHETA)U(A,R(J)) , J=1,2,...,N.
134
C When MBDCND has any other value, BDA is a dummy variable.
137
C A one-dimensional array of length N that specifies the boundary
138
C values of the solution at THETA = B. When MBDCND = 1, 4, or 5,
140
C BDB(J) = U(B,R(J)) , J=1,2,...,N.
142
C When MBDCND = 2,3, or 6,
144
C BDB(J) = (d/dTHETA)U(B,R(J)) , J=1,2,...,N.
146
C When MBDCND has any other value, BDB is a dummy variable.
149
C The range of R , i.e. C .LE. R .LE. D.
150
C C must be less than D. C must be non-negative.
153
C The number of unknowns in the interval (C,D). The unknowns in
154
C the R-direction are given by R(J) = C + (J-0.5)DR,
155
C J=1,2,...,N, where DR = (D-C)/N. N must be greater than 4.
158
C Indicates the type of boundary conditions at R = C
161
C = 1 If the solution is specified at R = C and R = D.
163
C = 2 If the solution is specified at R = C and the derivative
164
C of the solution with respect to R is specified at
165
C R = D. (See note 1 below)
167
C = 3 If the derivative of the solution with respect to R is
168
C specified at R = C and R = D.
170
C = 4 If the derivative of the solution with respect to R is
171
C specified at R = C and the solution is specified at
174
C = 5 If the solution is unspecified at R = C = 0 (See note 2
175
C below) and the solution is specified at R = D.
177
C = 6 If the solution is unspecified at R = C = 0 (See note 2
178
C below) and the derivative of the solution with respect to R
179
C is specified at R = D.
181
C NOTE 1: If C = 0 and MBDCND = 3,6,8 or 9, the system of equations
182
C to be solved is singular. The unique solution is
183
C determined by extrapolation to the specification of
184
C U(THETA(1),C). But in these cases the right side of the
185
C system will be perturbed by the constant PERTRB.
187
C NOTE 2: NBDCND = 5 or 6 cannot be used with MBDCND = 1, 2, 4, 5,
188
C or 7 (the former indicates that the solution is
189
C unspecified at R = 0; the latter indicates that the
190
C solution is specified). Use instead NBDCND = 1 or 2.
193
C A one dimensional array of length M that specifies the boundary
194
C values of the solution at R = C. When NBDCND = 1 or 2,
196
C BDC(I) = U(THETA(I),C) , I=1,2,...,M.
198
C When NBDCND = 3 or 4,
200
C BDC(I) = (d/dR)U(THETA(I),C), I=1,2,...,M.
202
C When NBDCND has any other value, BDC is a dummy variable.
205
C A one-dimensional array of length M that specifies the boundary
206
C values of the solution at R = D. When NBDCND = 1 or 4,
208
C BDD(I) = U(THETA(I),D) , I=1,2,...,M.
210
C When NBDCND = 2 or 3,
212
C BDD(I) = (d/dR)U(THETA(I),D) , I=1,2,...,M.
214
C When NBDCND has any other value, BDD is a dummy variable.
217
C The constant LAMBDA in the modified Helmholtz equation. If
218
C LAMBDA is greater than 0, a solution may not exist. However,
219
C HSTCSP will attempt to find a solution.
222
C A two-dimensional array that specifies the values of the right
223
C side of the modified Helmholtz equation. For I=1,2,...,M and
226
C F(I,J) = F(THETA(I),R(J)) .
228
C F must be dimensioned at least M X N.
231
C The row (or first) dimension of the array F as it appears in the
232
C program calling HSTCSP. This parameter is used to specify the
233
C variable dimension of F. IDIMF must be at least M.
236
C A one-dimensional array that must be provided by the user for
237
C work space. With K = INT(log2(N))+1 and L = 2**(K+1), W may
238
C require up to (K-2)*L+K+MAX(2N,6M)+4(N+M)+5 locations. The
239
C actual number of locations used is computed by HSTCSP and is
240
C returned in the location W(1).
243
C * * * * * * On Output * * * * * *
246
C Contains the solution U(I,J) of the finite difference
247
C approximation for the grid point (THETA(I),R(J)) for
248
C I=1,2,...,M, J=1,2,...,N.
251
C If a combination of periodic, derivative, or unspecified
252
C boundary conditions is specified for a Poisson equation
253
C (LAMBDA = 0), a solution may not exist. PERTRB is a con-
254
C stant, calculated and subtracted from F, which ensures
255
C that a solution exists. HSTCSP then computes this
256
C solution, which is a least squares solution to the
257
C original approximation. This solution plus any constant is also
258
C a solution; hence, the solution is not unique. The value of
259
C PERTRB should be small compared to the right side F.
260
C Otherwise, a solution is obtained to an essentially different
261
C problem. This comparison should always be made to insure that
262
C a meaningful solution has been obtained.
265
C An error flag that indicates invalid input parameters.
266
C Except for numbers 0 and 10, a solution is not attempted.
270
C = 1 A .LT. 0 or B .GT. PI
274
C = 3 MBDCND .LT. 1 or MBDCND .GT. 9
280
C = 6 NBDCND .LT. 1 or NBDCND .GT. 6
284
C = 8 NBDCND = 5 or 6 and MBDCND = 1, 2, 4, 5, or 7
286
C = 9 C .GT. 0 and NBDCND .GE. 5
294
C = 13 A = 0 and MBDCND =1,2,3,4,7 or 8
296
C = 14 B = PI and MBDCND .LE. 6
298
C = 15 A .GT. 0 and MBDCND = 5, 6, or 9
300
C = 16 B .LT. PI and MBDCND .GE. 7
302
C = 17 LAMBDA .NE. 0 and NBDCND .GE. 5
304
C Since this is the only means of indicating a possibly
305
C incorrect call to HSTCSP, the user should test IERROR after
309
C W(1) contains the required length of W. Also W contains
310
C intermediate values that must not be destroyed if HSTCSP
311
C will be called again with INTL = 1.
315
C * * * * * * * Program Specifications * * * * * * * * * * * *
317
C Dimension of BDA(N),BDB(N),BDC(M),BDD(M),F(IDIMF,N),
318
C Arguments W(See argument list)
323
C Subprograms HSTCSP,HSTCS1,BLKTRI,BLKTR1,INDXA,INDXB,INDXC,
324
C Required PROD,PRODP,CPROD,CPRODP,PPADD,PSGF,BSRH,PPSGF,
325
C PPSPF,COMPB,TEVLS,R1MACH
337
C Specialist Roland Sweet
341
C History Written by Roland Sweet at NCAR in May, 1977
343
C Algorithm This subroutine defines the finite-difference
344
C equations, incorporates boundary data, adjusts the
345
C right side when the system is singular and calls
346
C BLKTRI which solves the linear system of equations.
348
C Space 5269(decimal) = 12225(octal) locations on the
349
C Required NCAR Control Data 7600
351
C Timing and The execution time T on the NCAR Control Data
352
C Accuracy 7600 for subroutine HSTCSP is roughly proportional
353
C to M*N*log2(N), but depends on the input parameter
354
C INTL. Some values are listed in the table below.
355
C The solution process employed results in a loss
356
C of no more than FOUR significant digits for N and M
357
C as large as 64. More detailed information about
358
C accuracy can be found in the documentation for
359
C subroutine BLKTRI which is the routine that
360
C actually solves the finite difference equations.
363
C M(=N) INTL MBDCND(=NBDCND) T(MSECS)
364
C ----- ---- --------------- --------
371
C Portability American National Standards Institute Fortran.
372
C The machine accuracy is set using function R1MACH.
374
C Required COS,SIN,ABS,SQRT
378
C Reference Swarztrauber, P.N., 'A Direct Method For The
379
C Discrete Solution Of Separable Elliptic Equations,'
380
C SIAM J. Numer. Anal. 11(1974), pp. 1136-1150.
382
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
384
C***REFERENCES P. N. Swarztrauber and R. Sweet, Efficient Fortran
385
C subprograms for the solution of elliptic equations,
386
C NCAR TN/IA-109, July 1975, 138 pp.
387
C P. N. Swarztrauber, A direct method for the discrete
388
C solution of separable elliptic equations, SIAM Journal
389
C on Numerical Analysis 11, (1974), pp. 1136-1150.
390
C***ROUTINES CALLED HSTCS1, PIMACH
391
C***REVISION HISTORY (YYMMDD)
392
C 801001 DATE WRITTEN
393
C 890531 Changed all specific intrinsics to generic. (WRB)
394
C 890531 REVISION DATE from Version 3.2
395
C 891214 Prologue converted to Version 4.0 format. (BAB)
396
C 920501 Reformatted the REFERENCES section. (WRB)
397
C***END PROLOGUE HSTCSP
400
DIMENSION F(IDIMF,*) ,BDA(*) ,BDB(*) ,BDC(*) ,
402
C***FIRST EXECUTABLE STATEMENT HSTCSP
405
C CHECK FOR INVALID INPUT PARAMETERS
408
IF (A.LT.0. .OR. B.GT.PI) IERROR = 1
409
IF (A .GE. B) IERROR = 2
410
IF (MBDCND.LT.1 .OR. MBDCND.GT.9) IERROR = 3
411
IF (C .LT. 0.) IERROR = 4
412
IF (C .GE. D) IERROR = 5
413
IF (NBDCND.LT.1 .OR. NBDCND.GT.6) IERROR = 6
414
IF (N .LT. 5) IERROR = 7
415
IF ((NBDCND.EQ.5 .OR. NBDCND.EQ.6) .AND. (MBDCND.EQ.1 .OR.
416
1 MBDCND.EQ.2 .OR. MBDCND.EQ.4 .OR. MBDCND.EQ.5 .OR.
419
IF (C.GT.0. .AND. NBDCND.GE.5) IERROR = 9
420
IF (IDIMF .LT. M) IERROR = 11
421
IF (M .LT. 5) IERROR = 12
422
IF (A.EQ.0. .AND. MBDCND.NE.5 .AND. MBDCND.NE.6 .AND. MBDCND.NE.9)
424
IF (B.EQ.PI .AND. MBDCND.LE.6) IERROR = 14
425
IF (A.GT.0. .AND. (MBDCND.EQ.5 .OR. MBDCND.EQ.6 .OR. MBDCND.EQ.9))
427
IF (B.LT.PI .AND. MBDCND.GE.7) IERROR = 16
428
IF (ELMBDA.NE.0. .AND. NBDCND.GE.5) IERROR = 17
429
IF (IERROR .NE. 0) GO TO 101
439
CALL HSTCS1 (INTL,A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC,BDD,
440
1 ELMBDA,F,IDIMF,PERTRB,IERR1,W,W(IWBM),W(IWCM),
441
2 W(IWAN),W(IWBN),W(IWCN),W(IWSNTH),W(IWRSQ),W(IWWRK))
442
W(1) = W(IWWRK)+IWWRK-1
added
removed
deps/openlibm/slatec/hstcsp.f
