relibc_openlibm/slatec/dpolvl.f

*DECK DPOLVL
      SUBROUTINE DPOLVL (NDER, XX, YFIT, YP, N, X, C, WORK, IERR)
C***BEGIN PROLOGUE  DPOLVL
C***PURPOSE  Calculate the value of a polynomial and its first NDER
C            derivatives where the polynomial was produced by a previous
C            call to DPLINT.
C***LIBRARY   SLATEC
C***CATEGORY  E3
C***TYPE      DOUBLE PRECISION (POLYVL-S, DPOLVL-D)
C***KEYWORDS  POLYNOMIAL EVALUATION
C***AUTHOR  Huddleston, R. E., (SNLL)
C***DESCRIPTION
C
C     Abstract -
C        Subroutine DPOLVL calculates the value of the polynomial and
C     its first NDER derivatives where the polynomial was produced by
C     a previous call to DPLINT.
C        The variable N and the arrays X and C must not be altered
C     between the call to DPLINT and the call to DPOLVL.
C
C     ******  Dimensioning Information *******
C
C     YP   must be dimensioned by at least NDER
C     X    must be dimensioned by at least N (see the abstract )
C     C    must be dimensioned by at least N (see the abstract )
C     WORK must be dimensioned by at least 2*N if NDER is .GT. 0.
C
C     *** Note ***
C       If NDER=0, neither YP nor WORK need to be dimensioned variables.
C       If NDER=1, YP does not need to be a dimensioned variable.
C
C
C     *****  Input parameters
C       ***  All TYPE REAL variables are DOUBLE PRECISION ***
C
C     NDER - the number of derivatives to be evaluated
C
C     XX   - the argument at which the polynomial and its derivatives
C            are to be evaluated.
C
C     N    - *****
C            *       N, X, and C must not be altered between the call
C     X    - *       to DPLINT and the call to DPOLVL.
C     C    - *****
C
C
C     *****  Output Parameters
C       ***  All TYPE REAL variables are DOUBLE PRECISION ***
C
C     YFIT - the value of the polynomial at XX
C
C     YP   - the derivatives of the polynomial at XX.  The derivative of
C            order J at XX is stored in  YP(J) , J = 1,...,NDER.
C
C     IERR - Output error flag with the following possible values.
C          = 1  indicates normal execution
C
C     ***** Storage Parameters
C
C     WORK  = this is an array to provide internal working storage for
C             DPOLVL.  It must be dimensioned by at least 2*N if NDER is
C             .GT. 0.  If NDER=0, WORK does not need to be a dimensioned
C             variable.
C
C***REFERENCES  L. F. Shampine, S. M. Davenport and R. E. Huddleston,
C                 Curve fitting by polynomials in one variable, Report
C                 SLA-74-0270, Sandia Laboratories, June 1974.
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   740601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891006  Cosmetic changes to prologue.  (WRB)
C   891006  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  DPOLVL
      INTEGER I,IERR,IM1,IZERO,K,KM1,KM1PI,KM2PN,KM2PNI,M,MM,N,NDR,NDER,
     * NMKP1,NPKM1
      DOUBLE PRECISION  C(*),FAC,PIONE,PITWO,PONE,PTWO,X(*),XK,XX,
     * YFIT,YP(*),WORK(*)
C***FIRST EXECUTABLE STATEMENT  DPOLVL
      IERR=1
         IF (NDER.GT.0) GO TO 10020
C
C     *****   CODING FOR THE CASE NDER = 0
C
      PIONE=1.0D0
      PONE=C(1)
      YFIT=PONE
      IF (N.EQ.1) RETURN
      DO 10010 K=2,N
      PITWO=(XX-X(K-1))*PIONE
      PIONE=PITWO
      PTWO=PONE+PITWO*C(K)
      PONE=PTWO
10010 CONTINUE
      YFIT=PTWO
      RETURN
C
C     *****   END OF NDER = 0 CASE
C
10020 CONTINUE
         IF (N.GT.1) GO TO 10040
      YFIT=C(1)
C
C     *****  CODING FOR THE CASE  N=1 AND NDER .GT. 0
C
      DO 10030 K=1,NDER
      YP(K)=0.0D0
10030 CONTINUE
      RETURN
C
C     *****  END OF THE CASE  N = 1 AND  NDER .GT. 0
C
10040 CONTINUE
         IF (NDER.LT.N) GO TO 10050
C
C     *****  SET FLAGS FOR NUMBER OF DERIVATIVES AND FOR DERIVATIVES
C            IN EXCESS OF THE DEGREE (N-1) OF THE POLYNOMIAL.
C
      IZERO=1
      NDR=N-1
         GO TO 10060
10050 CONTINUE
      IZERO=0
      NDR=NDER
10060 CONTINUE
      M=NDR+1
      MM=M
C
C     *****  START OF THE CASE NDER .GT. 0  AND N .GT. 1
C     *****  THE POLYNOMIAL AND ITS DERIVATIVES WILL BE EVALUATED AT XX
C
      DO 10070 K=1,NDR
      YP(K)=C(K+1)
10070 CONTINUE
C
C     *****  THE FOLLOWING SECTION OF CODE IS EASIER TO READ IF ONE
C            BREAKS WORK INTO TWO ARRAYS W AND V. THE CODE WOULD THEN
C            READ
C                W(1) = 1.
C                PONE = C(1)
C               *DO   K = 2,N
C               *   V(K-1) =  XX - X(K-1)
C               *   W(K)   =  V(K-1)*W(K-1)
C               *   PTWO   =  PONE + W(K)*C(K)
C               *   PONE   =  PWO
C
C               YFIT = PTWO
C
      WORK(1)=1.0D0
      PONE=C(1)
      DO 10080 K=2,N
      KM1=K-1
      NPKM1=N+K-1
      WORK(NPKM1)=XX-X(KM1)
      WORK(K)=WORK(NPKM1)*WORK(KM1)
      PTWO=PONE+WORK(K)*C(K)
      PONE=PTWO
10080 CONTINUE
      YFIT=PTWO
C
C     ** AT THIS POINT THE POLYNOMIAL HAS BEEN EVALUATED AND INFORMATION
C        FOR THE DERIVATIVE EVALUATIONS HAVE BEEN STORED IN THE ARRAY
C        WORK
         IF (N.EQ.2) GO TO 10110
      IF (M.EQ.N) MM=NDR
C
C     ***** EVALUATE THE DERIVATIVES AT XX
C
C                  ******  DO K=2,MM   (FOR MOST CASES, MM = NDER + 1)
C                  *  ******  DO I=2,N-K+1
C                  *  *       W(I) = V(K-2+I)*W(I-1) + W(I)
C                  *  *       YP(K-1) = YP(K-1) + W(I)*C(K-1+I)
C                  ******  CONTINUE
C
      DO 10090 K=2,MM
      NMKP1=N-K+1
      KM1=K-1
      KM2PN=K-2+N
      DO 10090 I=2,NMKP1
      KM2PNI=KM2PN+I
      IM1=I-1
      KM1PI=KM1+I
      WORK(I)=WORK(KM2PNI)*WORK(IM1)+WORK(I)
      YP(KM1)=YP(KM1)+WORK(I)*C(KM1PI)
10090 CONTINUE
         IF (NDR.EQ.1) GO TO 10110
      FAC=1.0D0
      DO 10100 K=2,NDR
      XK=K
      FAC=XK*FAC
      YP(K)=FAC*YP(K)
10100 CONTINUE
C
C     ***** END OF DERIVATIVE EVALUATIONS
C
10110 CONTINUE
      IF (IZERO.EQ.0) RETURN
C
C     *****  SET EXCESS DERIVATIVES TO ZERO.
C
      DO 10120 K=N,NDER
      YP(K)=0.0D0
10120 CONTINUE
      RETURN
      END