[37] | 1 | function blks(p,Q_in,v_in) |
---|
| 2 | |
---|
| 3 | % Start by placing blocks in separate matrices |
---|
| 4 | [v1,dummy1,r1,dummy3]=dmperm(Q_in{1}+eye(length(Q_in{1}))); |
---|
| 5 | for blocks = 1:length(r1)-1 |
---|
| 6 | i1 = r1(blocks); |
---|
| 7 | i2 = r1(blocks+1)-1; |
---|
| 8 | Q{blocks} = Q_in{1}(i1:i2,i1:i2); |
---|
| 9 | v{blocks} = v_in{1}(i1:i2); |
---|
| 10 | S{blocks} = eye(length(v{blocks})); |
---|
| 11 | end |
---|
| 12 | |
---|
| 13 | pdecomp = 0; |
---|
| 14 | F = set([]); |
---|
| 15 | obj = 0; |
---|
| 16 | |
---|
| 17 | for blocks = 1:length(Q) |
---|
| 18 | if length(Q{blocks})> 1 |
---|
| 19 | if blocks <= 100 |
---|
| 20 | [Si,index] = blkOne(Q{blocks}); |
---|
| 21 | else |
---|
| 22 | Si = eye(length(Q{blocks})); |
---|
| 23 | end |
---|
| 24 | else |
---|
| 25 | Si = 1; |
---|
| 26 | index = 1; |
---|
| 27 | end |
---|
| 28 | Qreduced{blocks} = sdpvar(size(Si,2)); |
---|
| 29 | F = F + set(Qreduced{blocks} >= 0); |
---|
| 30 | S{blocks} = S{blocks}*Si; |
---|
| 31 | vi = S{blocks}'*v{blocks}; |
---|
| 32 | pdecomp = pdecomp + vi'*Qreduced{blocks}*vi; |
---|
| 33 | obj = obj + trace(Qreduced{blocks}); |
---|
| 34 | end |
---|
| 35 | |
---|
| 36 | solvesdp(F+set(coefficients(p-pdecomp,recover(depends(p))) == 0),obj,sdpsettings('dualize',1)); |
---|
| 37 | |
---|
| 38 | for i = 1:length(Q) |
---|
| 39 | Q{i} = double(Qreduced{i}); |
---|
| 40 | end |
---|
| 41 | |
---|
| 42 | clear ss |
---|
| 43 | for i = 1:length(Q) |
---|
| 44 | ss(i) = norm(Q{i}) > 1e-4; |
---|
| 45 | end |
---|
| 46 | Q = {Q{ss}}; |
---|
| 47 | S = {S{ss}}; |
---|
| 48 | v = {v{ss}}; |
---|
| 49 | |
---|
| 50 | F |
---|
| 51 | |
---|
| 52 | |
---|
| 53 | function [S,index] = blkOne(T) |
---|
| 54 | |
---|
| 55 | s = 1; |
---|
| 56 | S = []; |
---|
| 57 | clear v |
---|
| 58 | for i = 1:size(T,1) |
---|
| 59 | for j = i+1:size(T,1) |
---|
| 60 | v(i,j) = (T(i,:)*T(j,:)'./(norm(T(i,:))*norm(T(j,:))));%var(T(i,:)./T(j,:))/abs(mean(T(i,:)./T(j,:))); |
---|
| 61 | end |
---|
| 62 | end |
---|
| 63 | v = abs((abs(v)-1)) < 1e-6 ; |
---|
| 64 | |
---|
| 65 | top = 1; |
---|
| 66 | S = []; |
---|
| 67 | remove = []; |
---|
| 68 | for i = 1:size(T,1) |
---|
| 69 | if any(v(:,i)) |
---|
| 70 | j = find(v(:,i)); |
---|
| 71 | for k = 1:length(j) |
---|
| 72 | nn = (norm(T(i,:))/norm(T(j(k),:))); |
---|
| 73 | s = (T(i,:)*T(j(k),:)'./(norm(T(i,:))*norm(T(j(k),:)))); |
---|
| 74 | if j(k) <= size(S,2) |
---|
| 75 | if S(j(k),j(k)) |
---|
| 76 | S(i,j(k)) = s*(round(10*nn)/10); |
---|
| 77 | S(i,j(k)) = round(S(i,j(k))*10)/10; |
---|
| 78 | remove = [remove i]; |
---|
| 79 | end |
---|
| 80 | end |
---|
| 81 | end |
---|
| 82 | else |
---|
| 83 | S(i,i) = 1; |
---|
| 84 | end |
---|
| 85 | end |
---|
| 86 | |
---|
| 87 | S = S(:,any(S,1)); |
---|
| 88 | index = setdiff(1:size(T,2),remove); |
---|
| 89 | |
---|
| 90 | TT = T; |
---|
| 91 | TT(remove,:) = []; |
---|
| 92 | TT(:,remove) = []; |
---|
| 93 | S*TT*S' - T |
---|
| 94 | |
---|
| 95 | |
---|
| 96 | function [Q,v] = prune(Q,v,Qtemp) |
---|
| 97 | |
---|
| 98 | for sosfun = 1:length(Q) |
---|
| 99 | keep = diag(Qtemp{sosfun})>tol; |
---|
| 100 | Qtemp(:,find(~keep)) = []; |
---|
| 101 | Qtemp(find(~keep),:) = []; |
---|
| 102 | |
---|
| 103 | m{sosfun} = m{sosfun}(find(keep)); |
---|
| 104 | |
---|
| 105 | Qtemp(abs(Qtemp) < tol) = 0; |
---|
| 106 | [v1,dummy1,r1,dummy3]=dmperm(Qtemp+eye(length(Qtemp))); |
---|
| 107 | lengths{sosfun} = []; |
---|
| 108 | n{sosfun} = {}; |
---|
| 109 | for blocks = 1:length(r1)-1 |
---|
| 110 | i1 = r1(blocks); |
---|
| 111 | i2 = r1(blocks+1)-1; |
---|
| 112 | if i2>i1 |
---|
| 113 | n{sosfun}{blocks} = m{sosfun}(v1(i1:i2)); |
---|
| 114 | else |
---|
| 115 | n{sosfun}{blocks} = m{sosfun}(v1(i1)); |
---|
| 116 | end |
---|
| 117 | lengths{sosfun} = [lengths{sosfun} length(n{sosfun}{blocks})]; |
---|
| 118 | end |
---|
| 119 | lengths{sosfun} = sort(lengths{sosfun}); |
---|
| 120 | end |
---|
| 121 | |
---|