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