home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / MU_COMP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.9 KB  |  105 lines
  1. %function [bnd,dl,dri,rowp,sens,pinfo] = ...
  2. % mu_comp(mat,nblk,blk,blkp,widp,repp,fulp,rnblk,rblk,rblkp,opt,dleft,drighti)
  3. %
  4. %    Main computational routine for mu.  Not intended to be called by user.
  5. function [bnd,dl,dri,rowp,sens,pinfo] = ...
  6.   mu_comp(mat,nblk,blk,blkp,widp,repp,fulp,rnblk,rblk,rblkp,opt,dleft,drighti)
  7.  
  8. %              initialize some constants
  9.   [nr,nc]=size(mat);
  10.   restart = 0;
  11.   maxrestart = 2;
  12.   lbubok = 0;
  13.  
  14. % scale mat with the supplied d scalings
  15.  
  16.   if exist('drighti')
  17.     mmat = dleft*mat*drighti;
  18.      else
  19.     mmat = mat;
  20.     dleft = eye(nr); drighti = eye(nc);
  21.   end %if
  22.   bestup = norm(mmat);
  23.  
  24. % approximately minimize norm of dmd^(-1)
  25.  
  26.   if ~any(opt=='L')
  27.     [dl,dri,sm,sens] = ...
  28.        ap_up_bd(mmat,nblk,blk,blkp,repp,rnblk,rblk,rblkp,opt);
  29.     nsm = norm(sm);
  30.     if bestup <= nsm
  31.       sm = mmat;         dl = dleft;   dri= drighti;
  32.     else
  33.       bestup = nsm;    dl = dl*dleft;    dri = drighti*dri;
  34.     end
  35.     if any(opt=='c')
  36.       if length(repp) == 0
  37.         [sgdleft,sgdrighti] = mainsubg(sm,blk,3);
  38.       else
  39.         [sgdleft,sgdrighti] = mainsubg(sm,rblk,3);
  40.       end % if length
  41.         sm = sgdleft*sm*sgdrighti;
  42.         dl = sgdleft*dl;      dri = dri * sgdrighti;
  43.         bestup = norm(sm); 
  44.     end % if any(opt=='c')
  45.   else
  46.     dl = dleft;
  47.     dri = drighti;
  48.     sens = ones(1,nblk);
  49.     sm = mmat;
  50.   end
  51.   
  52.  while lbubok == 0
  53.    if ~any(opt=='l' | opt=='t' | opt=='L')
  54.      bnd = [bestup 0];
  55.      lbubok = 1;
  56.    else
  57.      [u,s,v] = svd(sm);
  58.      up = s(1,1);
  59.      if restart == 0 & ~any(opt=='R')
  60.        sv = v(:,1);
  61.      else
  62.        sv = crand(nc,1);
  63.        sv = (1/norm(sv))*sv;
  64.      end
  65.      sv = [sv sv];        % set up initial vectors for power iteration
  66.  
  67. %    LOWER BOUND POWER METHOD
  68.  
  69.      [lb,bw,az,actit,errflg] =  mup(sm,nblk,blk,blkp,repp,fulp,sv,opt); 
  70.      if errflg == 0        % CONVERGED
  71.        rowp = mk_pert(az,bw,nblk,blk,blkp,lb);
  72.        lbubok = 1;
  73.      else
  74. %      DIDN'T CONVERGE: GET LOWER BOUND FROM RHO(M*PERT)
  75.        rowp = mk_pert(az,bw,nblk,blk,blkp,1);
  76.        [ncp,nrp] = size(sm);
  77.        [pert,fnor] = pert_mat(rowp,blk,nrp,ncp);   % this could be way off 
  78. %      PERT has norm=1, and ROWP must be divided by FNOR to get PERT
  79.        evals = eig(mat*pert);
  80.        [seval,eorder] = sort(abs(evals));
  81.        lbeig = evals(eorder(length(eorder)));
  82.        [chlb,chfp,chrp] = cheapl(mat,blk,blkp);
  83.        if abs(lbeig) > chlb
  84.          rowp = rowp/fnor/lbeig;
  85.          lbubok = 1;
  86.          lb = abs(lbeig);
  87.        elseif restart == maxrestart
  88.          rowp = chrp;
  89.          lb = chlb;
  90.          lbubok = 1;
  91.        else
  92.          restart = restart + 1;
  93.        end
  94.      end
  95.    end
  96.  
  97. %  BOUNDS, AND ERROR/CONVERGENCE DATA
  98.  
  99.    bnd = [up lb];
  100.    pinfo = [errflg actit];
  101.        
  102. end %if ~any(opt=='l' | opt =='t' | opt=='L')
  103. %
  104. % Copyright MUSYN INC 1991,  All Rights Reserved
  105.