Shifted Kronecker Product Systems

Matlab Codes

KPShiftSolve

  function x = KPShiftSolve(A,n,b,lambda,alpha)
% Solves the shifted system
%   (alpha x Ap x ... x A2 x A1 - lambda*eye(N,N)) x = b
%
% INPUTS:
%       n : px1 vector of positive integers
%       A : px1 cell array, where A{k} is an n(k)-by-n(k) upper
%           quasi-tringular matrix, k=1:p
%       b : Nx1 column vector, where N = n(1)*...*n(p).
%  lambda : scalar
%   alpha : 1-by-1 or 2-by-2 matrix [optional argument]
%            If nargin == 4 then alpha is set to 1
%
% OUTPUTS:
%       x : Nx1 vector which solves the above shifted system

p = length(n); 
if nargin==4
   alpha = 1;
end
 
if length(alpha)==1
   
   A{p} = alpha*A{p};
   if p==1
      % if Hessenburg sys - call Hess solver
      if all(diag(A{1},-1))
        [v,U]=HessLU(A{1}-lambda*eye(n(1)));
        y=LBiDiSol(v,b);
        x=UTriSol(U,y);
      else  %not Hess - use quasi-triangular solver
        x = QuasiSolve(A{1},b,lambda);
      end
   else
      N = prod(n);
      x = zeros(N,1);
      mp = N/n(p);
      i = n(p);
      while (i>=1)
         if (i>1) & (abs(A{p}(i,i-1)) > 0 )
            %2-by-2 bump
            alpha = A{p}(i-1:i,i-1:i);
            idx = 1+(i-2)*mp:i*mp;
            x(idx) = KPShiftSolve(A,n(1:p-1),b(idx),lambda,alpha);
            
            z1 = KPvecprod(A,n(1:p-1),x(1+(i-2)*mp:(i-1)*mp));
            z2 = KPvecprod(A,n(1:p-1),x(1+(i-1)*mp:i*mp));
            for j=1:i-2
               jdx = 1+(j-1)*mp:j*mp;
               b(jdx) = b(jdx) - A{p}(j,i-1)*z1 - A{p}(j,i)*z2;
            end
            i = i-2;
         else
            %1-by-1 bump
            alpha = A{p}(i,i);
            idx = 1+(i-1)*mp:i*mp;

            x(idx) = KPShiftSolve(A,n(1:p-1),b(idx),lambda,alpha);
            
            z = KPvecprod(A,n(1:p-1),x(idx));
            for j=1:i-1
               jdx = 1+(j-1)*mp:j*mp;

               b(jdx) = b(jdx) - A{p}(j,i)*z;
            end
            i = i-1;
         end       
      end 
   end
else
   % alpha is a 2-by-2 matrix with complex conjugate eigenvalues
   [Q,T] = schur(alpha); 
   [Q,T] = rsf2csf(Q,T);

   % c = (Q kron I)'b
   m = length(b)/2;
   c = reshape(reshape(b,m,2)*conj(Q),2*m,1);
   
   % Solve (T kron A{p} kron ... kron A{1} - lambda I)z = c
   A{p+1} = T; n(p+1) = 2;
   z = KPShiftSolve(A,n,c,lambda);
   
   % x = (Q kron I)z
   x = reshape(reshape(z,m,2)*Q.',2*m,1);
   if norm(imag(x))<=100*eps
      x = real(x);
   end
end

 
HessLU

  function [v,U] = HessLU(A)
% [v,U] = HessLU(A)
%
% Computes LU factorization of upper Hessenburg, where U
% is upper triangular and L is lower unit triangular
%
% Inputs:
%   A : nxn upper Hessenburg matrix
%
% Outputs:
%   U : nxn upper triangular 
%   v : nx1 column vector with the property that 
%       L = I + diag(v(2:n),-1)

[n,n] = size(A);
v = zeros(n,1);
for k=1:n-1
   v(k+1) = A(k+1,k)/A(k,k);
   A(k+1,k:n) = A(k+1,k:n) - v(k+1)*A(k,k:n);
end
U = triu(A);

LBiDiSol

function x = LBiDiSol(l,b)
% Solves Lx=b for x when L lower bi-diagonal
%
% Inputs:
%   l : nx1 column vector representing a lower bi-diagonal
%       matrix L, L = I + diag(l(2:n),-1)
%   b : nx1 column vector
%
% Outputs:
%   x : nx1 column vector, solution to Lx=b

n = length(b); 
x = zeros(n,1);
x(1) = b(1);
for i=2:n
   x(i) = b(i) - l(i)*x(i-1);
end

UTriSol

  function x = UTriSol(U,b)
% Solves Ux=b for x when U upper triangular
%
% Inputs:
%   U : nxn upper triangular 
%   b : nx1 column vector
%
% Outputs: 
%   x : nx1 column vector, solution to Ux=b

n = length(b);
x = zeros(n,1);
for j=n:-1:2
   x(j) = b(j)/U(j,j);
   b(1:j-1) = b(1:j-1) - x(j)*U(1:j-1,j);
end
x(1) = b(1)/U(1,1);
   

QuasiSolve

  function x = QuasiSolve(T,b,lambda)
% Solves (T - lambda*I)x=b for x
%
% Inputs:
%  lambda : scalar
%       T : nxn upper quasi-triangular, (T - lambda*I nonsingular)
%       b : nx1 column vector
%
% note: if nargin == 2 then lambda is set to zero
%
% Outputs:
%   x : nx1 column vector, solution to (T - lambda*I)x=b

n = length(b);
x = zeros(n,1);
if nargin==2
   lambda = 0;
end

k = n;
while (k>=1)
   if k>1 & (abs(T(k,k-1)) > 0)
      % 2-by-2
      x(k-1:k) = (T(k-1:k,k-1:k) - lambda*eye(2,2))b(k-1:k);
      b(1:k-2) = b(1:k-2) - T(1:k-2,k-1:k)*x(k-1:k);
      k = k-2;
   else
      %1-by-1
      x(k) = b(k)/(T(k,k)-lambda);
      b(1:k-1) = b(1:k-1) - T(1:k-1,k)*x(k);
      k = k-1;
   end
end

function y = KPvecprod(A,n,x)
% Computes y = (A{p} kron ... kron A{1})*x
% in 2N*(n(1) + ... + n(p)) flops
%
% Inputs:
%   n : px1 vector of positive integers
%   A : px1 cell array, each A{i} is n(i)-by-n(i)
%
% Outputs:
%   x : Nx1 column vector, N=prod(n)

p = length(n);
N = prod(n);
m = N./n;
    
for i=1:p
   x = (A{i}*reshape(x,n(i),m(i))).';     
end
y = reshape(x,N,1);