function [fun,grad,hess] = fghbarrier(v, pars);
%  Compute the function, gradient and Hessian of the barrier function
%  for the robust spectral abscissa, i.e.,
%    gam - epsln log det(2 gam P - (P A(x) + A(x)'P))
%        - epsln log det(I - P) - epsln log det(P -  delta I)
%        - epsln log (xbound^2 - x'x)
%  wrt variables gam, x and svec(P)  (svec maps symmetric matrices to vectors)
%  Here epsln is the fixed barrier parameter,
%  delta is the fixed robustness parameter
%  and xbound is a fixed bound on the norm of x.
%  Input:
%     v collects the 3 variables together:  v = [gam; x; svec(P)]
%     pars is a structure with the necessary parameters in its fields
%       pars.m   length of x (number of variables in affine function)
%       pars.n   dimension of symmetric matrix P
%       pars.epsln  barrier parameter, positive
%       pars.delta  robust parameter, in [0,1)   
%       pars.A    cell array defining the affine function, A{m+1} is const term
%  Output:
%     fun and grad are the barrier function and gradient 
%     hess is the Hessian of the barrier function.  It is only computed if
%       it is requested, so the routine may also be used to efficiently compute
%       only fun and grad.
%  Supporting mfiles needed: svec.m, skron.m, sextract.m (which needs smat.m)
%    svec maps symmetric matrices of dimension m to vectors of length n(n+1)/2
%    skron is the symmetrized Kronecker product 
%      (see Alizadeh, Haeberly, Overton, SIAM J Optim 8 (1998), pp. 46-768)
%    sextract extracts gam, x and P from v = [gam; x; svec(P)]
%    smat is the inverse of svec
%  Reference: Burke, Lewis, Overton, 
%    "Two Numerical Methods for Optimizing Matrix Stability" (Appendix)
%     Lin Alg Appl (special issue on linear algebra in control), to appear in 2002.
%     Note: the derivation in the Appendix is quite compact.
%     The correspondence between that derivation and this m-file is not immediate, 
%     but it has been very carefully checked.
%  Copyright J.V. Burke, A.S. Lewis, M.L. Overton, 2001.

[gam, x, P] = sextract(v, pars.m, pars.n);   %  separate the variables

epsln = pars.epsln; delta = pars.delta;      %  access the parameters
A = pars.A; m = pars.m; n = pars.n;  xbndsq = pars.xbound^2;

%  evaluate the affine matrix function at x, using 
%    Abig, which stores the A{k} as row vectors)
%    however, note that the A{k} are still explicitly accessed below
B = pars.A0 + reshape(pars.Abig'*x, pars.n, pars.n);  % B = getaff(x,A);

%  evaluate the barrier function, checking the matrices are positive definite
M = 2*gam*P - P*B - B'*P;
I = eye(size(B));
N1 = P - delta*I;     % corresponds to inv(Delta) in Appendix
N2 = I - P;           % corresponds to inv(Lambda) in Appendix
N3 = xbndsq - x'*x;   % corresponds to 1/beta in Appendix
% first check positive definiteness via Cholesky
[L0,ind0] = chol(M);
[L1,ind1] = chol(N1);
[L2,ind2] = chol(N2);
if ind0 > 0 | ind1 > 0 | ind2 > 0 ...   % if any of the Cholesky's failed...
 | N3 <= 0                               % or if x'x >= xbndsq
   fun = inf;
   grad = [];
   hess = [];
else
   detprod = det(L0)*det(L1)*det(L2);
   fun = gam - epsln*(2*log(detprod) + log(N3));
% gradient of log barrier function wrt gam, x and P
   Q = inv(M);    % grad log det...
   C = Q*P*Q;
   gradgam = 1 - 2*epsln*inprod(P,Q);
   N1inv = inv(N1);
   N2inv = inv(N2);
   gradP = - svec(epsln*(2*gam*Q - B*Q - Q*B' + N1inv - N2inv), n);
   gradx = zeros(m,1);
   for k = 1:m
      T{k} = P*A{k}+A{k}'*P;      % called P{k} in Appendix
      gradx(k) = epsln*(inprod(Q,T{k}) + 2*x(k)/N3);
   end
   grad = [gradgam; gradx; gradP];
% now the Hessian, only if requested
   if nargout > 2
      Hgamgam = 4*inprod(C,P);
      HgamP = 2*svec(2*gam*C - B*C - C*B' - Q, n)';
      Hgamx = zeros(1,m);
      HxP = zeros(m,n*(n+1)/2);
      Hxx = (2/N3)*(eye(m) + (2/N3)*x*x');   % Hessian of log(xbndsq - x'x)
      for k = 1:m
         D{k} = Q*T{k}*Q;            % called C{k} in original notes
         Hgamx(k) = -2*inprod(C,T{k});
         HxP(k,:) = svec(-2*gam*D{k} + B*D{k} + D{k}*B' +  A{k}*Q+Q*A{k}', n)';
         for j = 1:k
            Hxx(k,j) = Hxx(k,j) + inprod(D{k},T{j});   % lower triangular part
         end
      end
      Hxx = Hxx + tril(Hxx,-1)';   % upper triangular part is same as lower
      BQ = B*Q;
      HPP = 4*gam^2*skron(Q,Q)-4*gam*skron(BQ+BQ',Q)+2*skron(BQ,BQ')...
        + 2*skron(BQ*B',Q) + skron(N1inv,N1inv) + skron(N2inv,N2inv);
%  Note: a careless reading of the derivation in the paper suggests wrongly that 
%  the second 4 in HPP should be a 2 - it really is a 4!
      hess = epsln*[Hgamgam  Hgamx  HgamP
                    Hgamx'   Hxx    HxP
                    HgamP'   HxP'   HPP ];
      hess = 0.5*(hess + hess');
   end
end
