0001 function v=v_dlyapsq(a,b)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028 [q,s]=schur(a');
0029 [q,s]=rsf2csf(q,s);
0030 [qd,r]=qr(b'*q,0);
0031
0032 r0=r;
0033 [m,n]=size(r);
0034 u=zeros(n,n);
0035 if m==1
0036 for i=1:n-1
0037 in=i+1:n;
0038 si=s(i,i);
0039 aa=sqrt(1-si*si');
0040 u(i,i)=r(1)/aa;
0041 u(i,in)=(u(i,i)*si'*s(i,in)+aa*r(2:end))/(eye(n-i)-si'*s(in,in));
0042 r=aa*(u(i,i)*s(i,in)+u(i,in)*s(in,in))-si*r(2:end);
0043 end
0044 u(n,n)=r/sqrt(1-s(n,n)*s(n,n)');
0045
0046 else
0047 w=zeros(m,1); w(m)=1;
0048 em=eye(m);
0049 for i=1:n-m
0050 in=i+1:n;
0051 si=s(i,i);
0052 aa=sqrt(1-si*si');
0053 u(i,i)=r(1,1)/aa;
0054 u(i,in)=(u(i,i)*si'*s(i,in)+aa*r(1,2:end))/(eye(n-i)-si'*s(in,in));
0055 vv=aa*(u(i,i)*s(i,in)+u(i,in)*s(in,in))-si*r(1,2:end);
0056 rr=zeros(m,n-i);
0057 rr(1:m-1,:)=r(2:end,2:end);
0058 [qq,r]=qrupdate(em,rr,w,vv');
0059 end
0060 for i=n-m+1:n-1
0061 in=i+1:n;
0062 si=s(i,i);
0063 aa=sqrt(1-si*si');
0064 u(i,i)=r(1,1)/aa;
0065 u(i,in)=(u(i,i)*si'*s(i,in)+aa*r(1,2:end))/(eye(n-i)-si'*s(in,in));
0066 vv=aa*(u(i,i)*s(i,in)+u(i,in)*s(in,in))-si*r(1,2:end);
0067 rr=zeros(n-i+1,n-i);
0068 rr(1:n-i,:)=r(2:end,2:end);
0069 [qq,rr]=qrupdate(eye(n-i+1),rr,w(m-n+i:end),vv');
0070 r=rr(1:n-i,:);
0071 end
0072
0073 u(n,n)=r/sqrt(1-s(n,n)*s(n,n)');
0074
0075 end
0076
0077 v=triu(qr(u*q'));
0078 dv=diag(v);
0079 ix=dv~=0;
0080 v(ix,:)=diag(abs(dv(ix))./dv(ix))*v(ix,:);
0081 if isreal(a) & isreal(b)
0082 v=real(v);
0083 end