% program tridiagonal_solver

  disp(' This program finds the solution and inverse of a tridiagonal matrix and ')
  disp(' the determinant of this matrix ')
  disp(' The form of the matrix is ')
  disp('   d1   u1   0    0   ')
  disp('   l1   d2   u2   0   ')
  disp('   0    l2   d3   u3  ')
  disp('   0     0   l3   d4  ')
  disp(' ')
  disp(' The file low must contain the lower diagonal entries. ')
  disp(' The file diagonal must contain the diagonal entries. ')
  disp(' The file up must contain the upper diagonal entries. ')
  disp(' The file rhs must contain the right hand side of the equation. ')

% allocate the memory
% l is the lower diagonal, d the diagonal, u the upper diagonal and w
% is for row switches
% row_switch locates the rows of the matrix to switch
% b is the right hand side (Tx = b)
% inverse will be the inverse of the tridiagonal

% input the data

pause()

   load low     % lower diagonal
   load diagonal   % diagonal
   load up         % upper diagonal
   
   load rhs     % right hand side of matrix equation b  (Tx = b)

  l = low;
  d = diagonal;
  n = length(diagonal);
  u(1:n-1) = up;
  u(n) = 0;
  b = rhs;
  
  n = length(b);
  
  w(1:n) = 0;           % stores the upper diagonal row switch
  rswitch(1:n) = 0;     % indicates a row switch was done
  invmat(1:n,1:n) = 0;   % the inverse of the tridiagonal
  
  % detsign is for recording row switches for the determinant
  % dmn gives the determinant of T the tridiagonal matrix
  
  dmn = 1;
  detsign = 1;
  
  % put the tridiagonal matrix in upper triagonal form

    for j = 1:n-1

      if (abs(d(j)) < abs(l(j)))  % do a row switch

        [d(j),l(j)] = swap(d(j),l(j));
        [u(j),d(j+1)] = swap(u(j),d(j+1));
        [w(j),u(j+1)] = swap(w(j),u(j+1));
        
        rswitch(j) = 1;
        detsign = -detsign;

      end
      if (abs(d(j)) < 10^(-12)) % a row switch does not help

        disp(' Your matrix is not invertible ')
        disp(' It''s determinant is zero ')
        disp(' bye ')
        return

      end

% if everything is ok do the gaussian elimination
 
      l(j) = -l(j)/d(j);
      d(j+1) = d(j+1) + l(j)*u(j);
      u(j+1) = u(j+1) + l(j)*w(j);

    end
   
    for j = 1:n    
      
      if (abs(d(j)) < 10^(-12)) % a row switch does not help

        disp(' Your matrix is not invertible ')
        disp(' It''s determinant is zero ')
        disp(' bye ')
        return

      end
      
      dmn = d(j)*dmn;
        
    end
    
    dmn = detsign*dmn;
    
    disp(' ')
    disp([' The determinant of your tridiagonal is ', num2str(dmn)])
    disp(' ')

% This row reduces the right hand side of a matrix for
% a tridiagonal system

    for j = 1:n-1

      if (rswitch(j) == 1)   % do this if a row switch was done
                                % during gaussian elimination

        [b(j),b(j+1)] = swap(b(j),b(j+1));

      end

      b(j+1) = b(j+1) + l(j)*b(j);  % do the gaussian elimination

    end 
    
    x(n) = b(n)/d(n);  % find the last entry
    x(n-1) = (b(n-1) - x(n)*u(n-1))/d(n-1); % find the next to last entry                                                              
                                                                
% find the remaining entries

    for i = n-2:-1:1

      x(i) = (b(i) - x(i+2)*w(i) - x(i+1)*u(i))/d(i);
      
    end
    disp(' The solution for your ')
    disp(b)
    disp(' is ')
    disp(x)
    disp(' ')
                
% 
% This determines the inverse of the tridiagonal matrix

for j = 1:n
    
    b(1:n) = 0;
    b(j) = 1;
    
    for k = 1:n-1
        
      if (rswitch(k) == 1)   % do this if a row switch was done
                                % during gaussian elimination

        [b(k),b(k+1)] = swap(b(k),b(k+1));

      end

      b(k+1) = b(k+1) + l(k)*b(k);  % do the gaussian elimination
      
    end

    invmat(n,j) = b(n)/d(n);  % find the last entry
    invmat(n-1,j) = (b(n-1) - invmat(n,j)*u(n-1))/d(n-1); 
    % find the next to last entry
    
% find the remaining entries

    for i = n-2:-1:1

      invmat(i,j) = (b(i) - invmat(i+2,j)*w(i) - ...
                     invmat(i+1,j)*u(i))/d(i);

    end 
    
end

disp(' The inverse of your tridiagonal matrix is ')
disp(invmat)

