home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / PWR_1.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.1 KB  |  78 lines
  1. function [bw,az,beta,ec] = pwr_1(mat,nblk,blk,blkp,repp,fulp,bw)
  2. %
  3. %  performs one iteration of the power iteration
  4. %    described in "A Power Method for the Structured
  5. %    Singular Value," Andy Packard, Michael Fan, and
  6. %    John Doyle, pg. 2137-2142, 1988 CDC proceedings,
  7. %    Austin, TX.
  8. %
  9. %    Not to be called by user.
  10. %
  11.     [nr nc] = size(mat);
  12.     ec = 0;
  13.     z = zeros(nr,1);
  14.     a = mat * bw(:,1);
  15.     beta(1) = norm(a);
  16.     if beta(1) < 1e-6
  17.       ec = 1;
  18.       return
  19.     end
  20.     a = (1.0/beta(1)) * a; 
  21.     for j=1:length(repp)
  22.       ww = bw(blkp(repp(j),1):blkp(repp(j)+1,1)-1,2);
  23.       aa = a(blkp(repp(j),2):blkp(repp(j)+1,2)-1);
  24.       wwaa = ww' * aa;
  25.     den = abs(wwaa);
  26.       if den < 10*eps
  27.     wwaa = wwaa + (den < 10*eps);    % set phase=1 for zero elements
  28.     den = abs(wwaa);
  29.       end %if den
  30.         z(blkp(repp(j),2):blkp(repp(j)+1,2)-1) = (wwaa/den) * ww;
  31.  
  32.     end %for
  33.     for j=1:length(fulp)
  34.       aa = a(blkp(fulp(j),2):blkp(fulp(j)+1,2)-1);
  35.       na = norm(aa);
  36.       if na < 10*eps
  37.        fix = 1;
  38.     else
  39.        fix = norm(bw(blkp(fulp(j),1):blkp(fulp(j)+1,1)-1,2)) / na;
  40.       end
  41.       z(blkp(fulp(j),2):blkp(fulp(j)+1,2)-1) = fix * aa;
  42.     end
  43.     az = [a z];
  44.  
  45.     bw(:,2) = mat' * z; 
  46.     beta(2) = norm(bw(:,2));
  47.     if beta(2) < 1e-6
  48.       ec = 1;
  49.       return
  50.     end
  51. %   check zero and return
  52.     bw(:,2) = (1.0/beta(2)) * bw(:,2);
  53.     for j=1:length(repp)
  54.       ww = bw(blkp(repp(j),1):blkp(repp(j)+1,1)-1,2);
  55.       aa = a(blkp(repp(j),2):blkp(repp(j)+1,2)-1);
  56.       aaww = aa' * ww;
  57.     den = abs(aaww);
  58.       if den < 10*eps
  59.     aaww = aaww + (den < 10*eps);    % set phase=1 for zero elements
  60.     den = abs(aaww);
  61.       end %if den
  62.     bw(blkp(repp(j),1):blkp(repp(j)+1,1)-1,1) =  (aaww/den)* aa;
  63.     end
  64.     for j=1:length(fulp)
  65.       ww = bw(blkp(fulp(j),1):blkp(fulp(j)+1,1)-1,2);
  66.       nw = norm(ww);
  67.       if nw < 10*eps
  68.        fix = 1;
  69.     else
  70.        fix = norm(a(blkp(fulp(j),2):blkp(fulp(j)+1,2)-1)) / nw ;
  71.       end
  72. %     check zero and return
  73.       bw(blkp(fulp(j),1):blkp(fulp(j)+1,1)-1,1) = fix * ww;
  74.     end
  75.  
  76. %
  77. % Copyright MUSYN INC 1991,  All Rights Reserved
  78.