Revision 2b3a85a87089d85cfbcc67cac5d6778f94719b8d authored by Roger Koenker on 15 July 2019, 13:00:03 UTC, committed by cran-robot on 15 July 2019, 13:00:03 UTC
1 parent 5071707
brute.r
``````# New Version of Brute Force Algorithm for the Powell Estimator
#
# Roger Koenker:  last revision: 12 February 2008
#
# Problem:  || min{Ax,c} - b ||_tau  = min!
#
# The matrix H is p by m matrix of the n choose p basis indices h on input.
# When called it is assumed (!!) that U = A[H[,1],]^{-1} and xh = U b[H[,1]].
subroutine brutpow(n,p,m,H,A,b,c,x,tau,U,xh,d,jminz,nflag)

integer n,p,m
double precision x(p),A(n,p),b(n),c(n)
double precision U(p,p),d(p),xh(p)
double precision zero, one,tau,pow,minz,z
integer H(p,m),k,findk,jminz,nflag

PARAMETER(zero = 0.0d0, one = 1.d0)

jminz = 1
minz  = pow(n,p,x,A,b,c,tau)
do j = 2,m {
k = findk(p,H(1,j),H(1,j-1))
if(k == 0)
{nflag = 4; return}
call pivot(n,p,H(1,j-1),H(k,j),H(k,j-1),A,U,d,xh,nflag)
if(nflag > 0)
return
do i = 1,p{
xh(i) = b(H(i,j))
}
call dgemv('N',p,p,one,U,p,xh,1,zero,x,1)
z = pow(n,p,x,A,b,c,tau)
if(z < minz) {
minz  = z
jminz = j
}
}
return
end

#####################################################

integer function findk(p,h,g)
integer p,k,h(p),g(p)
findk = 0
do k = 1,p{
if(h(k) != g(k))
{findk = k; break}
}
return
end

``````

Computing file changes ...