Product Triangular Systems With Shift

Matlab Codes

 TriProd Implicit Method, Triangular Case TriProdExp Explicit Method, Triangular Case QTriProd Implicit Method, Quasi-Triangular Case QTriProdExp Explicit Method, Quasi-Triangular Case Test Compares implicit and explicit methods GenProb Generates problems of specified size ErrBound Computes relative residual

TriProd

  function x = TriProd(R,lambda,b);
%
% x = TriProd(R,lambda,b)
% Computes the solution to
%
%               (R1*...*Rp - lambda*I)x = b
%
% where R is a px1 cell array containing
% n-by-n upper triangular matrices R{1},..,R{p},
% b is a n-by-1 column vector and lambda is a real scalar.
%
% Uses the implicit method.

p = length(R);
n = length(R{1});

% Calculate the first xc, i.e., x(n)

Rprod = R{1}(n,n);
for i=2:p
Rprod = Rprod*R{i}(n,n);
end
xc = b(n)/(Rprod-lambda);

% Create the first WC vectors.

WC(:,p) = xc;
for j=p-1:-1:1
WC(:,j) = R{j+1}(n,n)*WC(:,j+1);
end
% WC(:,j) holds the jth WC vector

for k=1:n-1
% Compute x(n-k)
gamma = b(n-k);
sigProd = 1;
for j=1:p
gamma   = gamma - sigProd*(R{j}(n-k,n-k+1:n)*WC(:,j));
sigProd = sigProd*R{j}(n-k,n-k);
end
gamma = gamma/(sigProd - lambda);

% Update xc
xc = [gamma;xc];

% Update the WC vectors
if k<n-1
omega(p) = gamma;
for j=p-1:-1:1
omega(j) = R{j+1}(n-k,n-k)*omega(j+1)+R{j+1}(n-k,n-k+1:n)*WC(:,j+1);
end
WC = [omega;WC];
end
end
x = xc;

TriProdExp
  function x = TriProdExp(R,lambda,b);
%
% x = TriProdExp(R,lambda,b)
% Computes the solution to
%
%               (R1*...*Rp - lambda*I)x = b
%
% where R is a px1 cell array containing
% n-by-n upper triangular matrices R{1},..,R{p},
% b is a n-by-1 column vector and lambda is a real scalar.
%
% Uses the explicit method.

p = length(R);
n = length(R{1});

% Explicitly form the matrix of coefficients
A = R{1};
for j=2:p
A = A*R{j};
end
A = A - lambda*eye(n,n);

% Apply back-substitution
x = A\b;


QTriProd

  function x = QTriProd(R,lambda,b)
%
% x = QTriProd(R,lambda,b)
% Computes the solution to
%
%               (R1*...*Rp - lambda*I)x = b
%
% where R is a px1 cell array containing
% n-by-n quasi-upper triangular matrix R{1}
% and n-by-n upper triangular matrices R{2},..,R{p}.
% b is a n-by-1 column vector and lambda is a real scalar.
%
% Uses the implicit method.

p = length(R);
n = length(R{1});

if abs(R{1}(n,n-1)) > 0
m = n-1;
else
m = n;
end

% Calculate the first xc, i.e., x(m:n)

Rprod = R{1}(m:n,m:n);
for i=2:p
Rprod = Rprod*R{i}(m:n,m:n);
end
xc = (Rprod-lambda*eye(n-m+1,n-m+1))\b(m:n);

% Create the first WC vectors.

WC(:,p) = xc;
for j=p-1:-1:1
WC(:,j) = R{j+1}(m:n,m:n)*WC(:,j+1);
end
% WC(:,j) holds the jth WC vector

k=n-m+1;
while k<=n-1

% Have the bottom k components of x, i.e., x(n-k+1:n).
% Do we go after the next component x(n-k) or the next two
% components x(n-k-1) and x(n-k).

m = n-k;
if k<n-1
if abs(R{1}(n-k,n-k-1))>0
m = n-k-1;
end
end

% Compute x(m:n-k)

gamma = b(m:n-k);
sigProd = eye(length(gamma),length(gamma));
for j=1:p
gamma = gamma - sigProd*(R{j}(m:n-k,n-k+1:n)*WC(:,j));
sigProd = sigProd*R{j}(m:n-k,m:n-k);
end
gamma = (sigProd - lambda*eye(n-k-m+1,n-k-m+1))\gamma;

% Update xc
xc = [gamma;xc];

% Update the WC vectors
if k<n-1
omega(1:length(gamma),p) = gamma;
for j=p-1:-1:1
omega(:,j) = R{j+1}(m:n-k,m:n-k)*omega(:,j+1)+R{j+1}(m:n-k,n-k+1:n)*WC(:,j+1);
end
WC = [omega;WC];
end
k = k + length(gamma);
end
x = xc;


QTriProdExp

   function x = QTriProdExp(R,lambda,b)
%
% x = QTriProdExp(R,lambda,b)
% Computes the solution to
%
%               (R1*...*Rp - lambda*I)x = b
%
% where R is a px1 cell array containing
% n-by-n upper quasi-triangular matrix R{1}
% and n-by-n upper triangular matrices R{2},..,R{p}.
% b is a n-by-1 column vector and lambda is a real scalar.
%
% Uses the explicit method.

p = length(R);
n = length(R{1});

% Explicitly form the matrix of coefficients
A = R{p};
for j=p-1:-1:1
A = TriMult(R{j},A);
end
A = A - lambda*eye(n,n);

% Apply back-substitution

x = zeros(n,1);
i = n;
while (i>=1)
j=i;
if i>1
if abs(A(i-1,i))>0
j=i-1;
end
end

% We solve for x(j:i).
% Note, j=i-1 means that we have encountered a 2-by-2 bump.

x(j:i) = A(j:i,j:i)\b(j:i);
if j-1>=1
b(1:j-1) = b(1:j-1) - A(1:j-1,j:i)*x(j:i);
end
i = j-1;
end



Test

   % Tests Error Bounds for QTriProdExp and QTriProd
home
randn('seed',0)
n = 20:20:100;
p = [2 5 10 15];
smallnn = [.01 .0001 .000001 .00000001];
m = 100;
for k = 0:length(smallnn)
disp(' ')
if k==0
disp('Well-Conditioned')
else
disp(sprintf('smallnn = %8.2e',smallnn(k)))
end

Ratio12 = zeros(length(n),length(p));
Ratio21 = zeros(length(n),length(p));
for i = 1:length(n)
for j = 1:length(p)
for sample = 1:m
if k==0
[R,lambda,b] = GenProb(n(i),p(j));
else
[R,lambda,b] = GenProb(n(i),p(j),smallnn(k));
end
x1 = QTriProd(R,lambda,b);    E1 = ErrBound(R,lambda,b,x1);
x2 = QTriProdExp(R,lambda,b); E2 = ErrBound(R,lambda,b,x2);

Ratio12(i,j) = max(Ratio12(i,j),E1/E2);
Ratio21(i,j) = max(Ratio21(i,j),E2/E1);
end
end
end

disp(' ')
Ratio12 = Ratio12
disp(' ')

disp(' ')
Ratio21 = 1./Ratio21
disp(' ')

end


GenProb
  function [R,lambda,b] = GenProb(n,p,smallnn)
%
% n and p are positive integers.
% R is a p-by-1 cell array with R{1} is an n-ny-n
% random upper quasi-triangular matrix and R{2},...,R{p}
% is an n-by-n random upper triangular matrix.
% b is an n-by-1 random vector.
% If nargin == 2 then lambda is chosen so that
%
%      A = R{1}*...*R{n} - lambda*I
%
% is well-conditioned. Otherwise, it is chosen so that
% the (n,n) entry of A equals smallnn.

b = randn(n,1);
R = cell(p,1);
for j = 1:p
R{j} = triu(randn(n,n));
end

% Make R{1} quasi-triangular.
for i=1:2:n-2
R{1}(i+1,i) = 1;
end

% Determine lambda
if nargin==2
lambda = 3*n;
else
prodnn = 1;
for i=1:p
prodnn = prodnn*R{i}(n,n);
end
lambda = prodnn - smallnn;
end


ErrBound
  function e = ErrBound(R,lambda,b,x)
% Computes the relative residual for a computed solution to
%
%               (R1*...*Rp - lambda*I)x = b
%
% where R is a px1 cell array containing
% n-by-n quasi-upper triangular matrix R{1}
% and n-by-n upper triangular matrices R{2},..,R{p}.
% b is a n-by-1 column vector and lambda is a real scalar.

p = length(R);
n = length(R{1});

% Form z =  R{1}* ... *R{p}x
z = x;
%d = 1;

for j=p:-1:1
z = R{j}*z;

end
% Compute residual vector
res = b - (z - lambda*x);
e = (norm(res)/norm(x));