1 | function [x_equ,H,A_equ,b_equ,factors] = solveequalities(F_struc,K,unitary) |
---|
2 | %SOLVEEQUALITIES Internal function remove equality constraints |
---|
3 | |
---|
4 | % Author Johan Löfberg |
---|
5 | % $Id: solveequalities.m,v 1.12 2006/05/14 17:21:09 joloef Exp $ |
---|
6 | |
---|
7 | % Extract the inequalities |
---|
8 | A_equ = F_struc(1:K.f,2:end); |
---|
9 | b_equ = -F_struc(1:K.f,1); |
---|
10 | |
---|
11 | if nargin<3 |
---|
12 | unitary = 1; |
---|
13 | end |
---|
14 | |
---|
15 | factors = []; |
---|
16 | |
---|
17 | % Improve numerics by removing obviously |
---|
18 | % redundant constraints |
---|
19 | hash = 1+rand(1+size(A_equ,2),1); |
---|
20 | [i,j] = unique([A_equ b_equ]*hash); |
---|
21 | remove = setdiff(1:size(A_equ,1),j); |
---|
22 | if ~isempty(remove) |
---|
23 | A_equ(remove,:) = []; |
---|
24 | b_equ(remove,:) = []; |
---|
25 | end |
---|
26 | |
---|
27 | if ~unitary |
---|
28 | % Just use a crappy basis derived from A |
---|
29 | % FIX : fix dependent row case |
---|
30 | if 0 |
---|
31 | [L,U,P] = lu(A_equ'); |
---|
32 | s = find(diag(U)); |
---|
33 | r = setdiff(1:length(U),s); |
---|
34 | n = size(s); |
---|
35 | [i,j] = find(P'); |
---|
36 | H1 = A_equ(:,i(s)); |
---|
37 | H2 = A_equ(:,i(r)); |
---|
38 | x_equ = P'*(L'\(U'\b_equ)); |
---|
39 | % FIX : use L and U stupid! |
---|
40 | H = P'*[-H1\H2;eye(size(H2,2))]; |
---|
41 | |
---|
42 | else |
---|
43 | [L,U,P] = lu(A_equ'); |
---|
44 | n = max(find(diag(U))); |
---|
45 | [i,j] = find(P'); |
---|
46 | H1 = A_equ(:,i(1:n)); |
---|
47 | H2 = A_equ(:,i(n+1:end)); |
---|
48 | try |
---|
49 | x_equ = P'*linsolve(full(L'),linsolve(full(U'),full(b_equ),struct('LT',1==1)),struct('UT',1==1)); |
---|
50 | catch |
---|
51 | x_equ = A_equ\b_equ; |
---|
52 | end |
---|
53 | % FIX : use L and U stupid! |
---|
54 | H = P'*[-H1\H2;eye(size(H2,2))]; |
---|
55 | end |
---|
56 | else |
---|
57 | % Use unitary basis |
---|
58 | try |
---|
59 | [Q,R,E] = qr(full(A_equ)'); |
---|
60 | catch |
---|
61 | [Q,R,E] = qr(A_equ'); % Ouch, that big! |
---|
62 | end |
---|
63 | n = max(find(sum(abs(R),2)>1e-14*size(R,2))); |
---|
64 | |
---|
65 | Q1 = Q(:,1:n); |
---|
66 | R = R(1:n,:); |
---|
67 | x_equ = Q1*(R'\E'*b_equ); |
---|
68 | H = Q(:,n+1:end); % New basis |
---|
69 | end |
---|