[566] | 1 | function [LD, D] = ld2ld(LD, strs, strt) |
---|
| 2 | % convert L'DL decomposition described by a source structure strs |
---|
| 3 | % to another corresponding to a target structure strt |
---|
| 4 | % LD = ld2ld(LD, strs, strt) |
---|
| 5 | % LD = ld2ld(L,D) composes the factors L and D into the matrix LD |
---|
| 6 | % [L,D] = ld2ld(LD) decomposes the matrix LD into the factors L and D |
---|
| 7 | % |
---|
| 8 | % LD : L'DL decomposition of a positive definite matrix |
---|
| 9 | % LD ... D placed instead of the unit diagonal of L |
---|
| 10 | % strs: source structure |
---|
| 11 | % strt: target structure |
---|
| 12 | % |
---|
| 13 | % Designed : M. Karny, P. Nedoma |
---|
| 14 | % Updated : April 2000, MK July 2004 |
---|
| 15 | % Project : post-ProDaCTool |
---|
| 16 | |
---|
| 17 | nPsi = length(LD); |
---|
| 18 | |
---|
| 19 | % call [L,D] = ld2ld(LD) |
---|
| 20 | if nargin==1, D = diag(diag(LD)); |
---|
| 21 | LD = LD - D + eye(nPsi); |
---|
| 22 | return; |
---|
| 23 | end |
---|
| 24 | |
---|
| 25 | % call LD = ld2ld(L,D) |
---|
| 26 | if nargin==2, LD = LD + strs - eye(nPsi); |
---|
| 27 | return; |
---|
| 28 | end |
---|
| 29 | |
---|
| 30 | % call LD = ld2ld(LD, strs, strt) |
---|
| 31 | if isempty(strt), LD=LD(1,1); return; end |
---|
| 32 | dim = size(strt,2); |
---|
| 33 | lens = size(strs,2); |
---|
| 34 | if lens<length(LD)-1 | lens>length(LD) |
---|
| 35 | error('incompatible source structure and LD'); |
---|
| 36 | end |
---|
| 37 | % add a fictive output channel (1000) |
---|
| 38 | if lens<length(LD) |
---|
| 39 | strs = [ [1000;0], strs ]; |
---|
| 40 | strt = [ [1000;0], strt ]; |
---|
| 41 | lens = lens + 1; |
---|
| 42 | dim = dim + 1; |
---|
| 43 | end |
---|
| 44 | for i=dim:-1:1 |
---|
| 45 | chn = strt(1,i); |
---|
| 46 | td = strt(2,i); |
---|
| 47 | j = find(strs(1,:)==chn & strs(2,:)==td); |
---|
| 48 | %chn,td,strs, strt, j, dim, |
---|
| 49 | if isempty (j), error('target structure item is not in source'); end |
---|
| 50 | LD = ldperm(LD,j); % column j to 1th position |
---|
| 51 | strs = [strs(:,j),strs(:,1:j-1),strs(:,j+1:lens)]; % rebuild strs |
---|
| 52 | %strs, pause |
---|
| 53 | end |
---|
| 54 | % cut and transform L'DL if target is shorter |
---|
| 55 | if lens>dim |
---|
| 56 | for i=2:dim |
---|
| 57 | LD = ldperm(LD,dim); |
---|
| 58 | end |
---|
| 59 | for i=dim+1:lens |
---|
| 60 | LD = ldperm(LD,lens); |
---|
| 61 | end |
---|
| 62 | LD = LD(lens-dim+1:lens, lens-dim+1:lens); |
---|
| 63 | LD = ldperm(LD, dim); |
---|
| 64 | end |
---|
| 65 | |
---|
| 66 | if nargout==2 |
---|
| 67 | [LD, D] = ld2ld(LD); |
---|
| 68 | end |
---|