1C Output from Public domain Ratfor, version 1.01 2 subroutine reg (sr, nobs, nnull, q, nxi, y, method, alpha, varht, 3 *score, dc, mchpr, v, mu, jpvt, wk, rkv, info) 4 double precision sr(nobs,*), q(nxi,*), y(*), alpha, varht, score, 5 *dc(*), mchpr, v(nnull+nxi,*), mu(*), wk(*) 6 integer nobs, nnull, nxi, method, jpvt(*), rkv, info 7 double precision ddot, dasum, rss, trc, dum 8 integer i, j, nn, idamax, infowk, idum 9 info = 0 10 nn = nnull + nxi 11 i=1 1223000 if(.not.(i.le.nn))goto 23002 13 mu(i) = ddot (nobs, sr(1,i), 1, y, 1) 14 j=i 1523003 if(.not.(j.le.nn))goto 23005 16 v(i,j) = ddot (nobs, sr(1,i), 1, sr(1,j), 1) 17 if(i.gt.nnull)then 18 v(i,j) = v(i,j) + q(i-nnull,j-nnull) 19 endif 2023004 j=j+1 21 goto 23003 2223005 continue 2323001 i=i+1 24 goto 23000 2523002 continue 26 infowk = 0 27 i=1 2823008 if(.not.(i.le.nn))goto 23010 29 infowk = infowk + jpvt(i) 3023009 i=i+1 31 goto 23008 3223010 continue 33 call dchdc (v, nn, nn, wk, jpvt, 1, rkv) 34 j = idamax (rkv-infowk, v(infowk+1,infowk+1), nn+1) 3523011 if(v(rkv,rkv).lt.v(infowk+j,infowk+j)*dsqrt(mchpr))then 36 rkv = rkv - 1 37 goto 23011 38 endif 3923012 continue 40 i=rkv+1 4123013 if(.not.(i.le.nn))goto 23015 42 v(i,i) = v(j,j) 43 call dset (i-rkv-1, 0.d0, v(rkv+1,i), 1) 4423014 i=i+1 45 goto 23013 4623015 continue 47 call dcopy (nn, mu, 1, dc, 1) 48 call dprmut (dc, nn, jpvt, 0) 49 call dtrsl (v, nn, nn, dc, 11, infowk) 50 call dset (nn-rkv, 0.d0, dc(rkv+1), 1) 51 call dtrsl (v, nn, nn, dc, 01, infowk) 52 call dprmut (dc, nn, jpvt, 1) 53 if(method.eq.4)then 54 return 55 endif 56 i=1 5723018 if(.not.(i.le.nobs))goto 23020 58 wk(i) = y(i) - ddot (nn, sr(i,1), nobs, dc, 1) 5923019 i=i+1 60 goto 23018 6123020 continue 62 if(method.eq.5)then 63 wk(nobs+1) = ddot (nobs, wk, 1, wk, 1) / dfloat (nobs) 64 i=1 6523023 if(.not.(i.le.nobs))goto 23025 66 call dcopy (nn, sr(i,1), nobs, mu, 1) 67 call dprmut (mu, nn, jpvt, 0) 68 call dtrsl (v, nn, nn, mu, 11, infowk) 69 wk(i) = ddot (nn, mu, 1, mu, 1) 7023024 i=i+1 71 goto 23023 7223025 continue 73 return 74 endif 75 if(method.eq.3)then 76 rss = ddot (nobs, y, 1, wk, 1) 77 if(nnull.gt.0)then 78 call dqrdc (sr, nobs, nobs, nnull, wk, idum, dum, 0) 79 i=1 8023030 if(.not.(i.le.nxi))goto 23032 81 call dqrsl (sr, nobs, nobs, nnull, wk, sr(1,nnull+i), dum, sr(1,nn 82 *ull+i), dum, dum, dum, 01000, infowk) 8323031 i=i+1 84 goto 23030 8523032 continue 86 endif 87 call dcopy (nxi, q, nxi+1, wk, 1) 88 i=1 8923033 if(.not.(i.le.nxi))goto 23035 90 j=i 9123036 if(.not.(j.le.nxi))goto 23038 92 q(i,j) = q(i,j) + ddot (nobs-nnull, sr(nnull+1,nnull+i), 1, sr(nnu 93 *ll+1,nnull+j), 1) 9423037 j=j+1 95 goto 23036 9623038 continue 9723034 i=i+1 98 goto 23033 9923035 continue 100 i=1 10123039 if(.not.(i.le.nxi))goto 23041 102 j=i 10323042 if(.not.(j.le.nxi))goto 23044 104 sr(i,j) = q(i,j) 105 sr(j,i) = q(i,j) 106 q(i,j) = q(j,i) 10723043 j=j+1 108 goto 23042 10923044 continue 11023040 i=i+1 111 goto 23039 11223041 continue 113 call dcopy (nxi, wk, 1, q, nxi+1) 114 call dsyev ('n', 'u', nxi, sr, nobs, mu, wk, 3*nxi, info) 115 trc = 0.d0 116 i=1 11723045 if(.not.(i.le.rkv-nnull))goto 23047 118 trc = trc + dlog (mu(nxi-i+1)) 11923046 i=i+1 120 goto 23045 12123047 continue 122 call dsyev ('n', 'u', nxi, q, nxi, mu, wk, 3*nxi, info) 123 i=1 12423048 if(.not.(i.le.rkv-nnull))goto 23050 125 trc = trc - dlog (mu(nxi-i+1)) 12623049 i=i+1 127 goto 23048 12823050 continue 129 score = rss / dfloat (nobs) * dexp (trc/dfloat(nobs-nnull)) 130 varht = rss / dfloat (nobs-nnull) 131 else 132 rss = ddot (nobs, wk, 1, wk, 1) / dfloat (nobs) 133 i=1 13423051 if(.not.(i.le.nobs))goto 23053 135 call dcopy (nn, sr(i,1), nobs, mu, 1) 136 call dprmut (mu, nn, jpvt, 0) 137 call dtrsl (v, nn, nn, mu, 11, infowk) 138 wk(i) = ddot (nn, mu, 1, mu, 1) 13923052 i=i+1 140 goto 23051 14123053 continue 142 trc = dasum (nobs, wk, 1) / dfloat (nobs) 143 if(method.eq.2)then 144 score = rss / (1.d0-alpha*trc)**2 145 varht = rss / (1.d0-trc) 146 else 147 score = rss + 2.d0 * varht * alpha * trc 148 endif 149 endif 150 wk(1) = rss 151 wk(2) = trc 152 return 153 end 154 subroutine regaux (v, nn, jpvt, rkv, r, nr, sms, nnull, wk) 155 double precision v(nn,*), r(nn,*), sms(nnull,*), wk(nn,*) 156 integer nn, jpvt(*), rkv, nr, nnull 157 double precision ddot 158 integer i, j, infowk 159 i=1 16023056 if(.not.(i.le.nr))goto 23058 161 call dprmut (r(1,i), nn, jpvt, 0) 162 call dtrsl (v, nn, nn, r(1,i), 11, infowk) 163 if(nn-rkv.gt.0)then 164 call dset (nn-rkv, 0.d0, r(rkv+1,i), 1) 165 endif 166 call dtrsl (v, nn, nn, r(1,i), 01, infowk) 167 call dprmut (r(1,i), nn, jpvt, 1) 16823057 i=i+1 169 goto 23056 17023058 continue 171 call dset (nn*nnull, 0.d0, wk, 1) 172 call dset (nnull, 1.d0, wk, nn+1) 173 i=1 17423061 if(.not.(i.le.nnull))goto 23063 175 call dtrsl (v, nn, nn, wk(1,i), 11, infowk) 17623062 i=i+1 177 goto 23061 17823063 continue 179 i=1 18023064 if(.not.(i.le.nnull))goto 23066 181 j=i 18223067 if(.not.(j.le.nnull))goto 23069 183 sms(i,j) = ddot (nn, wk(1,i), 1, wk(1,j), 1) 184 sms(j,i) = sms(i,j) 18523068 j=j+1 186 goto 23067 18723069 continue 18823065 i=i+1 189 goto 23064 19023066 continue 191 return 192 end 193