Try out NAG Library functions

Explore NAG maths and stats routines with interactive demos
Function ID
E02BEF
Name
nagf_fit_1dspline_auto
Description
Least squares cubic spline curve fit, automatic knot placement, one variable
Keywords
B-splines | cubic splines | least squares | smoothing
This example reads in a set of data values, followed by a set of values of S. For each value of S it calls E02BEF to compute a spline approximation, and prints the values of the knots and the B-spline coefficients ci.
The example data reflects that shown in the "Example" section of the routine documentation. You can change this here to try alternative inputs. The formatting will need to be kept as it is here, otherwise the program is likely to fail to run correctly.

Please note that incompatible data will however cause the example output to display an error message. These error messages are fully explained in the Routine document
    Program e02befe

!     E02BEF Example Program Text

!     Mark 26 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: e02bbf, e02bef, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: fp, s, txr
      Integer                          :: ifail, ioerr, j, lwrk, m, n, nest, r
      Character (1)                    :: start
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: c(:), lamda(:), sp(:), w(:), wrk(:), &
                                          x(:), y(:)
      Integer, Allocatable             :: iwrk(:)
!     .. Executable Statements ..
      Write (nout,*) 'E02BEF Example Program Results'

!     Skip heading in data file
      Read (nin,*)

!     Input the number of data points, followed by the data points (X),
!     the function values (Y) and the weights (W).

      Read (nin,*) m
      nest = m + 4
      lwrk = 4*m + 16*nest + 41
      Allocate (x(m),y(m),w(m),iwrk(nest),lamda(nest),wrk(lwrk),c(nest),       &
        sp(2*m-1))

      Do r = 1, m
        Read (nin,*) x(r), y(r), w(r)
      End Do

      start = 'C'

!     Read in successive values of S until end of data file.

data: Do
        Read (nin,*,Iostat=ioerr) s

        If (ioerr<0) Then
          Exit data
        End If

!       Determine the spline approximation.

        ifail = 0
        Call e02bef(start,m,x,y,w,s,nest,n,lamda,c,fp,wrk,lwrk,iwrk,ifail)

!       Evaluate the spline at each X point and midway between
!       X points, saving the results in SP.

        Do r = 1, m

          ifail = 0
          Call e02bbf(n,lamda,c,x(r),sp((r-1)*2+1),ifail)

        End Do

        Do r = 1, m - 1
          txr = (x(r)+x(r+1))/2.0E0_nag_wp

          ifail = 0
          Call e02bbf(n,lamda,c,txr,sp(r*2),ifail)

        End Do

!       Output the results.

        Write (nout,*)
        Write (nout,99999) 'Calling with smoothing factor S =', s
        Write (nout,*)
        Write (nout,*) '                                           B-Spline'
        Write (nout,*)                                                         &
          '             J       Knot LAMDA(J+2)   Coefficient C(J)'
        Write (nout,99998) 1, c(1)

        Do j = 2, n - 5
          Write (nout,99997) j, lamda(j+2), c(j)
        End Do

        Write (nout,99998) n - 4, c(n-4)
        Write (nout,*)
        Write (nout,99999) 'Weighted sum of squared residuals FP =', fp

        If (fp==0.0E0_nag_wp) Then
          Write (nout,*) '(The spline is an interpolating spline)'
        Else If (n==8) Then
          Write (nout,*)                                                       &
            '(The spline is the weighted least squares cubic polynomial)'
        End If

        Write (nout,*)
        start = 'W'
      End Do data

99999 Format (1X,A,1P,E12.3)
99998 Format (11X,I4,20X,F16.4)
99997 Format (11X,I4,2F18.4)
    End Program e02befe
The NAG Library
The world’s largest collection of robust, documented, tested and maintained numerical algorithms.
Learn more here or contact us for purchasing information