home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l450 / 1.ddi / MATPAK.1 / arc / BACKSUB.M < prev    next >
Encoding:
Text File  |  1989-09-07  |  1.1 KB  |  48 lines
  1. function x = backsub(LR, b)
  2. %BACKSUB Backsubstitution to solve triangular systems of equations.
  3. %    X = BACKSUB(L,B) performs backsubstitution to solve
  4. %    the equation L*X = B, where L is lower triangular.
  5. %    BACKSUB is provided for pedagogical purposes; in most cases
  6. %    it is faster to use X = L\B.
  7.  
  8. %    Copyright (C) 1988 the MathWorks, Inc.
  9.  
  10. % Check if any diagonal elements are essentially 0.
  11. if any(abs(diag(LR)) < eps*norm(LR))
  12.     error('Matrix is singular.')
  13. end
  14. [m,n] = size(LR);
  15. if m ~= n
  16.     error('Matrix must be square.')
  17. end
  18.  
  19. % Find out if LR is L or R
  20. lower = 0;
  21. if norm(LR - tril(LR),1) == 0
  22.     lower = 1;
  23. elseif norm(LR - triu(LR),1) == 0
  24.     lower = 0;
  25. else
  26.     error('Matrix must be triangular')
  27. end
  28.  
  29. % Reorder rows of b and LR if system is upper triangular
  30. if ~lower
  31.     b = b(n:-1:1);
  32.     LR = LR(n:-1:1,n:-1:1);
  33. end
  34.  
  35. % initialize x to force it to be a column vector.
  36. x = zeros(n,1);
  37.  
  38. % perform backsubstitution
  39. x(1) = b(1) / LR(1,1);
  40. for i = 2:n
  41.     x(i) = (b(i) - LR(i,1:i-1) * x(1:i-1)) / LR(i,i);
  42. end
  43.  
  44. % reorder answer if system is upper triangular
  45. if ~lower
  46.     x = x(n:-1:1);
  47. end
  48.