%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code demonstrates how to call symtrisolve.m for
% solving symmetric tridiagonal linear systems Tx = y.
% The symtrisolve.m code computes the LBL^T factorization of a
% symmetric tridiagonal matrix using the modified Bunch-Marcia
% pivoting strategy.
%
% Author: Roummel F. Marcia
% Date: July 10, 2014
%
% Based on the paper, "A simplified pivoting strategy for
% symmetric tridiagonal matrices", by James R. Bunch and
% Roummel F. Marcia, Numer. Linear Algebra Appl., vol. 13,
% p. 865-867 (2006).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
format compact
% Matrix size
n = 10000;
% Create random symmetric tridiagonal matrix
% TD is a vector corresponding to the diagonal entries
% TS is a vector corresponding to the sub- and super-diagonal entries
TD = randn(n,1);
TS = randn(n-1,1);
% Create random solution x
x = randn(n,1);
% Create RHS y
y = zeros(n,1);
y(1) = TD(1)*x(1) + TS(1)*x(2);
for j = 2:n-1
y(j) = TS(j-1)*x(j-1) + TD(j)*x(j) + TS(j)*x(j+1);
end
y(n) = TS(n-1)*x(n-1) + TD(n)*x(n);
% Solve using symmetric tridiagonal solver
x0 = symtrisolve(TD, TS, y);
% Verify that the solution is accurate by computing the
% relative residual norm
y0 = zeros(n,1);
y0(1) = TD(1)*x0(1) + TS(1)*x0(2);
for j = 2:n-1
y0(j) = TS(j-1)*x0(j-1) + TD(j)*x0(j) + TS(j)*x0(j+1);
end
y0(n) = TS(n-1)*x0(n-1) + TD(n)*x0(n);
fprintf('\nProblem size = %12d\n', n)
fprintf('Relative residual norm = %12.4e\n\n', norm(y-y0)/norm(y))