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