From 37feb7c484462fc15a31044cffd7ce2ddaac807a Mon Sep 17 00:00:00 2001 From: EdwardEisenhauer Date: Sun, 14 Mar 2021 18:40:53 +0100 Subject: [PATCH] Refactored comments --- .../Alg11.m | 5 ++- .../Alg1_outer_product_gaussian_elimination.m | 32 ++++++++------ ...ssian_elimination_with_complete_pivoting.m | 44 ++++++++++--------- .../Alg3_forward_substitution.m | 37 ++++++++++------ .../Alg4_back_substitution.m | 17 ++++--- .../Alg5_gauss_jordan_elimination.m | 4 +- .../Alg6_RREF.m | 42 +++++++++--------- .../README.md | 8 ++++ 8 files changed, 112 insertions(+), 77 deletions(-) create mode 100644 Direct Methods for Solving Linear Systems/README.md diff --git a/Direct Methods for Solving Linear Systems/Alg11.m b/Direct Methods for Solving Linear Systems/Alg11.m index 3fd774b..b33c7c4 100644 --- a/Direct Methods for Solving Linear Systems/Alg11.m +++ b/Direct Methods for Solving Linear Systems/Alg11.m @@ -1,5 +1,5 @@ function [Q R] = Alg11(A) -% Algorithm 11: QR factorization via Housholder algorithm. +% Algorithm 11: QR factorization via Householder algorithm. [m, n] = size(A); @@ -17,5 +17,8 @@ end R = triu(A) +H = tril(A, -1) + + end 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 31648c1..babce75 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,21 +1,25 @@ function A = Alg1_outer_product_gaussian_elimination(A) -% Algorithm 1: Outer Product Gaussian Elimination (Golub, Loan, 3.2.1) +% Algorithm 1: Outer Product Gaussian Elimination % Performs a gaussian eliminaion on a square matrix A. - [m, n] = size(A); +[m, n] = size(A); - if m ~= n - error('Matrix is not squared!') - end - -% if det(A) == 0 -% error('Matrix is not nonsingular!') -% end - - 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); +if m ~= n + error('Matrix is not square!') +end + +for k = 1:n-1 + if det(A(1:k, 1:k)) < eps + error('Matrix is not nonsingular!') end +end + +% The following algorithm is based on the Algrotihm 3.2.1 from [2]. + +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 end \ No newline at end of file diff --git a/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m b/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m index 7c01c8e..43f9cfe 100644 --- a/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m +++ b/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m @@ -1,43 +1,47 @@ function [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A) % Algorithm 2: Gaussian Elimination with Complete Pivoting. -% [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A) +% [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A) +% computes the complete pivoting factorization PAQ = LU. -[n, m] = size(A); +[m, n] = size(A); -if n ~= m - error('Matrix is not squared!') +if m ~= n + error('Matrix is not square!') end - -% if det(A) == 0 -% error('Matrix is not nonsingular!') -% end -p = 1:n; -q = 1:n; +% p and q are permutation vectors – respectively rows and columns +p = 1:m; +q = 1:m; -for k = 1 : n-1 - i = k:n; - j = k:n; +% The following algorithm is based on the Algrotihm 3.4.2 from [2]. + +for k = 1 : m-1 + i = k:m; + j = k:m; + % Find the maximum entry to be the next pivot [max_val, rows_of_max_in_col] = max(abs(A(i, j))); [~, max_col] = max(max_val); max_row = rows_of_max_in_col(max_col); - % Assign value of mu and lambda in respect to the main A matrix + % Assign value of mu and lambda in respect to the main matrix A [mi, lm] = deal(max_row+k-1, max_col+k-1); - A([k mi], 1:n) = deal(A([mi k], 1:n)); - A(1:n, [k lm]) = deal(A(1:n, [lm k])); + % Interchange the rows and columns of matrix A... + A([k mi], 1:m) = deal(A([mi k], 1:m)); + A(1:m, [k lm]) = deal(A(1:m, [lm k])); + % ...and respective permutation vectors entries. p([k, mi]) = p([mi, k]); q([k, lm]) = q([lm, k]); - % Perform Gaussian elimination with the greatest pivot if A(k, k) ~= 0 - rows = k+1 : n; + 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 end -I = eye(n); +I = eye(m); U = triu(A); L = tril(A, -1) + I; P = I(p, :); -Q = I(:, q); \ No newline at end of file +Q = I(:, q); + +end \ No newline at end of file diff --git a/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m b/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m index a0618c0..b793584 100644 --- a/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m +++ b/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m @@ -1,18 +1,29 @@ function b = Alg3_forward_substitution(L, b) -% Algorithm 3: Forward Substitution (Golub, Loan, Alg. 3.1.1) +% Algorithm 3: Forward Substitution +% b = Alg3_forward_substitution(L, b) overwrites b with the solution to +% Lx = b. - [m, n] = size(L); - if m ~= n - error('Matrix is not squar!') - end +[m, n] = size(L); + +if m ~= n + error('Matrix is not square!') +end + +if length(b) ~= m + error('Vector b has wrong length!') +end + +if det(L) < eps + error("Matrix is not nonsingular!") +end + +% The following algorithm is based on the Algrotihm 3.1.1 from [2]. + +% b(m, :) so that matrices are also accepted +b(1, :) = b(1, :)/L(1,1); +for i = 2:m + b(i, :) = (b(i, :) - L(i, 1:i-1)*b(1:i-1, :))/L(i, i); +end - if length(b) ~= m - error('Vector b has wrong length!') - end - - b(1, :) = b(1, :)/L(1,1); - for i = 2:m - b(i, :) = (b(i, :) - L(i, 1:i-1)*b(1:i-1, :))/L(i, i); - end end diff --git a/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m b/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m index 9f25ebf..ca4185e 100644 --- a/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m +++ b/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m @@ -1,11 +1,12 @@ function b = Alg4_back_substitution(U,b) -% Argorithm 4: Back Substitution (Golub, Loan, Alg. 3.1.2) -% Returns vetor b with solution to he Ux = b. +% Argorithm 4: Back Substitution +% b = Alg4_back_substitution(U,b) returns vetor b with solution to the +% Ux = b. [m, n] = size(U); if U ~= triu(U) - error('Matrix is not upper triangular!') + error('Matrix U is not upper triangular!') end if m ~= n @@ -16,12 +17,14 @@ if length(b) ~= m error('Vector b has wrong length!') end -% if det(U) < 0.001 -% error('Matrix is not nonsingular!') -% end +if det(U) < eps + error('Matrix is not nonsingular!') +end + + +% The following algorithm is based on the Algrotihm 3.1.2 from [2]. % b(m, :) so that matrices are also accepted - b(m, :) = b(m, :)/U(m, m); for i = m-1:-1:1 b(i, :) = (b(i, :) - U(i, i+1 : m)*b(i+1 : m, :))/U(i, i); diff --git a/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m b/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m index f510427..28f598d 100644 --- a/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m +++ b/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m @@ -1,6 +1,7 @@ function A = Alg5_gauss_jordan_elimination(A) % Algorithm 5: Gauss-Jordan Elimination -% Argument A is an augmented matrix +% A = Alg5_gauss_jordan_elimination(A) performs Gauss-Jordan elimination +% on an augmented matrix A. [m, n] = size(A); @@ -16,4 +17,5 @@ for k = 1 : m end end end + end \ No newline at end of file diff --git a/Direct Methods for Solving Linear Systems/Alg6_RREF.m b/Direct Methods for Solving Linear Systems/Alg6_RREF.m index 5d055db..c57bb83 100644 --- a/Direct Methods for Solving Linear Systems/Alg6_RREF.m +++ b/Direct Methods for Solving Linear Systems/Alg6_RREF.m @@ -2,38 +2,38 @@ function A = Alg6_RREF(A) % Algorithm 6: Reduced Row Echelon Form (RREF) % A = Alg6_RREF(A) returns RREF of matrix A. -[M, N] = size(A); +[m, n] = size(A); -n = 0; +j = 0; -for m = 1 : M - n = n + 1; - if n > N +for k = 1 : m + j = j + 1; + if j > n break end % We want the left-most coefficient to be 1 (pivot) - row = A(m, :); - if row(m) == 0 - n = n + 1; + row = A(k, :); + if row(k) == 0 + j = j + 1; end - row = row/row(n); - A(m, :) = row; + row = row/row(j); + A(k, :) = row; - for i = 1 : M - if i ~= m - A(i, :) = A(i, :)-(A(i, n))*row; + for i = 1 : m + if i ~= k + A(i, :) = A(i, :)-(A(i, j))*row; end end - for i = m + 1 : M - A(i:end, m+1:end); % Partial matrix (in which we are looking for non-zero pivots) - A(i:end, m+1); % Left-most column - if ~any(A(i:end, m+1)) % If the left-most column has only zeros check the next one - m = m + 1; + for i = k + 1 : m + A(i:end, k+1:end); % Partial matrix (in which we are looking for non-zero pivots) + A(i:end, k+1); % Left-most column + if ~any(A(i:end, k+1)) % If the left-most column has only zeros check the next one + k = k + 1; end - A(i:end, m+1:end); - if A(i, m+1) == 0 - non_zero_row = find(A(i:end,m+1), 1); + A(i:end, k+1:end); + if A(i, k+1) == 0 + non_zero_row = find(A(i:end,k+1), 1); if isempty(non_zero_row) continue end diff --git a/Direct Methods for Solving Linear Systems/README.md b/Direct Methods for Solving Linear Systems/README.md new file mode 100644 index 0000000..6119207 --- /dev/null +++ b/Direct Methods for Solving Linear Systems/README.md @@ -0,0 +1,8 @@ +# Direct Methods for Solving Linear Systems + +## Bibliography + +References in the code comments are represented by `[number]`: + + - 1 + - 2 – Golub, van Loan, Matrix Computations, 3rd edition \ No newline at end of file