[571] | 1 | function [LD,D]=ldform(AD,D0) |
---|
| 2 | % transform general decomposition A 'D0 A into L'DL decomposition with |
---|
| 3 | % lower triangular(trapeziod) L with unites on diagonal by Householder transformation. |
---|
| 4 | % Both input and output can be given either in a "packed form" |
---|
| 5 | % as one matrix AD resp. LD or as a pair of matrices (A,D0) resp. (L,D). |
---|
| 6 | % The new decomposition fullfills: A'*D0*A = L'*D*L |
---|
| 7 | % |
---|
| 8 | % [L, D] = ldform(A,D0) |
---|
| 9 | % [L, D] = ldform(AD) ( [A,D0]=ld2ld(AD) ) |
---|
| 10 | % [LD] = ldform(A,D0) |
---|
| 11 | % [LD] = ldform(AD) ( [A,D0]=ld2ld(AD) ) |
---|
| 12 | % |
---|
| 13 | % AD = an arbitrary dimension matrix with squares of row weights on its |
---|
| 14 | % diagonal. |
---|
| 15 | % A = an arbitrary dimension matrix |
---|
| 16 | % D0 = diagonal matrix containing squares of weights of corresponding |
---|
| 17 | % rows of L0 |
---|
| 18 | % L = lower triangular (trapeziod) matrix normalized to unit diagonal |
---|
| 19 | % D = diagonal matrix |
---|
| 20 | % LD = lower triangular matrix with unit diagonal replaced by diagonal |
---|
| 21 | % D |
---|
| 22 | % |
---|
| 23 | % ----- -- |
---|
| 24 | % | | --> | \ |
---|
| 25 | % | | | \ |
---|
| 26 | % ------ ------ |
---|
| 27 | % |
---|
| 28 | % Design : J. Bohm |
---|
| 29 | % Project: post ProDaCTools |
---|
| 30 | |
---|
| 31 | % Calls : ld2ld |
---|
| 32 | |
---|
| 33 | if(nargin<1) error('At least one input must be given'); end; |
---|
| 34 | |
---|
| 35 | [m,n]=size(AD); |
---|
| 36 | |
---|
| 37 | if(nargin<2) |
---|
| 38 | if(m~=n) |
---|
| 39 | error('packed input has sense only for square matrices'); |
---|
| 40 | end; |
---|
| 41 | [A,D0]=ld2ld(AD); |
---|
| 42 | else |
---|
| 43 | A=AD; |
---|
| 44 | end; |
---|
| 45 | |
---|
| 46 | if((nargout==1)&(m~=n)) |
---|
| 47 | error('packed output has sense only for square matrices'); |
---|
| 48 | end; |
---|
| 49 | |
---|
| 50 | if(size(D0,1)~=m) error('matrix D0 must have the same row count as A'); end; |
---|
| 51 | |
---|
| 52 | v=zeros(1,m); |
---|
| 53 | |
---|
| 54 | dd=sqrt(D0); |
---|
| 55 | L=dd*A; |
---|
| 56 | D=D0; |
---|
| 57 | mn=min(m,n); |
---|
| 58 | |
---|
| 59 | for i=n:-1:n-mn+1, |
---|
| 60 | v=zeros(1,m); |
---|
| 61 | sum=L(1:m-n+i,i)'*L(1:m-n+i,i); |
---|
| 62 | if(L(m-n+i,i)==0) |
---|
| 63 | beta=sqrt(sum); |
---|
| 64 | else |
---|
| 65 | beta= L(m-n+i,i)+sign(L(m-n+i,i))*sqrt(sum); |
---|
| 66 | end; |
---|
| 67 | |
---|
| 68 | if abs(beta)<eps,D(m-n+i,m-n+i)=0;L(m-n+i,i)=1;continue; end |
---|
| 69 | |
---|
| 70 | v(1:m-n+i-1)=L(1:m-n+i-1,i)/beta; |
---|
| 71 | v(m-n+i)=1; |
---|
| 72 | w=zeros(1,n); |
---|
| 73 | sum=v(1:m-n+i)*v(1:m-n+i)'; |
---|
| 74 | |
---|
| 75 | pom=-2/sum; |
---|
| 76 | for j=i:-1:1 |
---|
| 77 | w(j)=v(1:m-n+i)*L(1:m-n+i,j)*pom; |
---|
| 78 | end |
---|
| 79 | |
---|
| 80 | |
---|
| 81 | L(1:m-n+i,1:i-1)=L(1:m-n+i,1:i-1)+v(1:m-n+i)'*w(1:i-1); |
---|
| 82 | L(1:m-n+i-1,i)=0; |
---|
| 83 | L(m-n+i,i)=L(m-n+i,i)+v(m-n+i)*w(i); |
---|
| 84 | |
---|
| 85 | D(m-n+i,m-n+i)=L(m-n+i,i)^2; |
---|
| 86 | L(m-n+i,:)= L(m-n+i,:)/L(m-n+i,i); |
---|
| 87 | end |
---|
| 88 | |
---|
| 89 | for(i=1:m-n) |
---|
| 90 | D(i,i)=0; |
---|
| 91 | end; |
---|
| 92 | |
---|
| 93 | if(nargout==1) |
---|
| 94 | LD=ld2ld(L,D); |
---|
| 95 | else |
---|
| 96 | LD=L; |
---|
| 97 | end; |
---|