FUNCTION HR_LDP,G,h,phi
;+z=LDP(G,h,phi) implements Lawson and Hanson
;algorithm 23.27 to minimize ||z|| subject
;to Gz>=h. If the inequality constraints
;are consistent then phi=TRUE and a value will
;be returned for vector z. The value of
;phi=total(abs(z)) is returned. A very small
;value of phi indicates a possible inconsistent
;constraint set or a failure of the algorithm to
;converge to a good answer.
;
;H. Rhody
;May, 2004
;-

sg=size(G,/dim)
m=sg[1]
n=sg[0]
sh=size(h,/dim)
IF sh[1] NE m then Message,'G and h have incompatible dimensions'
E=[[transpose(G)],[transpose(h)]]
f=transpose([replicate(0.0,n),1.0])
;Use NNLS to minimize ||Eu-f||
u=nnls(E,f)

;E1=Transpose(E)
;f1=Transpose(f)
;se=size(E,/dim)
;BND=[replicate(0,1,se[0]),replicate(1e6,1,se[0])]
;bvls,E1,f1,BND,u
r=E##u-f
;Test for consistency
;IF Total(Abs(r)) GT 0 THEN BEGIN
;	phi=1
;	return,transpose(-r[0:n-1]/r[n])
;ENDIF ELSE BEGIN
;	phi=0
;	return,-1
;ENDELSE
phi=total(abs(r))
xans=transpose(-r[0:n-1]/r[n])

;Test using nnls_hessi
;se=size(E,/dim)
;NNLS_HESSI,E,se[1],se[0],f,x,rnorm,w,indx,mode

return,transpose(-r[0:n-1]/r[n])
END