Directory of image this file is from
This file as a plain text file

C
C     ..................................................................
C
C        SUBROUTINE DSG13
C
C        PURPOSE
C           TO COMPUTE A VECTOR OF SMOOTHED FUNCTION VALUES GIVEN
C           VECTORS OF ARGUMENT VALUES AND CORRESPONDING FUNCTION
C           VALUES.
C
C        USAGE
C           CALL DSG13(X,Y,Z,NDIM,IER)
C
C        DESCRIPTION OF PARAMETERS
C           X     -  GIVEN VECTOR OF DOUBLE PRECISION ARGUMENT VALUES
C                    (DIMENSION NDIM)
C           Y     -  GIVEN VECTOR OF DOUBLE PRECISION FUNCTION VALUES
C                    CORRESPONDING TO X (DIMENSION NDIM)
C           Z     -  RESULTING VECTOR OF DOUBLE PRECISION SMOOTHED
C                    FUNCTION VALUES (DIMENSION NDIM)
C           NDIM  -  DIMENSION OF VECTORS X,Y,AND Z
C           IER   -  RESULTING ERROR PARAMETER
C                    IER = -1  - NDIM IS LESS THAN 3
C                    IER =  0  - NO ERROR
C
C        REMARKS
C           (1)  IF IER=-1 THERE HAS BEEN NO COMPUTATION.
C           (2)   Z CAN HAVE THE SAME STORAGE ALLOCATION AS X OR Y.  IF
C                 X OR Y IS DISTINCT FROM Z, THEN IT IS NOT DESTROYED.
C
C        SUBROUTINES AND SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           EXCEPT AT THE ENDPOINTS X(1) AND X(NDIM), EACH SMOOTHED
C           VALUE Z(I) IS OBTAINED BY EVALUATING AT X(I) THE LEAST-
C           SQUARES POLYNOMIAL OF DEGREE 1 RELEVANT TO THE 3 SUCCESSIVE
C           POINTS (X(I+K),Y(I+K)) K = -1,0,1.(SEE HILDEBRAND, F.B.,
C           INTRODUCTION  TO NUMERICAL ANALYSIS, MC GRAW-HILL, NEW YORK/
C           TORONTO/LONDON, 1956, PP.258-311.)
C
C     ..................................................................
C
      SUBROUTINE DSG13(X,Y,Z,NDIM,IER)
C
C
      DIMENSION X(1),Y(1),Z(1)
      DOUBLE PRECISION X,Y,Z,XM,YM,T1,T2,T3,H
C
C        TEST OF DIMENSION
      IF(NDIM-3)7,1,1
C
C        START LOOP
    1 DO 6 I=3,NDIM
      XM=.33333333333333333D0*(X(I-2)+X(I-1)+X(I))
      YM=.33333333333333333D0*(Y(I-2)+Y(I-1)+Y(I))
      T1=X(I-2)-XM
      T2=X(I-1)-XM
      T3=X(I)-XM
      XM=T1*T1+T2*T2+T3*T3
      IF(XM)3,3,2
    2 XM=(T1*(Y(I-2)-YM)+T2*(Y(I-1)-YM)+T3*(Y(I)-YM))/XM
C
C        CHECK FIRST POINT
    3 IF(I-3)4,4,5
    4 H=XM*T1+YM
    5 Z(I-2)=H
    6 H=XM*T2+YM
C        END OF LOOP
C
C        UPDATE LAST TWO COMPONENTS
      Z(NDIM-1)=H
      Z(NDIM)=XM*T3+YM
      IER=0
      RETURN
C
C        ERROR EXIT IN CASE NDIM IS LESS THAN 3
    7 IER=-1
      RETURN
      END