root/library/utia_legacy/ticket_12/ldform.m @ 653

Revision 571, 2.3 kB (checked in by smidl, 15 years ago)

matlab test for #12

Line 
1function [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
33if(nargin<1) error('At least one input must be given'); end;
34
35[m,n]=size(AD);
36
37if(nargin<2)
38    if(m~=n)
39        error('packed input has sense only for square matrices');
40    end;
41    [A,D0]=ld2ld(AD);
42else
43    A=AD;
44end;
45
46if((nargout==1)&(m~=n))
47        error('packed output has sense only for square matrices');
48end;
49
50if(size(D0,1)~=m) error('matrix D0 must have the same row count as A'); end;
51   
52v=zeros(1,m);
53
54dd=sqrt(D0);
55L=dd*A;
56D=D0;
57mn=min(m,n);
58
59for 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);
87end
88
89for(i=1:m-n)
90    D(i,i)=0;
91end;
92
93if(nargout==1)
94    LD=ld2ld(L,D);
95else
96    LD=L;
97end;
Note: See TracBrowser for help on using the browser.