FIXME : Add documentation
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real, | intent(in), | dimension(:) | :: | x | ||
real, | intent(in) | :: | avh | |||
real, | intent(in), | dimension(:) | :: | dy | ||
integer, | intent(in) | :: | n | |||
real, | intent(in) | :: | rho | |||
real, | intent(out) | :: | p | |||
real, | intent(out) | :: | q | |||
real, | intent(out) | :: | fun | |||
real, | intent(in) | :: | var | |||
real, | intent(out), | dimension(6) | :: | stat | ||
real, | intent(in), | dimension(:) | :: | a | ||
real, | intent(in), | dimension(ic, 3) | :: | c | ||
integer, | intent(in) | :: | ic | |||
real, | intent(inout), | dimension(0:n+1, 3) | :: | r | ||
real, | intent(in), | dimension(0:n+1, 2) | :: | t | ||
real, | intent(out), | dimension(0:n+1) | :: | u | ||
real, | intent(out), | dimension(0:n+1) | :: | v |
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, intent(in) :: ic,n
! all variables should be double precision:
real, dimension(:), intent(in) :: x, dy, a
real, dimension(6), intent(out) :: stat
real, dimension(ic, 3), intent(in) :: c
real, dimension(0:n+1, 3), intent(in out) :: r
real, dimension(0:n+1, 2), intent(in) :: t
real, dimension(0:n+1), intent(out) :: u, v
real, intent(in) :: var, avh, rho
real, intent(out) :: fun, p, q
!
!---local variables---
integer :: i
real :: e,f,g,h,rho1
real, parameter :: zero = 0., one = 1., 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==one) p = zero
if (rho1==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>=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