cubgcv Subroutine

  subroutine cubgcv(x,f,df,n,y,c,ic,var,job,se,wk,ier)
!
!---specifications for arguments---
    integer, intent(in) :: n,ic,job
    integer, intent(out) :: ier
!
! All real variables should be double precision:
!
    real, dimension(:), intent(in) :: x, f
    real, dimension(:), intent(in out) :: df
    real, dimension(:), intent(out) :: y, se
    real, dimension(:, :), intent(out) :: c
    real, dimension(0:n+1, 7), intent(out) :: wk
    real, intent(in out) :: var
!
!---specifications for local variables---
    real :: delta,err,gf1,gf2,gf3,gf4,r1,r2,r3,r4,avh,avdf,avar,stat(6),p,q
    real, parameter :: ratio = 2., tau = 1.618033989, zero = 0., one = 1.
    integer :: i
!
!---initialize---
    ier = 133
    if (job.lt.0 .or. job.gt.1) go to 140
    call spint1(x,avh,f,df,avdf,n,y,c,ic,wk,wk(0,4),ier)
    if (ier.ne.0) go to 140
    avar = var
    if (var.gt.zero) avar = var*avdf*avdf
!
!---check for zero variance---
    if (var.ne.zero) go to 10
    r1 = zero
    go to 90
!
!---find local minimum of gcv or the expected mean square error---
10  r1 = one
    r2 = ratio*r1
    call spfit1(x,avh,df,n,r2,p,q,gf2,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
20  call spfit1(x,avh,df,n,r1,p,q,gf1,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
    if (gf1.gt.gf2) go to 30
!
!---exit if p zero---
    if (p.le.zero) go to 100
    r2 = r1
    gf2 = gf1
    r1 = r1/ratio
    go to 20

30  r3 = ratio*r2
40  call spfit1(x,avh,df,n,r3,p,q,gf3,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
    if (gf3.gt.gf2) go to 50
!
!---exit if q zero---
    if (q.le.zero) go to 100
    r2 = r3
    gf2 = gf3
    r3 = ratio*r3
    go to 40

50  r2 = r3
    gf2 = gf3
    delta = (r2-r1)/tau
    r4 = r1 + delta
    r3 = r2 - delta
    call spfit1(x,avh,df,n,r3,p,q,gf3,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
    call spfit1(x,avh,df,n,r4,p,q,gf4,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
!
!---golden section search for local minimum---
60  if (gf3.gt.gf4) go to 70
    r2 = r4
    gf2 = gf4
    r4 = r3
    gf4 = gf3
    delta = delta/tau
    r3 = r2 - delta
    call spfit1(x,avh,df,n,r3,p,q,gf3,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
    go to 80

70  r1 = r3
    gf1 = gf3
    r3 = r4
    gf3 = gf4
    delta = delta/tau
    r4 = r1 + delta
    call spfit1(x,avh,df,n,r4,p,q,gf4,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
80  err = (r2-r1)/ (r1+r2)
    if (err*err+one.gt.one .and. err.gt.1.0e-6) go to 60
    r1 = (r1+r2)*0.5
!
!---calculate spline coefficients---
90  call spfit1(x,avh,df,n,r1,p,q,gf1,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
100 call spcof1(x,avh,f,df,n,p,q,y,c,ic,wk(:,6),wk(:,7))
!
!---optionally calculate standard error estimates---
    if (var.ge.zero) go to 110
    avar = stat(6)
    var = avar/ (avdf*avdf)
110 if (job.eq.1) call sperr1(x,avh,df,n,wk,p,avar,se)
!
!---unscale df---
    df = df * avdf
!
!--put statistics in wk---
    do i = 0,5
       wk(i,1) = stat(i+1)
    end do
    wk(5,1) = stat(6)/ (avdf*avdf)
    wk(6,1) = avdf*avdf
    go to 150
!
!---check for error condition---
140 continue
!     if (ier.ne.0) continue
150 return
  end subroutine cubgcv