demo_linsys.m

SFLOP simulates a floating point operation x op y where op = +, -, *, /

clear all;

format long e


disp(' ')
disp(' Create an n by n  singular upper trigiagonal matrix A    ')
disp(' and a consistent right hand side b.                      ')
disp(' In exact arithmetic A x = b has infinitely many solutions')
disp(' Try to solve the linear system                           ')
disp(' ')
 
n = input(' n = ');
disp(' ')

% Create a random upper triangular matrix
A = zeros(n,n);
for i = 1:n
    A(1:i,i) = rand(i,1);
end

% Set one diagonal entry equal to zero to make A singular
n2 = fix(n/2);
A(n2,n2) = 0;

% Create a consistent right hand side
xex = ones(n,1);
b   = A*xex;


% forward solve
x = b;
for i = n:-1:1
    if( A(i,i) == 0 & x(i) ~= 0 )
        disp([' The ', int2str(i),'th diagonal of A is zero,',...
              ' but x(',int2str(i),') = ', num2str(x(i)) ])     
        disp(' The system is not solvable ')
        break
    end
    x(i) = x(i) / A(i,i);
    x(1:i-1) = x(1:i-1) - A(1:i-1,i)*x(i);
end


disp(' ')
disp(' ')
disp(' ')
disp(' Create an n by n  singular trigiagonal matrix A      ')
disp(' Compute its LU decomposition and see if U is singular')
disp(' ')
 
n = input(' n = ');
disp(' ')

% Create a random upper triangular matrix
A = 100*rand(n,n);

% Set one column of A equal to a linear combination 
% of the previous columns to make A singular
n2 = fix(n/2);
A(:,n2) = A(:,1:n2-1)*rand(n2-1,1);

[L, U] = lu(A);

disp(' Diagonal entries of U:')
disp( diag(U) )