diff --git a/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m b/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m index 1e47f17..565a867 100644 --- a/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m +++ b/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m @@ -1,9 +1,9 @@ function A = Alg1_outer_product_gaussian_elimination(A) % ADDME Algorithm 1: Outer Product Gaussian Elimination (Golub, Loan, 3.2.1) +% Matrix A has to be a square singular matrix. - -[n, m] = size(A); -if n ~= m +[m, n] = size(A); +if m ~= n error('Matrix is not squared!') end @@ -11,8 +11,8 @@ end % error('Matrix is not nonsingular!') % end - for k = 1 : n-1 - rows = k + 1 : n; + for k = 1 : m-1 + rows = k + 1 : m; A(rows, k) = A(rows, k)/A(k, k); A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows); end diff --git a/Direct Methods for Solving Linear Systems/Alg7.m b/Direct Methods for Solving Linear Systems/Alg7.m index ecb620a..6e8e9ca 100644 --- a/Direct Methods for Solving Linear Systems/Alg7.m +++ b/Direct Methods for Solving Linear Systems/Alg7.m @@ -5,9 +5,30 @@ function [L, U] = Alg7(A) [m, n] = size(A); -A = Alg1_outer_product_gaussian_elimination(A); +% This solution is mathematically correct, however computationaly +% inefficient. +% +% A = Alg1_outer_product_gaussian_elimination(A); +% +% U = triu(A); +% L = tril(A, -1) + eye(m); -U = triu(A); -L = tril(A, -1) + eye(m); +% Instead, we should use the Crout's algorithm +% It returns an Unit Upper Triangular matrix and a Lower Traingular matrix +% (in oppose to the Gaussian Elimination). +A + +L = zeros(m, n) +U = eye(m, n) + +for k = 1:n + for i = k:n + L(i, k) = A(i, k) - dot(L(i,1:k-1), U(1:k-1, k)); + end + + for j = k : n + U(k, j) = (A(k, j) - dot(L(k, 1:k-1), U(1:k-1, j))) / L(k, k); + end +end end \ No newline at end of file