polint.f 2.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162
  1. *DECK POLINT
  2. SUBROUTINE POLINT (N, X, Y, C)
  3. C***BEGIN PROLOGUE POLINT
  4. C***PURPOSE Produce the polynomial which interpolates a set of discrete
  5. C data points.
  6. C***LIBRARY SLATEC
  7. C***CATEGORY E1B
  8. C***TYPE SINGLE PRECISION (POLINT-S, DPLINT-D)
  9. C***KEYWORDS POLYNOMIAL INTERPOLATION
  10. C***AUTHOR Huddleston, R. E., (SNLL)
  11. C***DESCRIPTION
  12. C
  13. C Written by Robert E. Huddleston, Sandia Laboratories, Livermore
  14. C
  15. C Abstract
  16. C Subroutine POLINT is designed to produce the polynomial which
  17. C interpolates the data (X(I),Y(I)), I=1,...,N. POLINT sets up
  18. C information in the array C which can be used by subroutine POLYVL
  19. C to evaluate the polynomial and its derivatives and by subroutine
  20. C POLCOF to produce the coefficients.
  21. C
  22. C Formal Parameters
  23. C N - the number of data points (N .GE. 1)
  24. C X - the array of abscissas (all of which must be distinct)
  25. C Y - the array of ordinates
  26. C C - an array of information used by subroutines
  27. C ******* Dimensioning Information *******
  28. C Arrays X,Y, and C must be dimensioned at least N in the calling
  29. C program.
  30. C
  31. C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston,
  32. C Curve fitting by polynomials in one variable, Report
  33. C SLA-74-0270, Sandia Laboratories, June 1974.
  34. C***ROUTINES CALLED XERMSG
  35. C***REVISION HISTORY (YYMMDD)
  36. C 740601 DATE WRITTEN
  37. C 861211 REVISION DATE from Version 3.2
  38. C 891214 Prologue converted to Version 4.0 format. (BAB)
  39. C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
  40. C 920501 Reformatted the REFERENCES section. (WRB)
  41. C***END PROLOGUE POLINT
  42. DIMENSION X(*),Y(*),C(*)
  43. C***FIRST EXECUTABLE STATEMENT POLINT
  44. IF (N .LE. 0) GO TO 91
  45. C(1)=Y(1)
  46. IF(N .EQ. 1) RETURN
  47. DO 10010 K=2,N
  48. C(K)=Y(K)
  49. KM1=K-1
  50. DO 10010 I=1,KM1
  51. C CHECK FOR DISTINCT X VALUES
  52. DIF = X(I)-X(K)
  53. IF (DIF .EQ. 0.0) GO TO 92
  54. C(K) = (C(I)-C(K))/DIF
  55. 10010 CONTINUE
  56. RETURN
  57. 91 CALL XERMSG ('SLATEC', 'POLINT', 'N IS ZERO OR NEGATIVE.', 2, 1)
  58. RETURN
  59. 92 CALL XERMSG ('SLATEC', 'POLINT',
  60. + 'THE ABSCISSAS ARE NOT DISTINCT.', 2, 1)
  61. RETURN
  62. END