spfit1 Subroutine

  subroutine spfit1(x,avh,dy,n,rho,p,q,fun,var,stat,a,c,ic,r,t,u,v)
!
! fits a cubic smoothing spline to data with relative
! weighting dy for a given value of the smoothing parameter
! rho using an algorithm based on that of c.h. reinsch (1967),
! numer. math. 10, 177-183.
!
! the trace of the influence matrix is calculated using an
! algorithm developed by m.f.hutchinson and f.r.de hoog (numer.
! math., in press), enabling the generalized cross validation
! and related statistics to be calculated in order n operations.
!
! the arrays a, c, r and t are assumed to have been initialized
! by the subroutine spint1.  overflow and underflow problems are
! avoided by using p=rho/(1 + rho) and q=1/(1 + rho) instead of
! rho and by scaling the differences x(i+1) - x(i) by avh.
!
! the values in df are assumed to have been scaled so that the
! sum of their squared values is n.  the value in var, when it is
! non-negative, is assumed to have been scaled to compensate for
! the scaling of the values in df.
!
! the value returned in fun is an estimate of the true mean square
! when var is non-negative, and is the generalized cross validation
! when var is negative.
!
!---specifications for arguments---
    integer ic,n
! all variables should be double precision:
    real, dimension (:) :: x, dy, a
    real, dimension (6) :: stat
    real, dimension (ic,3) :: c
    real, dimension (0:n+1,3) :: r
    real, dimension (0:n+1,2) :: t
    real, dimension (0:n+1) :: u, v
    real :: fun, var, avh, p, q, rho
!
!---local variables---
    integer i
    real :: e,f,g,h,rho1
    real :: zero = 0.
    real :: one = 1.
    real :: two = 2.

!
!---use p and q instead of rho to prevent overflow or underflow---
    rho1 = one + rho
    p = rho/rho1
    q = one/rho1
    if (rho1.eq.one) p = zero
    if (rho1.eq.rho) q = zero
!
!---rational cholesky decomposition of p*c + q*t---
    f = zero
    g = zero
    h = zero
    do i = 0,1
       r(i,1) = zero
    end do
    do i = 2,n - 1
       r(i-2,3) = g*r(i-2,1)
       r(i-1,2) = f*r(i-1,1)
       r(i,1) = one/ (p*c(i,1)+q*t(i,1)-f*r(i-1,2)-g*r(i-2,3))
       f = p*c(i,2) + q*t(i,2) - h*r(i-1,2)
       g = h
       h = p*c(i,3)
    end do
!
!---solve for u---
    u(0) = zero
    u(1) = zero
    do i = 2,n - 1
       u(i) = a(i) - r(i-1,2)*u(i-1) - r(i-2,3)*u(i-2)
    end do
    u(n) = zero
    u(n+1) = zero
    do i = n - 1,2,-1
       u(i) = r(i,1)*u(i) - r(i,2)*u(i+1) - r(i,3)*u(i+2)
    end do
!
!---calculate residual vector v---
    e = zero
    h = zero
    do i = 1,n - 1
       g = h
       h = (u(i+1)-u(i))/ ((x(i+1)-x(i))/avh)
       v(i) = dy(i)* (h-g)
       e = e + v(i)*v(i)
    end do
    v(n) = dy(n)* (-h)
    e = e + v(n)*v(n)
!
!---calculate upper three bands of inverse matrix---
    r(n,1) = zero
    r(n,2) = zero
    r(n+1,1) = zero
    do i = n - 1,2,-1
       g = r(i,2)
       h = r(i,3)
       r(i,2) = -g*r(i+1,1) - h*r(i+1,2)
       r(i,3) = -g*r(i+1,2) - h*r(i+2,1)
       r(i,1) = r(i,1) - g*r(i,2) - h*r(i,3)
    end do
!
!---calculate trace---
    f = zero
    g = zero
    h = zero
    do i = 2,n - 1
       f = f + r(i,1)*c(i,1)
       g = g + r(i,2)*c(i,2)
       h = h + r(i,3)*c(i,3)
    end do
    f = f + two* (g+h)
!
!---calculate statistics---
    stat(1) = p
    stat(2) = f*p
    stat(3) = n*e/ (f*f)
    stat(4) = e*p*p/n
    stat(6) = e*p/f
    if (var.ge.zero) go to 80
    stat(5) = stat(6) - stat(4)
    fun = stat(3)
    go to 90

80  stat(5) = amax1(stat(4)-two*var*stat(2)/n+var,zero)
    fun = stat(5)
90  return
  end subroutine spfit1