% [a0, b0, c0, d0, Q0, gam] = evalguptri2 (a, b, c, d, Q, kstr)
% [a0, b0] = evalguptri2 (a, b, kstr)
% evalguptri2 (kstr)

function [a0, b0, c0, d0, Q0, gam] = evalguptri2 (a, b, c, d, Q, kstr, mode)

gam = 0;

tol = 1e-10;

if nargin < 7, mode = 0; end % normal, default mode
if mode ~= 0, tol = -0.001; end

if nargin < 1,
  disp ('usage: [a0, b0, c0, d0, Q0, gam] = evalguptri2 (a, b, c, d, Q, kstr)');
  return;
end
if nargin < 2, kstr = a;
elseif nargin < 4, kstr = c; end

k0 = 1;
items = 0;

if nargin > 2,
  a1 = a;
  b1 = b;
end

if nargin > 5,
  c1 = c;
  d1 = d;
  Q1 = Q;
end


% L and J blocks

for k = k0:size(kstr,2)
  if kstr(1,k) >= 0,
    n = kstr(1,k) - kstr(2,k);
    if n > 0,
      if items > 0, fprintf (', '); else items = 1; end;
      if n > 1, fprintf ('%d', n); end
      j = k-k0;
      if mode == 0,
        fprintf ('L%d', j);
        for i = 1:n,
          if nargin > 5,
            [a1, b1, c1, d1, Q1, gam] = reduce (a1, b1, c1, d1, Q1, gam, j, j+1);
           end
        end
      else
        fprintf ('L%d''', j);
      end
    end
    if kstr(1,k+1) >= 0,
      n = kstr(2,k) - kstr(1,k+1);
      if n > 0,
        if items > 0, fprintf (', '); else items = 1; end;
        if n > 1, fprintf ('%d', n); end
        j = k-k0+1;
        fprintf ('J%d(0)', j);
        for i = 1:n,
          if nargin > 5 && mode == 0,
            [a1, b1, c1, d1, Q1, gam] = reduce (a1, b1, c1, d1, Q1, gam, j, j);
          end
        end
      end
    end
  else
    k0 = k+1;
    break;
  end
end


% N blocks

for k = k0:size(kstr,2)
  if kstr(1,k) >= 0,
    n = kstr(1,k) - kstr(2,k);
    if n > 0,
      % L blocks should not appear here
      error ('Internal error in evalguptri2: unexpected L block');
    end
    if kstr(1,k+1) >= 0,
      n = kstr(2,k) - kstr(1,k+1);
      if n > 0,
        if items > 0, fprintf (', '); else items = 1; end;
        if n > 1, fprintf ('%d', n); end
        j = k-k0+1;
        fprintf ('N%d', j);
        for i = 1:n,
          if nargin > 5 && mode == 0,
            [a1, b1, c1, d1, Q1, gam] = reduce (a1, b1, c1, d1, Q1, gam, j, j);
          end
        end
      end
    end
  else
    k0 = k+1;
    break;
  end
end


% first indices of J() block
nm = sum (kstr(:,kstr(1,1:k0-1) > 0), 2);

% L' blocks

for k = k0:size(kstr,2)
  if kstr(1,k) >= 0,
    n = kstr(1,k) - kstr(2,k);
    if n > 0,
      if items > 0, fprintf (', '); else items = 1; end;
      if n > 1, fprintf ('%d', n); end
      j = k - k0;
      if mode == 0,
        fprintf ('L%d''', j);
      else
        fprintf ('L%d', j);
        if nargin > 5,
          for i = 1:n,
            [a1, b1, c1, d1, Q1, gam] = reduce (a1, b1, c1, d1, Q1, gam, j, j+1);
          end
        end
     end
    end
    if kstr(1,k+1) >= 0,
      n = kstr(2,k) - kstr(1,k+1);
      if n > 0,
        error ('Internal error in evalguptri2: unexpected N block');
      end
    end
  else
    k0 = k+1;
    break;
  end
end

if mode ~= 0,
  nm = sum (kstr(:,kstr(1,1:k0-1) > 0), 2) - nm;
end
% the third and last block should be a single square block
% J(?) blocks

k = k0;
if kstr(1,k) >= 0,
  if kstr(1,k) == 0,
    error ('kstr(1,k) == 0 not expected');
  end
  if kstr(2,k) == 0,
    error ('kstr(2,k) == 0 not expected');
  end
  if kstr(2,k) ~= kstr(1,k)
    error ('kstr(1,k) == kstr(2,k) expected');
  end
end

eigs = [];
n = max (0, kstr (1, k0));
if n > 0,
    if nargin > 2,
    % We assume that the eigenvalues are ordered so that the left hand
    % plane poles come first

    % The N blocks are located just after the square regular
    % eigenvalue blocks. Sort these to the front

      eigs = eig (a(nm(2)+(1:n),nm(1)+(1:n)), b(nm(2)+(1:n),nm(1)+(1:n)));
      [~, i] = sort (real(eigs));
      eigs = eigs(i);
      n = sum (real(eigs) < tol);
      % check that truncated eigenvalues are negative 
      if n ~= sum (real (eig (a(nm(2)+(1:n),nm(1)+(1:n)), ...
                              b(nm(2)+(1:n),nm(1)+(1:n)))) < tol),
        error ('eigenvalue order not correct');
      end
    end
  if ~isempty (eigs),
    nneg = sum (real (eigs) < -abs(tol));
    npos = sum (real (eigs) > abs(tol));
    nimag = length(eigs) - nneg -npos;
    if nneg > 0, 
      if items > 0, fprintf (', '); else items = 1; end;
      if nneg > 1, fprintf ('%d', nneg); end
      fprintf ('J1(-)');
    end
    if npos > 0, 
      if items > 0, fprintf (', '); else items = 1; end;
      if npos > 1, fprintf ('%d', npos); end
      fprintf ('J1(+)');
    end
    if nimag > 0, 
      if items > 0, fprintf (', '); else items = 1; end;
      if nimag > 1, fprintf ('%d', nimag); end
      fprintf ('J1(±j)');
    end
    fprintf ('\n');
    disp (eigs');
  else
    if items > 0, fprintf (', '); else items = 1; end;
    if n > 1, fprintf ('%d', n); end
    fprintf ('J1(?)');
    fprintf ('\n');
  end
else
  fprintf ('\n');
end

if nargin > 5,
  % Note that we cannot truncate the complex eigenvalues one
  % by one while computing gamma, since they come in complex
  % pairs. Instead we compute the minimum value of gamma first
  % for the imaginary eigenvalues and then truncate all
  % eigenvalues in the left hand plane in one batch.
  for i = 1:n,
    mu = eigs(i);
    if abs (real (mu)) < tol,
      N = null ((a1-mu*b1)');
      cn = c1'*N;
      dn = triu (qr (d1'*N));
      % disp ([mu max(svd(cn/dn))])
      gam = max (gam, max (svd (cn/dn)));
    end
  end
  % truncate real(mu) < tol
  j = find (real(eigs) >= tol);
  if ~isempty(j), j = j(1); else j = length(eigs)+1; end
  a1 = a1(j:end, j:end);
  b1 = b1(j:end, j:end);
  c1 = c1(j:end, :);
  d1 = d1(j:end, :);
  Q1 = Q1(:, j:end);
end

if nargin > 5,
  a0 = a1;
  b0 = b1;
  c0 = c1;
  d0 = d1;
  Q0 = Q1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [a, b, c, d, Q, gam] = reduce (a0, b0, c0, d0, Q0, gam0, j, k)
if j == k && (abs(a0(1,1)) > 100*eps || abs(b0(1,1)) > 100*eps),
  if abs(b0(1,1)) > eps,
    mu = a0(1,1)/b0(1,1);
    N = null ((a0-mu*b0)');
  else
    N = null (b0');
  end
  cn = c0'*N;
  dn = d0'*N;
  gam = max (gam0, max (svd (cn*pinv(dn))));
else
  gam = gam0;
end

a = a0(j+1:end, k+1:end);
b = b0(j+1:end, k+1:end);
c = c0(j+1:end, :);
d = d0(j+1:end, :);
Q = Q0(:, k+1:end);
