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