1 | function M = blockdiagmoment(obj,M) |
---|
2 | |
---|
3 | rt = []; |
---|
4 | |
---|
5 | r_old = []; |
---|
6 | go_on = 1; |
---|
7 | while go_on |
---|
8 | exponent_p = getexponentbase(obj+sum(sum(rt)),recover(depends(obj))); |
---|
9 | Msparsity = zeros(size(M{end},1)); |
---|
10 | vars = recover(depends(obj)); |
---|
11 | for i = 1:size(M{end},1) |
---|
12 | for j = i:size(M{end},1) |
---|
13 | exponent_Mij = getexponentbase(M{end}(i,j),vars); |
---|
14 | if ~isempty(findrows(exponent_p,exponent_Mij)) & sum(exponent_Mij)~=0 |
---|
15 | Msparsity(i,j) = 1; |
---|
16 | Msparsity(j,i) = 1; |
---|
17 | end |
---|
18 | end |
---|
19 | end |
---|
20 | |
---|
21 | [p,q,r,s] = dmperm(Msparsity+eye(length(Msparsity))); |
---|
22 | if isequal(r,r_old) |
---|
23 | go_on = 0; |
---|
24 | else |
---|
25 | Mblocked = []; |
---|
26 | Mrs = M{end}(p,p); |
---|
27 | MM = ones(length(M{end})); |
---|
28 | blocks = zeros(1,length(M{end})); |
---|
29 | Ms = Msparsity(p,p); |
---|
30 | for i = 1:length(r)-1 |
---|
31 | blocks(r(i+1)-r(i)) = blocks(r(i+1)-r(i)) + 1; |
---|
32 | Mtest = Ms(r(i):(r(i+1)-1),r(i):(r(i+1)-1)); |
---|
33 | if nnz(Mtest)==0 |
---|
34 | else |
---|
35 | MM = blkdiag(MM,ones(r(i+1)-r(i))); |
---|
36 | Mblocked = blkdiag(Mblocked,Mrs(r(i):(r(i+1)-1),r(i):(r(i+1)-1))); |
---|
37 | end |
---|
38 | end |
---|
39 | string = 'Blocks : '; |
---|
40 | for i = 1:length(blocks) |
---|
41 | if blocks(i)>0 |
---|
42 | string = [string num2str(i) 'x' num2str(i) '(' num2str(blocks(i)) ') ' ]; |
---|
43 | end |
---|
44 | end |
---|
45 | disp(string) |
---|
46 | rt = Mblocked(:); |
---|
47 | r_old = r; |
---|
48 | end |
---|
49 | end |
---|
50 | M{end} = Mblocked; |
---|
51 | |
---|
52 | % |
---|
53 | % mt = yalmip('monomtable'); |
---|
54 | % hash = randn(size(mt,2),1); |
---|
55 | % for i = 1:length(u{end}) |
---|
56 | % for j = 1:length(u{end}) |
---|
57 | % if i==1 |
---|
58 | % v1 = zeros(1,size(mt,2)); |
---|
59 | % else |
---|
60 | % vi = getvariables(u{end}(i)); |
---|
61 | % v1 = mt(vi,:); |
---|
62 | % end |
---|
63 | % if j==1 |
---|
64 | % v2 = zeros(1,size(mt,2)); |
---|
65 | % else |
---|
66 | % vj = getvariables(u{end}(j)); |
---|
67 | % v2 = mt(vj,:); |
---|
68 | % end |
---|
69 | % v = v1+v2; |
---|
70 | % vh = v*hash; |
---|
71 | % uiuj(i,j) = vh; |
---|
72 | % end |
---|
73 | % end |
---|
74 | |
---|