1
2#:::::::::
3#   reg
4#:::::::::
5
6subroutine  reg (sr, nobs, nnull, q, nxi, y, method, alpha, varht,
7                 score, dc, mchpr, v, mu, jpvt, wk, rkv, info)
8
9double precision  sr(nobs,*), q(nxi,*), y(*), alpha, varht, score, dc(*),
10                  mchpr, v(nnull+nxi,*), mu(*), wk(*)
11integer  nobs, nnull, nxi, method, jpvt(*), rkv, info
12
13double precision  ddot, dasum, rss, trc, dum
14integer  i, j, nn, idamax, infowk, idum
15
16info = 0
17nn = nnull + nxi
18
19#   form linear system
20for (i=1;i<=nn;i=i+1) {
21    mu(i) = ddot (nobs, sr(1,i), 1, y, 1)
22    for (j=i;j<=nn;j=j+1) {
23        v(i,j) = ddot (nobs, sr(1,i), 1, sr(1,j), 1)
24        if (i>nnull)  v(i,j) = v(i,j) + q(i-nnull,j-nnull)
25    }
26}
27#   Cholesky decomposition
28infowk = 0
29for (i=1;i<=nn;i=i+1)  infowk = infowk + jpvt(i)
30call  dchdc (v, nn, nn, wk, jpvt, 1, rkv)
31j = idamax (rkv-infowk, v(infowk+1,infowk+1), nn+1)
32while (v(rkv,rkv)<v(infowk+j,infowk+j)*dsqrt(mchpr))  rkv = rkv - 1
33for (i=rkv+1;i<=nn;i=i+1) {
34    v(i,i) = v(j,j)
35    call  dset (i-rkv-1, 0.d0, v(rkv+1,i), 1)
36}
37#   obtain coefficients
38call  dcopy (nn, mu, 1, dc, 1)
39call  dprmut (dc, nn, jpvt, 0)
40call  dtrsl (v, nn, nn, dc, 11, infowk)
41call  dset (nn-rkv, 0.d0, dc(rkv+1), 1)
42call  dtrsl (v, nn, nn, dc, 01, infowk)
43call  dprmut (dc, nn, jpvt, 1)
44#   early return
45if (method==4)  return
46#   calculate residuals
47for (i=1;i<=nobs;i=i+1)  wk(i) = y(i) - ddot (nn, sr(i,1), nobs, dc, 1)
48#   diagonal of smoothing matrix
49if (method==5) {
50    wk(nobs+1) = ddot (nobs, wk, 1, wk, 1) / dfloat (nobs)
51    for (i=1;i<=nobs;i=i+1) {
52        call  dcopy (nn, sr(i,1), nobs, mu, 1)
53        call  dprmut (mu, nn, jpvt, 0)
54        call  dtrsl (v, nn, nn, mu, 11, infowk)
55        wk(i) = ddot (nn, mu, 1, mu, 1)
56    }
57    return
58}
59if (method==3) {
60    #   GML
61    rss = ddot (nobs, y, 1, wk, 1)
62    #   determinant
63    if (nnull>0) {
64        call  dqrdc (sr, nobs, nobs, nnull, wk, idum, dum, 0)
65        for (i=1;i<=nxi;i=i+1) {
66            call  dqrsl (sr, nobs, nobs, nnull, wk, sr(1,nnull+i),
67                         dum, sr(1,nnull+i), dum, dum, dum, 01000, infowk)
68        }
69    }
70    call  dcopy (nxi, q, nxi+1, wk, 1)
71    for (i=1;i<=nxi;i=i+1) {
72        for (j=i;j<=nxi;j=j+1)
73            q(i,j) = q(i,j) + ddot (nobs-nnull, sr(nnull+1,nnull+i), 1,
74                                    sr(nnull+1,nnull+j), 1)
75    }
76    for (i=1;i<=nxi;i=i+1) {
77        for (j=i;j<=nxi;j=j+1) {
78            sr(i,j) = q(i,j)
79            sr(j,i) = q(i,j)
80            q(i,j) = q(j,i)
81        }
82    }
83    call  dcopy (nxi, wk, 1, q, nxi+1)
84#    call  rs (nobs, nxi, sr, mu, 0, dum, wk, y, info)
85    call  dsyev ('n', 'u', nxi, sr, nobs, mu, wk, 3*nxi, info)
86    trc = 0.d0
87    for (i=1;i<=rkv-nnull;i=i+1)  trc = trc + dlog (mu(nxi-i+1))
88#    call  rs (nxi, nxi, q, mu, 0, dum, wk, y, info)
89    call  dsyev ('n', 'u', nxi, q, nxi, mu, wk, 3*nxi, info)
90    for (i=1;i<=rkv-nnull;i=i+1)  trc = trc - dlog (mu(nxi-i+1))
91    #   return values
92    score = rss / dfloat (nobs) * dexp (trc/dfloat(nobs-nnull))
93    varht = rss / dfloat (nobs-nnull)
94}
95else {
96    #   GCV or Cp
97    rss = ddot (nobs, wk, 1, wk, 1) / dfloat (nobs)
98    #   trace
99    for (i=1;i<=nobs;i=i+1) {
100        call  dcopy (nn, sr(i,1), nobs, mu, 1)
101        call  dprmut (mu, nn, jpvt, 0)
102        call  dtrsl (v, nn, nn, mu, 11, infowk)
103        wk(i) = ddot (nn, mu, 1, mu, 1)
104    }
105    trc = dasum (nobs, wk, 1) / dfloat (nobs)
106    #   return values
107    if (method==2) {
108        score = rss / (1.d0-alpha*trc)**2
109        varht = rss / (1.d0-trc)
110    }
111    else  score = rss + 2.d0 * varht * alpha * trc
112}
113wk(1) = rss
114wk(2) = trc
115
116return
117end
118
119
120#::::::::::::
121#   regaux
122#::::::::::::
123
124subroutine  regaux (v, nn, jpvt, rkv, r, nr, sms, nnull, wk)
125
126double precision  v(nn,*), r(nn,*), sms(nnull,*), wk(nn,*)
127integer  nn, jpvt(*), rkv, nr, nnull
128
129double precision  ddot
130integer  i, j, infowk
131
132#   drcr
133for (i=1;i<=nr;i=i+1) {
134    call  dprmut (r(1,i), nn, jpvt, 0)
135    call  dtrsl (v, nn, nn, r(1,i), 11, infowk)
136    if (nn-rkv>0)  call  dset (nn-rkv, 0.d0, r(rkv+1,i), 1)
137    call  dtrsl (v, nn, nn, r(1,i), 01, infowk)
138    call  dprmut (r(1,i), nn, jpvt, 1)
139}
140#   sms
141call  dset (nn*nnull, 0.d0, wk, 1)
142call  dset (nnull, 1.d0, wk, nn+1)
143for (i=1;i<=nnull;i=i+1)
144    call  dtrsl (v, nn, nn, wk(1,i), 11, infowk)
145for (i=1;i<=nnull;i=i+1) {
146    for (j=i;j<=nnull;j=j+1) {
147        sms(i,j) = ddot (nn, wk(1,i), 1, wk(1,j), 1)
148        sms(j,i) = sms(i,j)
149    }
150}
151
152return
153end
154