Nothing
#:::::::::::
# ddeev
#:::::::::::
subroutine ddeev (vmu, nobs,_
q, ldqr, ldqc, n, nq, u, ldu, uaux, t, x,_ # inputs
theta, nlaht, score, varht,_
hes, ldh, gra,_ # outputs
hwk1, hwk2, gwk1, gwk2,_ # work arrays
kwk, ldk, work1, work2, work3,_
info)
# Acronym: Double precision DErivative EValuation.
# Purpose: This routine calculates the gradient and the Hessian of
# V(theta|lambda) or M(theta|lambda).
character*1 vmu
integer nobs, ldqr, ldqc, n, nq, ldu, ldh, ldk, info
double precision q(ldqr,ldqc,*), u(ldu,*), uaux(*), t(2,*), x(*),_
theta(*), nlaht, score, varht,_
hes(ldh,*), gra(*), hwk1(nq,*), hwk2(nq,*), gwk1(*), gwk2(*),_
kwk(ldk,ldk,*), work1(*), work2(*), work3(*)
# On entry:
# vmu 'v': GCV criterion.
# 'm': GML criterion.
# 'u': unbiased risk estimate.
# nobs the number of observations.
# q F_{2}^{T} Q_{i} F_{2}, of size (n,n,nq).
# n the size of q.
# nq the number of Q_{i}'s.
# u,uaux Householder vectors of U, of size (n-1,n-2),
# where U^{T}DU is tridiagonal.
# t U^{T} (D-n\lambda I) U in packed form, of size (2,n).
# x U^{T}z = U^{T}F_{2}^{T}y, of size (n).
# theta the current log(theta) for the D matrix, of dimension (nq).
# nlaht the estimated log10(n*lambda) in the current model.
# score the minimum GCV/GML score found at (theta, nlaht).
# varht the variance estimate at (theta, nlaht).
# On exit:
# hes Hessian at point (theta, nlaht), of size (nq,nq).
# gra gradient at point (theta, nlaht), of size (nq).
# info 0 : normal termination.
# -1 : dimension error.
# -2 : D !>= 0.
# -3 : tuning parameters are out of scope.
# Work arrays:
# hwk1,2 of sizes at least (nq,nq).
# gwk1,2 of sizes at least (nq).
# kwk of size at least (n,n,nq).
# work1-3 of sizes at least (n).
# Routines called directly:
# Fortran -- dble
# Blas -- daxpy, dcopy, ddot, dscal
# Blas2 -- dgemv
# Linpack -- dpbfa, dpbsl, dqrsl
# Rkpack -- dqrslm
# Other -- dset
# Written: Chong Gu, Statistics, Purdue, latest version 12/29/91.
double precision trc, det, dum, ddot
integer i, j, m
info = 0
call dset (nq, 0.d0, gra, 1)
call dset (nq*nq, 0.d0, hes, 1)
# check tuning parameters
if ( vmu != 'v' & vmu != 'm' & vmu != 'u' ) {
info = -3
return
}
# check dimension
if ( nobs < n | ldqr < n | ldqc < n | nq <= 0 | ldu < n-1 | ldh < nq | ldk < n ) {
info = -1
return
}
# compute K_{i} = U^{T}(\theta_{i}Q_{i})U
for (i=2;i<=nq;i=i+1) {
# from i=2 to nq
if ( theta(i) <= -25.d0 ) next
for (j=1;j<=n;j=j+1) {
call dcopy (n-j+1, q(j,j,i), 1, kwk(j,j,i), 1)
call dscal (n-j+1, 10.d0 ** theta(i), kwk(j,j,i), 1)
}
call dqrslm (u, ldu, n-1, n-2, uaux, kwk(2,2,i), n, 0, info, work1)
call dqrsl (u, ldu, n-1, n-2, uaux, kwk(2,1,i), dum, kwk(2,1,i),_
dum, dum, dum, 01000, info)
}
# compute K_{1} through the identity: U^{T}(\sum K_{i})U = T
call dcopy (n, t(2,1), 2, kwk(1,1,1), n+1)
call dcopy (n-1, t(1,2), 2, kwk(2,1,1), n+1)
for (j=1;j<n-1;j=j+1) call dset (n-j-1, 0.d0, kwk(j+2,j,1), 1)
for (i=2;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
for (j=1;j<=n;j=j+1)
call daxpy (n-j+1, -1.d0, kwk(j,j,i), 1, kwk(j,j,1), 1)
}
# fill the upper triangles of K_{i}
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
for (j=1;j<n;j=j+1) call dcopy (n-j, kwk(j+1,j,i), 1, kwk(j,j+1,i), n)
}
# decompose the tridiagonal matrix U^{T}DU
call dset (n, 10.d0 ** nlaht, work1, 1)
call daxpy (n, 1.d0, work1, 1, t(2,1), 2)
call dpbfa (t, 2, n, 1, info)
if ( info != 0 ) {
info = -2
return
}
# compute T^{-1}K_{i}
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
for (j=1;j<=n;j=j+1) call dpbsl (t, 2, n, 1, kwk(1,j,i))
}
#:::::::::: Compute the gradient and the Hessian ::::::::::
# compute -m x^{-T}T^{-m}K_{i}T^{-1}x: m = 2('v') or 1('m')
call dcopy (n, x, 1, work1, 1)
call dpbsl (t, 2, n, 1, work1)
if ( vmu != 'm' ) {
call dcopy (n, work1, 1, work2, 1)
call dscal (n, 2.d0, work2, 1)
}
else call dcopy (n, x, 1, work2, 1)
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
call dgemv ('t', n, n, 1.d0, kwk(1,1,i), n, work2, 1, 0.d0, work3, 1)
gwk1(i) = - ddot (n, work1, 1, work3, 1)
}
# compute - tr[T^{-m}K_{i}]: m = 2('v') or 1('m')
for (i=1;i<=nq;i=i+1) {
gwk2(i) = 0.d0
if ( theta(i) <= -25.d0 ) next
for (j=1;j<=n;j=j+1) {
if ( vmu != 'm' ) {
call dcopy (n, kwk(1,j,i), 1, work1, 1)
call dpbsl (t, 2, n, 1, work1)
gwk2(i) = gwk2(i) - work1(j)
}
else gwk2(i) = gwk2(i) - kwk(j,j,i)
}
}
if ( vmu != 'm' ) {
# compute 2 x^{T}T^{-1} [K_{i}T^{-2}K_{j}+T^{-1}K_{i}T^{-1}K_{j}
# +K_{i}T^{-1}K_{j}T^{-1}]T^{-1}x for 'v'
call dcopy (n, x, 1, work1, 1)
call dpbsl (t, 2, n, 1, work1)
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
call dgemv ('n', n, n, 1.d0, kwk(1,1,i), n, work1, 1, 0.d0, work2, 1)
for (j=1;j<=i;j=j+1) {
if ( theta(j) <= -25.d0 ) next
call dgemv ('n', n, n, 1.d0, kwk(1,1,j), n, work1, 1, 0.d0, work3, 1)
hwk1(i,j) = 2.d0 * ddot (n, work2, 1, work3, 1)
call dgemv ('t', n, n, 1.d0, kwk(1,1,j), n, work1, 1, 0.d0, work3, 1)
hwk1(i,j) = hwk1(i,j) + 2.d0 * ddot (n, work2, 1, work3, 1)
}
call dgemv ('t', n, n, 1.d0, kwk(1,1,i), n, work1, 1, 0.d0, work2, 1)
for (j=1;j<=i;j=j+1) {
if ( theta(j) <= -25.d0 ) next
call dgemv ('n', n, n, 1.d0, kwk(1,1,j), n, work1, 1, 0.d0, work3, 1)
hwk1(i,j) = hwk1(i,j) + 2.d0 * ddot (n, work2, 1, work3, 1)
}
}
}
else {
# compute 2 x^{T} T^{-1}K_{i}T^{-1}K_{j}T^{-1}x for 'm'
call dcopy (n, x, 1, work1, 1)
call dpbsl (t, 2, n, 1, work1)
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
call dgemv ('n', n, n, 1.d0, kwk(1,1,i), n, work1, 1, 0.d0, work2, 1)
for (j=1;j<=i;j=j+1) {
if ( theta(j) <= -25.d0 ) next
call dgemv ('t', n, n, 1.d0, kwk(1,1,j), n, x, 1, 0.d0, work3, 1)
hwk1(i,j) = 2.d0 * ddot (n, work2, 1, work3, 1)
}
}
}
# adjust diagonal
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
hwk1(i,i) = hwk1(i,i) + gwk1(i)
}
# compute m tr[T^{-m}K_{i}T^{-1}K_{j}]: m = 2('v') or 1('m')
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
for (m=1;m<=i;m=m+1) {
hwk2(i,m) = 0.d0
if ( theta(m) <= -25.d0 ) next
for (j=1;j<=n;j=j+1) {
if ( vmu != 'm' ) {
call dcopy (n, kwk(1,j,m), 1, work1, 1)
call dpbsl (t, 2, n, 1, work1)
hwk2(i,m) = hwk2(i,m) + 2.d0 * ddot (n, kwk(j,1,i), n, work1, 1)
}
else hwk2(i,m) = hwk2(i,m) + ddot (n, kwk(j,1,i), n, kwk(1,j,m), 1)
}
}
}
# adjust diagonal
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
hwk2(i,i) = hwk2(i,i) + gwk2(i)
}
# compute the gradient
if ( vmu == 'v' ) {
trc = dble (nobs) * 10.d0 ** (-nlaht) * varht / score
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
gra(i) = gwk1(i) / trc / trc - 2.d0 * score * gwk2(i) / trc / dble(nobs)
}
call dscal (nq, dble (nobs), gra, 1)
}
if ( vmu == 'u' ) {
dum = 10.d0 ** nlaht
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
gra(i) = dum * dum * gwk1(i) - 2.d0 * varht * dum * gwk2(i)
}
call dscal (nq, 1.d0/dble (n), gra, 1)
}
if ( vmu == 'm' ) {
det = 10.d0 ** (-nlaht) * varht / score
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
gra(i) = gwk1(i) / det - dble (nobs) / dble (n) * score * gwk2(i)
}
call dscal (nq, 1.d0 / dble (nobs), gra, 1)
}
# compute the Hessian
if ( vmu == 'v' ) {
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
for (j=1;j<=i;j=j+1) {
if ( theta(j) <= -25.d0 ) next
hes(i,j) = hwk1(i,j) / trc / trc - 2.d0 * gwk1(i) * gwk2(j) / trc ** 3_
- 2.d0 * gwk1(j) * gwk2(i) / trc ** 3 - 2.d0 * score * hwk2(i,j)_
/ trc / dble (nobs) + 6.d0 * score * gwk2(i) * gwk2(j)_
/ trc / trc / dble (nobs)
}
call dscal (i, dble (nobs), hes(i,1), ldh)
}
}
if ( vmu == 'u' ) {
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
for (j=1;j<=i;j=j+1) {
if ( theta(j) <= -25.d0 ) next
hes(i,j) = dum * dum * hwk1(i,j) - 2.d0 * varht * dum * hwk2(i,j)
}
call dscal (i, 1.d0/dble (n), hes(i,1), ldh)
}
}
if ( vmu == 'm' ) {
for (i=1;i<=nq;i=i+1) {
if ( theta(i) <= -25.d0 ) next
for (j=1;j<=i;j=j+1) {
if ( theta(j) <= -25.d0 ) next
hes(i,j) = hwk1(i,j) / det - gwk1(i) * gwk2(j) / det / dble (n)_
- gwk1(j) * gwk2(i) / det / dble (n) - dble (nobs)_
/ dble (n) * score * hwk2(i,j) + dble (nobs)_
/ dble (n) ** 2 * score * gwk2(i) * gwk2(j)
}
call dscal (i, 1.d0 / dble (nobs), hes(i,1), ldh)
}
}
return
end
#....................................................................................
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.