function [mu,d,nsteps,grad]=ssv(a,k,d,max,tol,verbose);
%SSV Structured singular value (at one frequency).   
%       [mu,d,nsteps,grad] = ssv(a,k,d,max,tol,verbose)
%                       mu = ssv(a,k)
% Calculates the structured singular value (mu) of a square complex
% matrix. Block perturbations are allowed, but not correlated.
% Input arguments:
%       a  nxn complex matrix whose SSV is to be computed.
%       k  vector of m block sizes, such that sum(k) == n.
%                                       ( Default: ones(1,n) )
%       d  vector of m initial scaling numbers. ( Default: ones(1,m) )
%     max  maximum number of hill-climbing steps. (Default: 100)
%     tol  convergence criterion. (Default: 1e-6)
% verbose  1 if diagnostic information to be output, 0 otherwise.
%                                                       (Default: 1)
% Output arguments:
%      mu  an upper bound for the structured singular value
%       d  vector of m optimal scaling numbers. ( d(1) == 1 always)
%  nsteps  number of hill-climbing steps taken
%    grad  gradient vector at termination

% Copyright (C) 1989 by Cambridge Control Ltd. J.M.Maciejowski,7.8.1989.
% Corrected by JMM, 8.9.89.
% Acknowledgement:
% This function is based on a file supplied by M.Steinbuch of Philips
% Research Lab., Eindhoven, The Netherlands, which was in turn based on
% a NASA MatrixX routine. It has been subsequently modified extensively
% by J.M.Maciejowski at Cambridge Control Ltd to allow non-scalar blocks,
% to improve the interface and error-checking, the speed, the termination
% criteria, and the documentation.

%%%%%  Check input arguments and assign defaults:
nargchk(nargin,1,6);
[n,nc]=size(a); if n ~= nc, error('Input matrix must be square'), end;
if nargin<6, verbose = 1; end;
if nargin<5, tol = 1e-6; end;
if nargin<4, max = 100; end; % NB 'max' is now a variable,not function!
if nargin<3, 
  if nargin<2, k = ones(1,n); end;
  d = ones(1,length(k));
end;
if any(k<1), error('Block sizes must be positive integers'), end;
if sum(k)~=n, error('Sum of block sizes must equal size of matrix'), end;
if length(d)~=length(k), 
  error('Length of scaling vector (d) must equal the number of blocks'),
end
if any(d==0), error('Zero elements not allowed in d vector'), end;
if any(imag(d)~=0) | any(real(d)<0),
  d = abs(d);
  if verbose, 
    disp('Warning: Imaginary or negative elements in d vector')
    disp('         have been replaced by their magnitudes.   ')
  end
end
%%%%% End of input argument checks
%
%%%%% Deal with the 1 block case: (This added 8.9.89.)
if length(k) == 1, 
  mu = svd(a);  mu = mu(1);
  d = 1;  nsteps = 0;  grad = [];
  if verbose, 
    disp('SSV: Only 1 block. No hill climbing required.')
    mu
  end;
  return;
end
%%%%% End of 1 block case
%
%%%%% Define some variables:
m = length(k);  % m is the number of blocks
d=d(:)/d(1);  % d is now a column vector with d(1)==1
maxsigo=10000;  % Initial upper bound for SSV
perf=100;  % Initial value of convergence indicator
continue=1;  icount=0;
if verbose,  disp('SSV:   Step     mu     norm(grad)     perf'); end;
%%%%% End of initial definitions
%
%%%%% Begin hill-climbing loop:
while continue,
  %%%%% Begin gradient computation:
  temp = [];  grad = [];  % Clear variables 
  for i = 1:m, temp = [temp; d(i)*ones(k(i),1)]; end;
  dd = diag(temp);  % dd is now an nxn diagonal matrix
  dda=dd*a;  addi=a/dd;  dadi=dd*addi;
  [u,svv,v]=svd(dadi);
  maxsig=svv(1);
  for kk=2:m,   % This FOR loop assumes svv(1)~=svv(2). 
		% See Freudenberg etal (Int.Jnl.Contr. 1982).
	index = [sum(k(1:kk-1))+1 : sum(k(1:kk))];
	onevec = ones(length(index),1);
	ddk = zeros(n,1); didk = ddk;
	ddk(index) = onevec;   didk(index) = -onevec*d(kk)^(-2);
	ddk = diag(ddk);  didk = diag(didk);
	gradk=real(u(:,1)'*(ddk*addi+dda*didk)*v(:,1));
	grad=[grad;gradk]; 
  end;
  %%%%% End of gradient computation
  %
  if verbose,
    fprintf('       %3.0f  %10.3e %10.3e',icount,maxsig,norm(grad));
    fprintf('  %10.3e\n',perf);
  end;
  icount = icount + 1;
  %
  %%%%% Begin Quasi-Newton step:
  %%%%% Begin determination of search direction:
  if icount==1, newhess = 1; maxsigo = maxsig; end;
  if newhess,   % Reset Hessian, take one steepest descent step.
    hh=eye(m-1);  newhess = 0; 
  else
    y=grad-grado;
    temp=norm(y);
    if temp < 0.1e-12   
      [dum1,dum2]=size(y);
      y=0.1e-12*ones(dum1,dum2);
    end;
    aa=(deld*deld')/(deld'*y);
    bb=-(hh*y*y'*hh)/(y'*hh*y);
    hh=hh+aa+bb;
  end; 
  delds=-hh*grad;
  grado=grad;
  %%%%% End of search direction determination
  %
  %%%%% Begin line search:
  kappa=1;
  % Try to avoid nonpositive weights:
  dm=d(2:m)+kappa*delds; % d(1) remains unchanged
  [mindm,h] = min(dm);
  if mindm <= 0,
	kappa=abs(.5*d(h+1)/delds(h)); % New step size 
  end;
  % Line search until maxsig is at least as small as before:
  tryagain = 1;
  while tryagain == 1,
    temp = [];
    dd=[1;d(2:m)+kappa*delds];
    for i = 1:m, temp = [temp; dd(i)*ones(k(i),1)]; end;
    dd = diag(temp);
    dda=dd*a; addi=a/dd; dadi=dd*addi;
    svv=svd(dadi);  maxsig=svv(1);
    if maxsig > maxsigo,
      kappa=kappa/2;
    else
      tryagain = 0;
    end;
  end;
  %%%%% End of line search
  %
  %%%%% Update scaling parameter vector d:
  deld=kappa*delds;
  d(2:m) = d(2:m)+deld;
  %
  %%%%% Check for convergence:
  perf=(maxsigo-maxsig)/maxsig;
  maxsigo=maxsig;
  if icount>m,
    if perf < tol,  % See 'Practical Optimization', Gill,Murray and
      if norm(grad) <= tol^(1/3)*(1+maxsig),     % Wright, page 306
	continue=0;
      else           % Line search can't find a much better point, but
	newhess = 1; % the gradient not yet small. So try a different
      end;           % direction - try one steepest descent step.
    end;
  end;
  if icount==max, continue=0; end;  % Terminate without convergence
end;
%%%%% End of hill-climbing loop
%
%%%%% Finish off:
mu = maxsig;
nsteps = icount;
if verbose,
  fprintf('       %3.0f  %10.3e     NaN   ',nsteps,mu);
  fprintf('  %10.3e\n',perf);
  if nsteps==max, disp('SSV: Not converged yet'),end
end;