FUNCTION LSIE,E,f,C,d,G,t,phi,TOL=tol,_EXTRA=extra
;+
;Minimize ||Ex-f|| subject to Cx=d and Gx>=h
;
;Use ortho_factor to find the C=HRK^T
;
;The parameter phi is returned by LDP, which is the
;search routine. If phi is very small, the solution may
;be faulty.
;
;H. Rhody
;June, 2004
;
;-

IF n_elements(tol) LE 0 THEN tol=-1e-3

;STEP 0: PARAMETERS AND ERROR CHECKING
se=size(E)
IF se[0] NE 2 THEN MESSAGE,'E array must be 2D'
m2=se[2]
n=se[1]
IF n_elements(f) NE m2 THEN MESSAGE, 'E and f of incompatible sizes'
f=reform(f,1,m2)
sc=size(C)
IF sc[1] NE n THEN MESSAGE, 'E and C must have same column dimension'
sg=size(G)
IF sg[1] NE n THEN MESSAGE, 'E and G must have same column dimension'


;STEP 1: FACTOR C
rc=Ortho_Factor(C,H,R,K,_EXTRA=extra)
TEMP=C##K
R1=R[0:rc-1,0:rc-1]
H_1=H[0:rc-1,*]

;STEP 2: COMPUTE y1
y1=invert(R1)##Transpose(H_1)##d

;IF rc=n then there are no free variables. Entire
;solution is determined by the equality constraints
;All there is to do is to check whether the inequality
;constraints are satisfied.
IF rc EQ n THEN BEGIN
	x=K##y1
	IF MIN(G##x-t) GE tol THEN RETURN,x ELSE RETURN,-1
ENDIF

;STEP 3: COMPUTE E##K^T=[E1,E2] partition
TEMP=E##K
E1=TEMP[0:rc-1,*]
E2=TEMP[rc:*,*]
f1=f-E1##y1

;STEP 4: COMPUTE G##K^T=[G1,G2] partition
TEMP=G##K
G1=TEMP[0:rc-1,*]
G2=TEMP[rc:*,*]
h1=t-G1##y1

;We now have an LSI problem to minimize ||E2##y2-f1||
;subject to G2##y2>=h1. We can reduce this problem
;to LDP.

;STEP 5: FACTOR E2
re=ortho_factor(E2,U,P,V,_EXTRA=extra)
P1=P[0:re-1,0:re-1]

;STEP 6: PARTITION THE SOLUTION
U1=U[0:re-1,*]
;U2=U[re:*,*]
V1=V[0:re-1,*]
;V2=V[re:*,*]
P1i=Invert(P1)

;STEP 7: SET UP LDP
G3=G2##V1##P1i##Transpose(U1)
h3=(h1-G2##V1##P1i##Transpose(U1)##f1)
s=HR_LDP(G3,h3,phi)
;IF phi EQ 0 THEN MESSAGE,'Inconsistent Inequality Constraints'

;STEP 8: UNWRAP THE SOLUTION
z=P1i##Transpose(U1)##(s+f1)

y2=V1##z
x=K##[[y1],[y2]]

RETURN,x
END