home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / CHEAPL.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.4 KB  |  60 lines
  1. function [lb,fulpert,rowp] = cheapl(mat,blk,blkp)
  2.  [nblk dum] = size(blk);
  3.  [nr nc] = size(mat);
  4.  new = [];
  5.  for i=1:nblk
  6.    new = [new ; mat(blkp(i,2):blkp(i+1,2)-1,1:nc)];
  7.    if blk(i,2) ~= 0
  8.      if blk(i,1) > blk(i,2)
  9.        zer = zeros(blk(i,1)-blk(i,2),nc);
  10.        new = [new ; zer];
  11.      end
  12.    end
  13.  end
  14.  [nr nc]=size(new);
  15.  out =[];
  16.  for i=1:nblk
  17.    out = [out  new(1:nr,blkp(i,1):blkp(i+1,1)-1)];
  18.    if blk(i,2) ~= 0
  19.      if blk(i,2) > blk(i,1)
  20.        zer = zeros(nr,blk(i,2)-blk(i,1));
  21.        out = [out zer];
  22.      end
  23.    end
  24.  end
  25.  eigv=eig(out);
  26.  lb = max(abs(eigv));
  27.  if lb == 0
  28.    disp(['LOWER BOUND FOR mu(M) is ZERO,']);
  29.    disp(['  OFFENDING PERTURBATION IS VERY']);
  30.    disp(['  LARGE TO INDICATE THIS']);
  31.    fulpert = randel(blk,1e12);
  32.    rowl = 0;
  33.    for i=1:nblk
  34.      if blk(i,2) == 0
  35.        rowl = rowl + 1;
  36.      else
  37.        rowl = rowl + blk(i,1)*blk(i,2);
  38.      end
  39.    end
  40.    rowp = 1e12*ones(1,rowl);
  41.  else
  42.    loc = find(lb==abs(eigv));
  43.    llb = 1/eigv(loc(1));
  44.    fulpert = [];
  45.    rowp = [];
  46.    for i=1:nblk
  47.      if blk(i,2) == 0
  48.        fulpert = daug(fulpert,(llb)*eye(blk(i,1)));
  49.        rowp = [rowp llb];
  50.      else
  51.        tmp = (llb)*eye(blk(i,1),blk(i,2));
  52.        fulpert = daug(fulpert,tmp);
  53.        tmpp = reshape(tmp.',blk(i,1)*blk(i,2),1);
  54.        rowp = [rowp tmpp.'];
  55.      end
  56.    end
  57.  end
  58. %
  59. % Copyright MUSYN INC 1991,  All Rights Reserved
  60.