# 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