% 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);