home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / AXXBC.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  991 b   |  39 lines
  1. % function [x]=axxbc(a,b,c)
  2. %   solves the equation a*x+x*b=c, using the function SYLV.
  3.  
  4. function [x]=axxbc(a,b,c)
  5. if nargin < 3
  6.   disp(['usage: x = axxbc(a,b,c)']);
  7.   return
  8. end
  9. [n,p]=size(a);
  10. [q,m]=size(b);if (p~=n)|(q~=m),error('dimensions wrong in axxbc');end;
  11. [p,q]=size(c);if (p~=n)|(q~=m),error('dimensions wrong in axxbc');end;
  12. [u,a]=schur(a');a=a';
  13. [v,b]=hess(b);
  14. c=u'*c*v;
  15. x=zeros(n,m);
  16. bl=1;bsz=[];
  17. for i=1:n-1,
  18.   if abs(a(i,i+1))<eps,bsz=[bsz,bl];bl=1;
  19.   else bl=2;end
  20. end
  21. bsz=[bsz,bl];
  22. [ss,s]=size(bsz);bst=zeros(1,s);
  23. bst(1)=1;
  24. for i=1:s-1,bst(i+1)=bst(i)+bsz(i);end
  25. if bsz(1)==1,x(1,:)=c(1,:)/(a(1,1)*eye(m)+b);...
  26. else y=sylv(a(1:2,1:2),b,c(1:2,:));x(1:2,:)=y;end
  27. for i=2:s,
  28.   t=bst(i);
  29.   if bsz(i)==1,x(t,:)=(c(t,:)-a(t,1:t-1)...
  30.     *x(1:t-1,:))/(a(t,t)*eye(m)+b);...
  31.   else y=sylv(a(t:t+1,t:t+1),b,...
  32.     c(t:t+1,:)-a(t:t+1,1:t-1)*x(1:t-1,:));
  33.     x(t:t+1,:)=y;
  34.   end
  35. end
  36. x=u*x*v';
  37. %
  38. % Copyright MUSYN INC 1991,  All Rights Reserved
  39.