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