https://github.com/cran/quantreg
Tip revision: de888d3d89a7022314982fb84db6825413d28dc6 authored by Roger Koenker on 02 October 2020, 03:20:02 UTC
version 5.73
version 5.73
Tip revision: de888d3
penalty.f
C Output from Public domain Ratfor, version 1.0
subroutine penalty(n,m,q,x,y,bnd,tlist,tlptr,tlend,rax,jax,ned,eps
*,ierr)
integer n,m,q,lp,lpl,ned,ierr
integer bnd(n),tlist(q),tlptr(q),tlend(n),n4(4),p4(4),jax(m)
double precision x(n),y(n),rax(m),eps
double precision x4(4),y4(4),g4(4)
logical orient
ned = 0
do23000 i=1,n
lpl = tlend(i)
lp = lpl
23002 continue
lp = tlptr(lp)
j = iabs(tlist(lp))
if(j .gt. i)then
n4(1) = i
n4(2) = j
call fadjs(n4,n,q,tlist,tlptr,tlend)
if(bnd(i)*bnd(j) .eq. 0)then
ned = ned + 1
do23009 k = 1,4
x4(k) = x(n4(k))
y4(k) = y(n4(k))
23009 continue
23010 continue
if(orient(x4,y4))then
call iswap(1,n4(3),1,n4(4),1)
call dswap(1,x4(3),1,x4(4),1)
call dswap(1,y4(3),1,y4(4),1)
endif
call ggap(x4,y4,g4,eps,ierr)
if(ierr .eq. 1)then
return
endif
call srtpai(n4,1,p4,1,4)
do23015 k = 1,4
rax((ned - 1)*4 + k) = g4(p4(k))
jax((ned - 1)*4 + k) = n4(p4(k))
23015 continue
23016 continue
if(ned*4 .gt. m)then
return
endif
endif
endif
if(lp .eq. lpl)then
goto 23004
endif
23003 goto 23002
23004 continue
23000 continue
23001 continue
return
end
logical function orient(x,y)
double precision x(4), y(4)
orient = (y(2) -y(1))*(x(3)-x(4))+(x(1)-x(2))*(y(3)-y(4)) .gt. 0
return
end
subroutine fadjs(n4,n,q,tlist,tlptr,tlend)
integer n,q,vp,vpl,v,v0,match
integer n4(4),tlist(q),tlptr(q),tlend(n)
match = 0
vpl = tlend(n4(1))
vp = vpl
k = 0
23021 continue
k = k+1
vp = tlptr(vp)
v = tlist(vp)
if(k.gt.1 .and. iabs(v) .eq. n4(2))then
n4(3) = iabs(v0)
match = 1
goto 23022
endif
if(match .gt. 0)then
n4(4) = iabs(v)
goto 23023
endif
v0 = v
23022 goto 23021
23023 continue
return
end
subroutine ggap(x,y,g,eps,ierr)
double precision x(4),y(4),g(4),w(2,4),h(2),d1,d2,eps
d1 = -x(2) * y(1) + x(3) * y(1) + x(1) * y(2) - x(3) * y(2) - x(1)
* * y(3) + x(2) * y(3)
d2 = -x(2) * y(1) + x(4) * y(1) + x(1) * y(2) - x(4) * y(2) - x(1)
* * y(4) + x(2) * y(4)
if(dabs(d1) .lt. eps .or. dabs(d2) .lt. eps)then
ierr = 1
return
endif
h(1) = -(y(1) - y(2))
h(2) = (x(1) - x(2))
w(1, 1) = (y(2) - y(3))/d1 - (y(2) - y(4))/d2
w(2, 1) = (x(3) - x(2))/d1 - (x(4) - x(2))/d2
w(1, 2) = (y(3) - y(1))/d1 - (y(4) - y(1))/d2
w(2, 2) = (x(1) - x(3))/d1 - (x(1) - x(4))/d2
w(1, 3) = (y(1) - y(2))/d1
w(2, 3) = (x(2) - x(1))/d1
w(1, 4) = (y(2) - y(1))/d2
w(2, 4) = (x(1) - x(2))/d2
do23030 i = 1,4
g(i) = h(1)*w(1,i)+h(2)*w(2,i)
23030 continue
23031 continue
ierr = 0
return
end