% Given a matrix pencil (A - lambda B) with lambda(A,B) = C (i.e.
% (A - lambda B) is a column singular pencil), cspsolve(A,B)computes
% a basis of polynomial solutions x_1(lambda), ... , x_p(lambda)
% for the eigenproblem (A - lambda B) x = 0. It uses the helper
% function computeHi.

% Output description:
%
% X is an (epsilon(p)+1) by p matrix, the jth column of which represents
% the jth polynomial solution x_j(lambda) in R^n[lambda]. All entries
% in column j after the first (epsilon(j)+1)n rows are zero. If the
% possibly nonzero rows of column j are partitioned into epsilon(j)+1
% column vectors x_0,j (row 1 to n), x_1,j (row n+1 to 2n), ... ,
% x_epsilon(j),j (row epsilon(j)n+1 to (epsilon(j)+1)n of length n,
% then the jth basis solution is
%
%                           epsilon(j)
%             x_j(lambda) =   sum      x_k,j lambda^k.
%                             k=0
%
% epsilon(j), j=1..p, is the degree of the jth basis solution. Here 
% epsilon(1) <= epsilon(2) <= ... <= epsilon(p).
% 
% p is the number of basis solutions.
%
% epsilonhat is a vector of the u unique degrees in epsilon. Here
% epsilonhat(1) < epsilonhat(2) < ... < epsilonhat(u).
%
% rho(i), i=1..u, is the multiplicity of the degree in epsilonhat(i).
% Note that p is the sum of the rho(i)'s.
%
% u is the number of unique degrees in epsilon.

function [X, epsilon, p, epsilonhat, rho, u] = cspsolve(A,B)

[m,n] = size(A);
minmn = min(m,n);
epsilon = zeros(1,minmn);
epsilonhat = zeros(1,minmn); rho = zeros(1,minmn);
Tk = []; X = [];

i = 1; k = -1; rk = -1; p = 0; q = 0;
while (1)
  while ((rk <= (k+1)*p-q) & (k <= min(m,n)))
    Tk = [Tk zeros((k+2)*m,n); zeros(m,(k+2)*n)];
    Tk((k+1)*m+1:(k+2)*m,(k+1)*n+1:(k+2)*n) = A;
    Tk((k+2)*m+1:(k+3)*m,(k+1)*n+1:(k+2)*n) = -B;
    k = k+1;
    rk = (k+1)*n-rank(Tk);
  end
  if (k > min(m,n)) break; end;
  epsilonhat(i) = k;
  rho(i) = rk-((k+1)*p-q);
  M = computeHi(i,n,p,q,X,rho,epsilonhat);
  [U,S,V] = svd(M);
  if (M ~= []) sigma1 = S(1,1); else sigma1 = 0; end
  rankM = sum(diag(S) >  max(size(M))*sigma1*eps);
  N = null(Tk);
  N1 = U(:,1:rankM);
  if (N1 ~= []) P = N-N1*N1'*N; else P = N; end
  [U,S,V] = svd(P);
  if (P ~= []) sigma1 = S(1,1); else sigma1 = 0; end
  rankP = sum(diag(S) >  max(size(P))*sigma1*eps);
  N2 = U(:,1:rankP);
  if (i > 1)
    X = [X; zeros((epsilonhat(i)-epsilonhat(i-1))*n,p)];
  end
  X = [X N2];
  epsilon(p+1:p+rho(i)) = k*ones(1,rho(i));
  p = p+rho(i);
  q = q+rho(i)*epsilonhat(i);
  i = i+1;
end

u = i-1;
epsilon = epsilon(1:p);
epsilonhat = epsilonhat(1:u); rho = rho(1:u);