relibc_openlibm/slatec/hwsssp.f

*DECK HWSSSP
      SUBROUTINE HWSSSP (TS, TF, M, MBDCND, BDTS, BDTF, PS, PF, N,
     +   NBDCND, BDPS, BDPF, ELMBDA, F, IDIMF, PERTRB, IERROR, W)
C***BEGIN PROLOGUE  HWSSSP
C***PURPOSE  Solve a finite difference approximation to the Helmholtz
C            equation in spherical coordinates and on the surface of the
C            unit sphere (radius of 1).
C***LIBRARY   SLATEC (FISHPACK)
C***CATEGORY  I2B1A1A
C***TYPE      SINGLE PRECISION (HWSSSP-S)
C***KEYWORDS  ELLIPTIC, FISHPACK, HELMHOLTZ, PDE, SPHERICAL
C***AUTHOR  Adams, J., (NCAR)
C           Swarztrauber, P. N., (NCAR)
C           Sweet, R., (NCAR)
C***DESCRIPTION
C
C     Subroutine HWSSSP solves a finite difference approximation to the
C     Helmholtz equation in spherical coordinates and on the surface of
C     the unit sphere (radius of 1):
C
C          (1/SIN(THETA))(d/dTHETA)(SIN(THETA)(dU/dTHETA))
C
C             + (1/SIN(THETA)**2)(d/dPHI)(dU/dPHI)
C
C             + LAMBDA*U = F(THETA,PHI)
C
C     Where THETA is colatitude and PHI is longitude.
C
C     * * * * * * * *    Parameter Description     * * * * * * * * * *
C
C             * * * * * *   On Input    * * * * * *
C
C     TS,TF
C       The range of THETA (colatitude), i.e., TS .LE. THETA .LE. TF.
C       TS must be less than TF.  TS and TF are in radians.  A TS of
C       zero corresponds to the north pole and a TF of PI corresponds to
C       the south pole.
C
C     * * * * * * * * * * * * * * IMPORTANT * * * * * * * * * * * * * *
C
C     If TF is equal to PI then it must be computed using the statement
C     TF = PIMACH(DUM). This insures that TF in the users program is
C     equal to PI in this program which permits several tests of the
C     input parameters that otherwise would not be possible.
C
C
C     M
C       The number of panels into which the interval (TS,TF) is
C       subdivided.  Hence, there will be M+1 grid points in the
C       THETA-direction given by THETA(I) = (I-1)DTHETA+TS for
C       I = 1,2,...,M+1, where DTHETA = (TF-TS)/M is the panel width.
C       M must be greater than 5.
C
C     MBDCND
C       Indicates the type of boundary condition at THETA = TS and
C       THETA = TF.
C
C       = 1  If the solution is specified at THETA = TS and THETA = TF.
C       = 2  If the solution is specified at THETA = TS and the
C            derivative of the solution with respect to THETA is
C            specified at THETA = TF (see note 2 below).
C       = 3  If the derivative of the solution with respect to THETA is
C            specified at THETA = TS and THETA = TF (see notes 1,2
C            below).
C       = 4  If the derivative of the solution with respect to THETA is
C            specified at THETA = TS (see note 1 below) and the
C            solution is specified at THETA = TF.
C       = 5  If the solution is unspecified at THETA = TS = 0 and the
C            solution is specified at THETA = TF.
C       = 6  If the solution is unspecified at THETA = TS = 0 and the
C            derivative of the solution with respect to THETA is
C            specified at THETA = TF (see note 2 below).
C       = 7  If the solution is specified at THETA = TS and the
C            solution is unspecified at THETA = TF = PI.
C       = 8  If the derivative of the solution with respect to THETA is
C            specified at THETA = TS (see note 1 below) and the
C            solution is unspecified at THETA = TF = PI.
C       = 9  If the solution is unspecified at THETA = TS = 0 and
C            THETA = TF = PI.
C
C       NOTES:  1.  If TS = 0, do not use MBDCND = 3,4, or 8, but
C                   instead use MBDCND = 5,6, or 9  .
C               2.  If TF = PI, do not use MBDCND = 2,3, or 6, but
C                   instead use MBDCND = 7,8, or 9  .
C
C     BDTS
C       A one-dimensional array of length N+1 that specifies the values
C       of the derivative of the solution with respect to THETA at
C       THETA = TS.  When MBDCND = 3,4, or 8,
C
C            BDTS(J) = (d/dTHETA)U(TS,PHI(J)), J = 1,2,...,N+1  .
C
C       When MBDCND has any other value, BDTS is a dummy variable.
C
C     BDTF
C       A one-dimensional array of length N+1 that specifies the values
C       of the derivative of the solution with respect to THETA at
C       THETA = TF.  When MBDCND = 2,3, or 6,
C
C            BDTF(J) = (d/dTHETA)U(TF,PHI(J)), J = 1,2,...,N+1  .
C
C       When MBDCND has any other value, BDTF is a dummy variable.
C
C     PS,PF
C       The range of PHI (longitude), i.e., PS .LE. PHI .LE. PF.  PS
C       must be less than PF.  PS and PF are in radians.  If PS = 0 and
C       PF = 2*PI, periodic boundary conditions are usually prescribed.
C
C     * * * * * * * * * * * * * * IMPORTANT * * * * * * * * * * * * * *
C
C     If PF is equal to 2*PI then it must be computed using the
C     statement PF = 2.*PIMACH(DUM). This insures that PF in the users
C     program is equal to 2*PI in this program which permits tests of
C     the input parameters that otherwise would not be possible.
C
C
C     N
C       The number of panels into which the interval (PS,PF) is
C       subdivided.  Hence, there will be N+1 grid points in the
C       PHI-direction given by PHI(J) = (J-1)DPHI+PS  for
C       J = 1,2,...,N+1, where DPHI = (PF-PS)/N is the panel width.
C       N must be greater than 4.
C
C     NBDCND
C       Indicates the type of boundary condition at PHI = PS and
C       PHI = PF.
C
C       = 0  If the solution is periodic in PHI, i.e.,
C            U(I,J) = U(I,N+J).
C       = 1  If the solution is specified at PHI = PS and PHI = PF
C            (see note below).
C       = 2  If the solution is specified at PHI = PS (see note below)
C            and the derivative of the solution with respect to PHI is
C            specified at PHI = PF.
C       = 3  If the derivative of the solution with respect to PHI is
C            specified at PHI = PS and PHI = PF.
C       = 4  If the derivative of the solution with respect to PHI is
C            specified at PS and the solution is specified at PHI = PF
C            (see note below).
C
C       NOTE:  NBDCND = 1,2, or 4 cannot be used with
C              MBDCND = 5,6,7,8, or 9 (the former indicates that the
C                       solution is specified at a pole, the latter
C                       indicates that the solution is unspecified).
C                       Use instead
C              MBDCND = 1 or 2  .
C
C     BDPS
C       A one-dimensional array of length M+1 that specifies the values
C       of the derivative of the solution with respect to PHI at
C       PHI = PS.  When NBDCND = 3 or 4,
C
C            BDPS(I) = (d/dPHI)U(THETA(I),PS), I = 1,2,...,M+1  .
C
C       When NBDCND has any other value, BDPS is a dummy variable.
C
C     BDPF
C       A one-dimensional array of length M+1 that specifies the values
C       of the derivative of the solution with respect to PHI at
C       PHI = PF.  When NBDCND = 2 or 3,
C
C            BDPF(I) = (d/dPHI)U(THETA(I),PF), I = 1,2,...,M+1  .
C
C       When NBDCND has any other value, BDPF 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, HWSSSP will
C       attempt to find a solution.
C
C     F
C       A two-dimensional array that specifies the value 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(THETA(I),PHI(J)).
C
C       On the boundaries F is defined by
C
C            MBDCND   F(1,J)            F(M+1,J)
C            ------   ------------      ------------
C
C              1      U(TS,PHI(J))      U(TF,PHI(J))
C              2      U(TS,PHI(J))      F(TF,PHI(J))
C              3      F(TS,PHI(J))      F(TF,PHI(J))
C              4      F(TS,PHI(J))      U(TF,PHI(J))
C              5      F(0,PS)           U(TF,PHI(J))   J = 1,2,...,N+1
C              6      F(0,PS)           F(TF,PHI(J))
C              7      U(TS,PHI(J))      F(PI,PS)
C              8      F(TS,PHI(J))      F(PI,PS)
C              9      F(0,PS)           F(PI,PS)
C
C            NBDCND   F(I,1)            F(I,N+1)
C            ------   --------------    --------------
C
C              0      F(THETA(I),PS)    F(THETA(I),PS)
C              1      U(THETA(I),PS)    U(THETA(I),PF)
C              2      U(THETA(I),PS)    F(THETA(I),PF)   I = 1,2,...,M+1
C              3      F(THETA(I),PS)    F(THETA(I),PF)
C              4      F(THETA(I),PS)    U(THETA(I),PF)
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
C     IDIMF
C       The row (or first) dimension of the array F as it appears in the
C       program calling HWSSSP.  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)+(16+INT(log2(N+1)))(M+1)
C       locations. The actual number of locations used is computed by
C       HWSSSP and is output in location W(1). INT( ) denotes the
C       FORTRAN integer function.
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 (THETA(I),PHI(J)),
C       I = 1,2,...,M+1,   J = 1,2,...,N+1  .
C
C     PERTRB
C       If one specifies a combination of periodic, derivative or
C       unspecified boundary conditions for a Poisson equation
C       (LAMBDA = 0), a solution may not exist.  PERTRB is a constant,
C       calculated and subtracted from F, which ensures that a solution
C       exists.  HWSSSP then computes this solution, which is a least
C       squares solution to the original approximation.  This solution
C       is not unique and is unnormalized. The value of PERTRB should
C       be small compared to the right side F. Otherwise , a solution
C       is obtained to an essentially different problem. This comparison
C       should always be made to insure that a meaningful solution has
C       been obtained.
C
C     IERROR
C       An error flag that indicates invalid input parameters.  Except
C       for numbers 0 and 8, a solution is not attempted.
C
C       = 0  No error
C       = 1  TS.LT.0 or TF.GT.PI
C       = 2  TS.GE.TF
C       = 3  MBDCND.LT.1 or MBDCND.GT.9
C       = 4  PS.LT.0 or PS.GT.PI+PI
C       = 5  PS.GE.PF
C       = 6  N.LT.5
C       = 7  M.LT.5
C       = 8  NBDCND.LT.0 or NBDCND.GT.4
C       = 9  ELMBDA.GT.0
C       = 10 IDIMF.LT.M+1
C       = 11 NBDCND equals 1,2 or 4 and MBDCND.GE.5
C       = 12 TS.EQ.0 and MBDCND equals 3,4 or 8
C       = 13 TF.EQ.PI and MBDCND equals 2,3 or 6
C       = 14 MBDCND equals 5,6 or 9 and TS.NE.0
C       = 15 MBDCND.GE.7 and TF.NE.PI
C
C       Since this is the only means of indicating a possibly incorrect
C       call to HWSSSP, the user should test IERROR after a call.
C
C     W
C       Contains intermediate values that must not be destroyed if
C       HWSSSP will be called again with INTL = 1. W(1) contains the
C       required length of W .
C
C *Long Description:
C
C     * * * * * * *   Program Specifications    * * * * * * * * * * * *
C
C     Dimension of   BDTS(N+1),BDTF(N+1),BDPS(M+1),BDPF(M+1),
C     Arguments      F(IDIMF,N+1),W(see argument list)
C
C     Latest         January 1978
C     Revision
C
C
C     Subprograms    HWSSSP,HWSSS1,GENBUN,POISD2,POISN2,POISP2,COSGEN,ME
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
C
C     Specialist     Paul Swarztrauber
C
C     Language       FORTRAN
C
C     History        Version 1 - September 1973
C                    Version 2 - April     1976
C                    Version 3 - January   1978
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
C     Required       CONTROL DATA 7600
C
C     Timing and        The execution time T on the NCAR Control Data
C     Accuracy       7600 for subroutine HWSSSP 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                    The machine dependent constant PI is defined in
C                    function PIMACH.
C
C     Required       SIN,COS
C     Resident
C     Routines
C
C     References     P. N. Swarztrauber,'The Direct Solution Of The
C                    Discrete Poisson Equation On The Surface Of a
C                    Sphere, SIAM J. Numer. Anal.,15(1974), pp 212-215
C
C                    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***REFERENCES  P. N. Swarztrauber and R. Sweet, Efficient Fortran
C                 subprograms for the solution of elliptic equations,
C                 NCAR TN/IA-109, July 1975, 138 pp.
C               P. N. Swarztrauber, The direct solution of the discrete
C                 Poisson equation on the surface of a sphere, SIAM
C                 Journal on Numerical Analysis 15 (1974), pp. 212-215.
C***ROUTINES CALLED  HWSSS1, PIMACH
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891009  Removed unreferenced variable.  (WRB)
C   891009  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  HWSSSP
C
      DIMENSION       F(IDIMF,*) ,BDTS(*)    ,BDTF(*)    ,BDPS(*)    ,
     1                BDPF(*)    ,W(*)
C***FIRST EXECUTABLE STATEMENT  HWSSSP
      PI = PIMACH(DUM)
      TPI = 2.*PI
      IERROR = 0
      IF (TS.LT.0. .OR. TF.GT.PI) IERROR = 1
      IF (TS .GE. TF) IERROR = 2
      IF (MBDCND.LT.1 .OR. MBDCND.GT.9) IERROR = 3
      IF (PS.LT.0. .OR. PF.GT.TPI) IERROR = 4
      IF (PS .GE. PF) IERROR = 5
      IF (N .LT. 5) IERROR = 6
      IF (M .LT. 5) IERROR = 7
      IF (NBDCND.LT.0 .OR. NBDCND.GT.4) IERROR = 8
      IF (ELMBDA .GT. 0.) IERROR = 9
      IF (IDIMF .LT. M+1) IERROR = 10
      IF ((NBDCND.EQ.1 .OR. NBDCND.EQ.2 .OR. NBDCND.EQ.4) .AND.
     1    MBDCND.GE.5) IERROR = 11
      IF (TS.EQ.0. .AND.
     1    (MBDCND.EQ.3 .OR. MBDCND.EQ.4 .OR. MBDCND.EQ.8)) IERROR = 12
      IF (TF.EQ.PI .AND.
     1    (MBDCND.EQ.2 .OR. MBDCND.EQ.3 .OR. MBDCND.EQ.6)) IERROR = 13
      IF ((MBDCND.EQ.5 .OR. MBDCND.EQ.6 .OR. MBDCND.EQ.9) .AND.
     1    TS.NE.0.) IERROR = 14
      IF (MBDCND.GE.7 .AND. TF.NE.PI) IERROR = 15
      IF (IERROR.NE.0 .AND. IERROR.NE.9) RETURN
      CALL HWSSS1 (TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS,BDPF,
     1             ELMBDA,F,IDIMF,PERTRB,W,W(M+2),W(2*M+3),W(3*M+4),
     2             W(4*M+5),W(5*M+6),W(6*M+7))
      W(1) = W(6*M+7)+6*(M+1)
      RETURN
      END