[37] | 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 |
---|