lu_tridiag.m

function [ dU, odL, odU ] = lu_tridiag( vector1, vector2, vector3 )
%Function makes a 3x3 matrix with three 1x3 vectors, then conducts an L U
%decomposition. Function also pivots matrix 
%   dU:  upper matrix diagonal
%   odU: upper matrix diagonals, order 4, 8, 7
%   odL: diagonal of lower matrix, order: 2, 6, 3

matrix = [ vector1; vector2; vector3];
pivoted = [];

zeros = matrix == 0;

row_one_zeros = find(zeros(:,1) == 1); % gives row # where zeros are in first column
row_two_zeros = find(zeros(:,2) == 1); % gives row # where zeros are in second column

if isempty(row_one_zeros)
    
    pivoted = matrix;
    
end
   
if isempty(row_two_zeros) == 0 && isempty(row_one_zeros) == 0
   
if sum(row_one_zeros == row_two_zeros) > 0  && length(row_one_zeros) > length(row_two_zeros) % begin pivoting sequence for a matrix with two zeros in a row
    
    % Next see if there is a row with zeros in both first and second column

    location = row_one_zeros == row_two_zeros; % logical vector where two vectors are equal
    third_row_number = sum(row_one_zeros' * location); % row number where two zeros exist
    pivoted(3, :) = matrix(third_row_number, :); % assigns row with two zeros to row three

    % Next look to see which rows have a zero only in the first row
    
    second_row_number = find(row_one_zeros ~= third_row_number); % finds whick row is not already used for row tree
    pivoted(2, :) = matrix(row_one_zeros(second_row_number),:); % assigns this row to row 2

    % Next assign last row to row one
    row_num = 1:3; 
    first_row_number = sum(row_num - row_num .* ismember(1:3,row_one_zeros')); % finds first row 
    pivoted(1,:) = matrix(first_row_number,:); % assigns top row to pivoted
    
end

end

if length(row_one_zeros) == 1 && matrix(1) == 0  % begin pivoting sequence for only zero in first row
    
    pivoted(1,:) = matrix(2,:);
    pivoted(2,:) = matrix(1,:);
    pivoted(3,:) = matrix(3,:);
        
end

if length(row_one_zeros) == 1 && matrix(1) ~= 0 % pivots if only one zero in first column
    
    if matrix(2) == 0
        pivoted(1,:) = matrix(1,:);
        pivoted(2,:) = matrix(3,:);
        pivoted(3,:) = matrix(2,:);
    end
    
    if matrix(3) == 0
        
        pivoted = matrix;
        
    end
    
end

if isempty(pivoted) && length(row_one_zeros) > 1  % brute force method for matrix with many zeros
    
    row_num = 1:3;
    
    first_row_number = sum(row_num - row_num .* ismember(1:3,row_one_zeros'));
    
    if first_row_number == 1
        
        pivoted = matrix;
        
    else 
        
    if length(row_two_zeros) > 0 
        
    second_row_number = sum(row_num - row_num .* ismember(1:3,row_two_zeros'));
    
    end
    
    if second_row_number == first_row_number 
        
        second_row_number = 2;
        
    end
    
    third_row_number = sum(1:3) - first_row_number - second_row_number;
    pivoted(1,:) = matrix(first_row_number,:);
    pivoted(2,:) = matrix(second_row_number,:);
    pivoted(3,:) = matrix(third_row_number,:);
    
    
    end
end

vector1 = pivoted(1,:);
vector2 = pivoted(2,:);
vector3 = pivoted(3,:);

% begins the l u decomposition

   l2 = vector2(1)/vector1(1);
   row2 = vector2 - l2*vector1; % row 2 of upper matrix
   l3 = vector3(1)/vector1(1);
   row3 = vector3 - l3*vector1;
   l4 = row3(2)/row2(2);
   row3 = row3 - l4*row2; % row 3 of upper matrix
   pivoted = [vector1; row2; row3]; % upper matrix
   dU = diag(pivoted)'; % diagonal of upper matrix
   L = [ 1 0 0; l2 1 0; l3 l4 1];
   verify = L*pivoted; 
   
   if verify ~= matrix % makes sure L U decomposition is correct, else displays an error message
       fprintf('An error has occured \n')
   end
   
   odL = [ l2 l4 l3 ]; % off diagonal of lower matrix, order: 2, 6, 3
   odU = [ diag(pivoted, 1)', diag(pivoted,2)]; % off diagonal of upper matrix, order: 4, 8, 7
   
end