1
2#:::::::::::
3#   dgold
4#:::::::::::
5
6subroutine  dgold (vmu, q, ldq, n, z, low, upp, nlaht, score, varht,_
7                   info, twk, work)
8
9#  Purpose:  To evaluate GCV/GML function based on tridiagonal form and to
10#      search minimum on an interval by golden section search.
11
12character*1       vmu
13integer           ldq, n, info
14double precision  q(ldq,*), z(*), low, upp, nlaht, score, varht, twk(2,*),_
15                  work(*)
16
17#  On entry:
18#      vmu        'v':  GCV criterion.
19#                 'm':  GML criterion.
20#                 'u':  unbiased risk estimate.
21#      q          tidiagonal matrix in diagonal and super diagonal.
22#      ldq        leading dimension of Q.
23#      n          size of the matrix.
24#      z          U^{T} F_{2}^{T} y.
25#      low        lower limit of log10(n*lambda).
26#      upp        upper limit of log10(n*lambda).
27#      varht      known variance if vmu=='u'.
28
29#  On exit:
30#      nlaht      the estimated log(n*lambda).
31#      score      the GCV/GML/URE score at the estimated lambda.
32#      varht      the variance estimate at the estimated lambda.
33#      info        0: normal termination.
34#                 -1: dimension error.
35#                 -2: tridiagonal form is not non-negative definite.
36#                 -3: vmu is none of 'v', 'm', or 'u'.
37
38#  Work arrays:
39#      twk        of size at least (2,n).
40#      work       of size at least (n).
41
42#  Routines called directly:
43#      Fortran -- dsqrt
44#      Blas    -- daxpy, dcopy
45#      Rkpack  -- dtrev
46#      Other   -- dset
47
48#  Written:  Chong Gu, Statistics, Purdue, latest version 12/29/91.
49
50double precision  ratio, mlo, mup, tmpl, tmpu
51
52ratio = ( dsqrt (5.d0) - 1.d0 ) / 2.d0
53
54info = 0
55
56#   interchange the boundaries if necessary
57if ( upp < low ) {
58    mlo = low
59    low = upp
60    upp = mlo
61}
62
63#   check vmu
64if ( vmu != 'v' & vmu != 'm' & vmu != 'u' ) {
65    info = -3
66    return
67}
68
69#   check dimension
70if ( n < 1 | n > ldq ) {
71    info = -1
72    return
73}
74
75#   initialize golden section search for scrht
76mlo = upp - ratio * (upp - low)
77call  dset (n, 10.d0 ** (mlo), twk(2,1), 2)
78call  daxpy (n, 1.d0, q, ldq+1, twk(2,1), 2)
79call  dcopy (n-1, q(1,2), ldq+1, twk(1,2), 2)
80twk(1,1) = 10.d0**mlo
81call  dtrev (vmu, twk, 2, n, z, tmpl, varht, info, work)
82if ( info != 0 ) {
83    info = -2
84    return
85}
86mup = low + ratio * (upp - low)
87call  dset (n, 10.d0 ** (mup), twk(2,1), 2)
88call  daxpy (n, 1.d0, q, ldq+1, twk(2,1), 2)
89call  dcopy (n-1, q(1,2), ldq+1, twk(1,2), 2)
90twk(1,1) = 10.d0**mup
91call  dtrev (vmu, twk, 2, n, z, tmpu, varht, info, work)
92if ( info != 0 ) {
93    info = -2
94    return
95}
96
97#   golden section search for estimate of lambda
98repeat {
99    if ( mup - mlo < 1.d-7 )  break
100    if ( tmpl < tmpu ) {
101        upp = mup
102        mup = mlo
103        tmpu = tmpl
104        mlo = upp - ratio * (upp - low)
105        call  dset (n, 10.d0 ** (mlo), twk(2,1), 2)
106        call  daxpy (n, 1.d0, q, ldq+1, twk(2,1), 2)
107        call  dcopy (n-1, q(1,2), ldq+1, twk(1,2), 2)
108        twk(1,1) = 10.d0**mlo
109        call  dtrev (vmu, twk, 2, n, z, tmpl, varht, info, work)
110        if ( info != 0 ) {
111            info = -2
112            return
113        }
114    }
115    else {
116        low = mlo
117        mlo = mup
118        tmpl = tmpu
119        mup = low + ratio * (upp - low)
120        call  dset (n, 10.d0 ** (mup), twk(2,1), 2)
121        call  daxpy (n, 1.d0, q, ldq+1, twk(2,1), 2)
122        call  dcopy (n-1, q(1,2), ldq+1, twk(1,2), 2)
123        twk(1,1) = 10.d0**mup
124        call  dtrev (vmu, twk, 2, n, z, tmpu, varht, info, work)
125        if ( info != 0 ) {
126            info = -2
127            return
128        }
129    }
130}
131
132#   compute the return value
133nlaht = ( mup + mlo ) / 2.d0
134call  dset (n, 10.d0 ** (nlaht), twk(2,1), 2)
135call  daxpy (n, 1.d0, q, ldq+1, twk(2,1), 2)
136call  dcopy (n-1, q(1,2), ldq+1, twk(1,2), 2)
137twk(1,1) = 10.d0**nlaht
138call  dtrev (vmu, twk, 2, n, z, score, varht, info, work)
139if ( info != 0 ) {
140    info = -2
141    return
142}
143
144return
145end
146
147#...............................................................................
148
149