%function [bnd,dl,dri,rowp,sens,pinfo] = ...
% mu_comp(mat,nblk,blk,blkp,widp,repp,fulp,rnblk,rblk,rblkp,opt,dleft,drighti)
%
%	Main computational routine for mu.  Not intended to be called by user.
function [bnd,dl,dri,rowp,sens,pinfo] = ...
  mu_comp(mat,nblk,blk,blkp,widp,repp,fulp,rnblk,rblk,rblkp,opt,dleft,drighti)

%  			initialize some constants
  [nr,nc]=size(mat);
  restart = 0;
  maxrestart = 2;
  lbubok = 0;

% scale mat with the supplied d scalings

  if exist('drighti')
	mmat = dleft*mat*drighti;
     else
	mmat = mat;
	dleft = eye(nr); drighti = eye(nc);
  end %if
  bestup = norm(mmat);

% approximately minimize norm of dmd^(-1)

  if ~any(opt=='L')
    [dl,dri,sm,sens] = ...
       ap_up_bd(mmat,nblk,blk,blkp,repp,rnblk,rblk,rblkp,opt);
    nsm = norm(sm);
    if bestup <= nsm
      sm = mmat;	     dl = dleft;   dri= drighti;
    else
      bestup = nsm;    dl = dl*dleft;    dri = drighti*dri;
    end
    if any(opt=='c')
      if length(repp) == 0
        [sgdleft,sgdrighti] = mainsubg(sm,blk,3);
      else
        [sgdleft,sgdrighti] = mainsubg(sm,rblk,3);
      end % if length
        sm = sgdleft*sm*sgdrighti;
        dl = sgdleft*dl;      dri = dri * sgdrighti;
        bestup = norm(sm); 
    end % if any(opt=='c')
  else
    dl = dleft;
    dri = drighti;
    sens = ones(1,nblk);
    sm = mmat;
  end
  
 while lbubok == 0
   if ~any(opt=='l' | opt=='t' | opt=='L')
     bnd = [bestup 0];
     lbubok = 1;
   else
     [u,s,v] = svd(sm);
     up = s(1,1);
     if restart == 0 & ~any(opt=='R')
       sv = v(:,1);
     else
       sv = crand(nc,1);
       sv = (1/norm(sv))*sv;
     end
     sv = [sv sv];		% set up initial vectors for power iteration

%    LOWER BOUND POWER METHOD

     [lb,bw,az,actit,errflg] =  mup(sm,nblk,blk,blkp,repp,fulp,sv,opt); 
     if errflg == 0		% CONVERGED
       rowp = mk_pert(az,bw,nblk,blk,blkp,lb);
       lbubok = 1;
     else
%      DIDN'T CONVERGE: GET LOWER BOUND FROM RHO(M*PERT)
       rowp = mk_pert(az,bw,nblk,blk,blkp,1);
       [ncp,nrp] = size(sm);
       [pert,fnor] = pert_mat(rowp,blk,nrp,ncp);   % this could be way off 
%      PERT has norm=1, and ROWP must be divided by FNOR to get PERT
       evals = eig(mat*pert);
       [seval,eorder] = sort(abs(evals));
       lbeig = evals(eorder(length(eorder)));
       [chlb,chfp,chrp] = cheapl(mat,blk,blkp);
       if abs(lbeig) > chlb
         rowp = rowp/fnor/lbeig;
         lbubok = 1;
         lb = abs(lbeig);
       elseif restart == maxrestart
         rowp = chrp;
         lb = chlb;
         lbubok = 1;
       else
         restart = restart + 1;
       end
     end
   end

%  BOUNDS, AND ERROR/CONVERGENCE DATA

   bnd = [up lb];
   pinfo = [errflg actit];
       
end %if ~any(opt=='l' | opt =='t' | opt=='L')
%
% Copyright MUSYN INC 1991,  All Rights Reserved
