https://github.com/cran/quantreg
Tip revision: 65f51b515659b5ccec38dd61f7ccda712180d272 authored by Roger Koenker on 08 August 1977, 00:00:00 UTC
version 3.05
version 3.05
Tip revision: 65f51b5
rqbr.f
subroutine rqbr(m,nn,m5,n3,n4,a,b,t,toler,ift,x,e,s,wa,wb,nsol,
&ndsol,sol,dsol,lsol,h,qn,cutoff,ci,tnmat,big,lci1)
integer i,j,k,kl,kount,kr,l,lsol,m,m1,m2,m3,m4,m5,ift
integer n,n1,n2,n3,nsol,ndsol,out,s(m),h(nn,nsol)
integer nn,n4,idxcf
logical stage,test,init,iend,lup
logical lci1,lci2,skip
double precision a1,aux,b1,big,d,dif,pivot,smax,t,t0,t1,tnt
double precision min,max,toler,zero,half,one,two
double precision b(m),sol(n3,nsol),a(m,nn),x(nn),wa(m5,n4),wb(m)
double precision sum,e(m),dsol(m,ndsol)
double precision qn(nn),cutoff,ci(4,nn),tnmat(4,nn),tnew,told,tn
data zero/0.d0/
data one/1.d0/
data two/2.d0/
n=nn
ift = 0
wa(m+2,nn+1) = one
if(.not.(m5.ne.m+5))goto 23000
ift = 3
23000 continue
if(.not.(n3.ne.n+3))goto 23002
ift = 4
23002 continue
if(.not.(n4.ne.n+4))goto 23004
ift = 5
23004 continue
if(.not.(m.le.zero.or.n.le.zero))goto 23006
ift = 6
23006 continue
if(.not.(ift.le.two))goto 23008
half = one/two
iend = .true.
lci2 = .false.
lup = .true.
skip = .false.
idxcf = 0
tnew = zero
tn = zero
m1 = m+1
n1 = n+1
n2 = n+2
m2 = m+2
m3 = m+3
m4 = m+4
do 23010 j = 1,n
x(j) = zero
23010 continue
do 23012 i = 1,m
e(i) = zero
23012 continue
if(.not.(t.lt.zero.or.t.gt.one))goto 23014
t0 = one/float(m)-toler
t1 = one-t0
t = t0
iend = .false.
lci1 = .false.
23014 continue
23016 continue
do 23019 i = 1,m
k = 1
do 23021 j = 1,nn
if(.not.(k.le.nn))goto 23023
if(.not.(j.eq.idxcf))goto 23025
skip = .true.
goto 23026
23025 continue
skip = .false.
23026 continue
if(.not.(.not.skip))goto 23027
wa(i,k) = a(i,j)
k = k+1
23027 continue
23023 continue
23021 continue
wa(i,n4) = n+i
wa(i,n2) = b(i)
if(.not.(idxcf .ne. 0))goto 23029
wa(i,n3) = tnew*a(i,idxcf)
goto 23030
23029 continue
wa(i,n3) = zero
23030 continue
wa(i,n1) = wa(i,n2)-wa(i,n3)
if(.not.(wa(i,n1).lt.zero))goto 23031
do 23033 j = 1,n4
wa(i,j) = -wa(i,j)
23033 continue
23031 continue
23019 continue
do 23035 j = 1,n
wa(m4,j) = j
wa(m2,j) = zero
wa(m3,j) = zero
do 23037 i = 1,m
aux = sign(one,wa(m4,j))*wa(i,j)
wa(m2,j) = wa(m2,j)+aux*(one-sign(one,wa(i,n4)))
wa(m3,j) = wa(m3,j)+aux*sign(one,wa(i,n4))
23037 continue
wa(m3,j) = two*wa(m3,j)
23035 continue
dif = zero
init = .false.
if(.not.(.not.lci2))goto 23039
do 23041 k = 1,n
wa(m5,k) = zero
do 23043 i = 1,m
wa(m5,k) = wa(m5,k)+a(i,k)
23043 continue
wa(m5,k) = wa(m5,k)/float(m)
23041 continue
23039 continue
lsol = 1
kount = 0
23045 continue
do 23048 j = 1,n
wa(m1,j) = wa(m2,j)+wa(m3,j)*t
23048 continue
if(.not.(.not.init))goto 23050
stage = .true.
kr = 1
kl = 1
go to 30
23050 continue
23052 continue
stage = .false.
23055 continue
max = -big
do 23058 j = kr,n
d = wa(m1,j)
if(.not.(d.lt.zero))goto 23060
if(.not.(d.gt.(-two)))goto 23062
goto 23058
23062 continue
d = -d-two
23060 continue
if(.not.(d.gt.max))goto 23064
max = d
in = j
23064 continue
23058 continue
if(.not.(max.le.toler))goto 23066
goto 23054
23066 continue
if(.not.(wa(m1,in).le.zero))goto 23068
do 23070 i = 1,m4
wa(i,in) = -wa(i,in)
23070 continue
wa(m1,in) = wa(m1,in)-two
wa(m2,in) = wa(m2,in)-two
23068 continue
23072 continue
k = 0
do 23075 i = kl,m
d = wa(i,in)
if(.not.(d.gt.toler))goto 23077
k = k+1
wb(k) = wa(i,n1)/d
s(k) = i
test = .true.
23077 continue
23075 continue
23079 continue
if(.not.(k.le.0))goto 23082
test = .false.
goto 23083
23082 continue
min = big
do 23084 i = 1,k
if(.not.(wb(i).lt.min))goto 23086
j = i
min = wb(i)
out = s(i)
23086 continue
23084 continue
wb(j) = wb(k)
s(j) = s(k)
k = k-1
23083 continue
if(.not.(.not.test.and.stage))goto 23088
goto 23081
23088 continue
if(.not.(.not.test))goto 23090
goto 23047
23090 continue
pivot = wa(out,in)
if(.not.(wa(m1,in)-pivot-pivot.le.toler))goto 23092
go to 10
23092 continue
do 23094 j = kr,n3
d = wa(out,j)
wa(m1,j) = wa(m1,j)-d-d
wa(m2,j) = wa(m2,j)-d-d
wa(out,j) = -d
23094 continue
wa(out,n4) = -wa(out,n4)
23080 goto 23079
23081 continue
do 23096 i = 1,m4
d = wa(i,kr)
wa(i,kr) = wa(i,in)
wa(i,in) = d
23096 continue
kr = kr+1
go to 20
10 do 23098 j = kr,n3
if(.not.(j.ne.in))goto 23100
wa(out,j) = wa(out,j)/pivot
23100 continue
23098 continue
do 23102 i = 1,m3
if(.not.(i.ne.out))goto 23104
d = wa(i,in)
do 23106 j = kr,n3
if(.not.(j.ne.in))goto 23108
wa(i,j) = wa(i,j)-d*wa(out,j)
23108 continue
23106 continue
23104 continue
23102 continue
do 23110 i = 1,m3
if(.not.(i.ne.out))goto 23112
wa(i,in) = -wa(i,in)/pivot
23112 continue
23110 continue
wa(out,in) = one/pivot
d = wa(out,n4)
wa(out,n4) = wa(m4,in)
wa(m4,in) = d
kount = kount+1
if(.not.(.not.stage))goto 23114
goto 23074
23114 continue
kl = kl+1
do 23116 j = kr,n4
d = wa(out,j)
wa(out,j) = wa(kount,j)
wa(kount,j) = d
23116 continue
20 if(.not.(kount+kr.eq.n1))goto 23118
goto 23057
23118 continue
30 max = -one
do 23120 j = kr,n
if(.not.(abs(wa(m4,j)).le.n))goto 23122
d = abs(wa(m1,j))
if(.not.(d.gt.max))goto 23124
max = d
in = j
23124 continue
23122 continue
23120 continue
if(.not.(wa(m1,in).lt.zero))goto 23126
do 23128 i = 1,m4
wa(i,in) = -wa(i,in)
23128 continue
23126 continue
23073 goto 23072
23074 continue
23056 goto 23055
23057 continue
23053 goto 23052
23054 continue
if(.not.(kr.eq.1))goto 23130
do 23132 j = 1,n
d = abs(wa(m1,j))
if(.not.(d.le.toler.or.two-d.le.toler))goto 23134
ift = 1
wa(m2,nn+1) = zero
go to 80
23134 continue
23132 continue
23130 continue
80 kount = 0
sum = zero
if(.not.(.not.lci2))goto 23136
do 23138 i = 1,kl-1
k = wa(i,n4)*sign(one,wa(i,n4))
x(k) = wa(i,n1)*sign(one,wa(i,n4))
23138 continue
23136 continue
do 23140 i = 1,n
kd = abs(wa(m4,i))-n
dsol(kd,lsol) = one+wa(m1,i)/two
if(.not.(wa(m4,i).lt.zero))goto 23142
dsol(kd,lsol) = one-dsol(kd,lsol)
23142 continue
if(.not.(.not.lci2))goto 23144
sum = sum + x(i)*wa(m5,i)
sol(i+3,lsol) = x(i)
h(i,lsol) = kd
23144 continue
23140 continue
do 23146 i = kl,m
kd = abs(wa(i,n4))-n
if(.not.(wa(i,n4).lt.zero))goto 23148
dsol(kd,lsol) = zero
23148 continue
if(.not.(wa(i,n4).gt.zero))goto 23150
dsol(kd,lsol) = one
23150 continue
23146 continue
if(.not.(.not.lci2))goto 23152
sol(1,lsol) = smax
sol(2,lsol) = sum
sum = zero
do 23154 j=kl,m
d = wa(j,n1)*sign(one,wa(j,n4))
sum = sum + d*(smax + half*(sign(one,d) - one))
23154 continue
sol(3,lsol) = sum
do 23156 i=1,m
dsol(i,lsol+1) = dsol(i,lsol)
23156 continue
23152 continue
if(.not.(lci2))goto 23158
a1 = zero
do 23160 i = 1,m
a1 = a1+a(i,idxcf)*(dsol(i,lsol)+t-one)
23160 continue
tn = a1/sqrt(qn(idxcf)*t*(one-t))
if(.not.(abs(tn).lt.cutoff))goto 23162
if(.not.(lup))goto 23164
smax = big
goto 23165
23164 continue
smax = -big
23165 continue
do 23166 i =1,kl-1
k = wa(i,n4)*sign(one,wa(i,n4))
sol(k,1) = wa(i,n2)*sign(one,wa(i,n4))
sol(k,2) = wa(i,n3)*sign(one,wa(i,n4))/tnew
23166 continue
do 23168 i = kl,m
a1 = zero
b1 = zero
k = wa(i,n4)*sign(one,wa(i,n4))-n
l = 1
do 23170 j = 1,n
if(.not.(j.eq.idxcf))goto 23172
l = l+1
23172 continue
a1 = a1 + a(k,l)*sol(j,1)
b1 = b1 + a(k,l)*sol(j,2)
l = l+1
23170 continue
tnt = (b(k)-a1)/(a(k,idxcf)-b1)
if(.not.(lup))goto 23174
if(.not.(tnt.gt.tnew))goto 23176
if(.not.(tnt.lt.smax))goto 23178
smax = tnt
out = i
23178 continue
23176 continue
goto 23175
23174 continue
if(.not.(tnt.lt.tnew))goto 23180
if(.not.(tnt.gt.smax))goto 23182
smax = tnt
out = i
23182 continue
23180 continue
23175 continue
23168 continue
if(.not.(lup))goto 23184
told = tnew
tnew = smax + toler
ci(3,idxcf) = told - toler
tnmat(3,idxcf) = tn
if(.not.(.not.(tnew .lt. big-toler)))goto 23186
ci(3,idxcf) = big
ci(4,idxcf) = big
tnmat(3,idxcf) = tn
tnmat(4,idxcf) = tn
lup = .false.
go to 70
23186 continue
goto 23185
23184 continue
told = tnew
tnew = smax - toler
ci(2,idxcf) = told + toler
tnmat(2,idxcf) = tn
if(.not.(.not.(tnew .gt. -big+toler)))goto 23188
ci(2,idxcf) = -big
ci(1,idxcf) = -big
tnmat(2,idxcf) = tn
tnmat(1,idxcf) = tn
lup = .true.
go to 60
23188 continue
23185 continue
do 23190 i = 1,m
wa(i,n3) = wa(i,n3)/told*tnew
wa(i,n1) = wa(i,n2) - wa(i,n3)
23190 continue
do 23192 j = kr,n3
d = wa(out,j)
wa(m1,j) = wa(m1,j) -d -d
wa(m2,j) = wa(m2,j) -d -d
wa(out,j) = -d
23192 continue
wa(out,n4) = -wa(out,n4)
init = .true.
goto 23163
23162 continue
if(.not.(lup))goto 23194
ci(4,idxcf) = tnew - toler
tnmat(4,idxcf) = tn
lup = .false.
go to 70
23194 continue
ci(1,idxcf) = tnew + toler
tnmat(1,idxcf) = tn
lup = .true.
go to 60
23195 continue
23163 continue
23158 continue
if(.not.((iend).and.(.not.lci2)))goto 23196
go to 40
23196 continue
if(.not.(.not.lci2))goto 23198
init = .true.
lsol = lsol+1
do 23200 i = 1,m
s(i) = zero
23200 continue
do 23202 j = 1,n
x(j) = zero
23202 continue
smax = two
do 23204 j = 1,n
b1 = wa(m3,j)
a1 = (-two-wa(m2,j))/b1
b1 = -wa(m2,j)/b1
if(.not.(a1.ge.t))goto 23206
if(.not.(a1.lt.smax))goto 23208
smax = a1
dif = (b1-a1)/two
23208 continue
23206 continue
if(.not.(b1.gt.t))goto 23210
if(.not.(b1.lt.smax))goto 23212
smax = b1
dif = (b1-a1)/two
23212 continue
23210 continue
23204 continue
tnt = smax+toler*(one+abs(dif))
if(.not.(tnt.ge.t1+toler))goto 23214
iend = .true.
23214 continue
t = tnt
if(.not.(iend))goto 23216
t = t1
23216 continue
23198 continue
23046 goto 23045
23047 continue
wa(m2,nn+1) = two
ift = 2
go to 50
40 if(.not.(lsol.gt.2))goto 23218
sol(1,1) = zero
sol(2,1) = zero
sol(3,1) = zero
sol(1,lsol) = one
sol(2,lsol) = zero
sol(3,lsol) = zero
do 23220 i = 1,m
dsol(i,1) = one
dsol(i,lsol) = zero
dsol(i,lsol+1) = zero
23220 continue
23218 continue
l = kl-1
do 23222 i = 1,l
if(.not.(wa(i,n1).lt.zero))goto 23224
do 23226 j = kr,n4
wa(i,j) = -wa(i,j)
23226 continue
23224 continue
23222 continue
50 sum = zero
if(.not.(.not.lci2))goto 23228
do 23230 i = kl,m
k = wa(i,n4)*sign(one,wa(i,n4))
d = wa(i,n1)*sign(one,wa(i,n4))
sum = sum+d*sign(one,d)*(half+sign(one,d)*(t-half))
k = k-n
e(k) = d
23230 continue
wa(m2,n2) = kount
wa(m1,n2) = n1-kr
wa(m1,n1) = sum
23228 continue
if(.not.(wa(m2,nn+1).eq.two))goto 23232
goto 23018
23232 continue
if(.not.(.not.lci1))goto 23234
goto 23018
23234 continue
if(.not.(.not.lci2))goto 23236
lci2 = .true.
n = nn-1
n1 = n+1
n2 = n+2
n3 = n+3
n4 = n+4
60 idxcf = idxcf+1
if(.not.(idxcf.gt.nn))goto 23238
goto 23018
23238 continue
70 if(.not.(lup))goto 23240
tnew = x(idxcf)+toler
told = tnew
ci(3,idxcf) = x(idxcf)
tnmat(3,idxcf) = zero
goto 23241
23240 continue
tnew = x(idxcf)-toler
told = tnew
ci(2,idxcf) = x(idxcf)
tnmat(2,idxcf) = zero
23241 continue
23236 continue
23017 goto 23016
23018 continue
do 23242 i=1,m
dsol(i,lsol) = dsol(i,lsol+1)
23242 continue
23008 continue
return
end