Added Algorithm 8 and refactored code

.gitignore

@@ -1 +1,2 @@
+# MATLAB
 *.asv

Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m

@@ -1,20 +1,21 @@
 function A = Alg1_outer_product_gaussian_elimination(A)
-% ADDME Algorithm 1: Outer Product Gaussian Elimination (Golub, Loan, 3.2.1)
+% Algorithm 1: Outer Product Gaussian Elimination (Golub, Loan, 3.2.1)
-% Matrix A has to be a square singular matrix.
+% 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!')
+    if m ~= n
-end
+        error('Matrix is not squared!')
+    end
     
 %  if det(A) == 0
 %     error('Matrix is not nonsingular!')
 %  end
  
- for k = 1 : m-1
+    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
  
 end

Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m

@@ -1,27 +1,27 @@
 function [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A)
-% ADDME Algorithm 2: Gaussian Elimination with Complete Pivoting.
+% Algorithm 2: Gaussian Elimination with Complete Pivoting.
 % [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A) 
 
 [n, m] = size(A);
+
 if n ~= m
     error('Matrix is not squared!')
 end
     
-if det(A) == 0
+% if det(A) == 0
-    error('Matrix is not nonsingular!')
+%     error('Matrix is not nonsingular!')
-end
+% end
 
 p = 1:n;
 q = 1:n;
 
-% for k = 1 : n-1
 for k = 1 : n-1
     i = k:n;
     j = k:n;
     [max_val, rows_of_max_in_col] = max(abs(A(i, j)));
-    [max_val, max_col] = max(max_val);
+    [~, max_col] = max(max_val);
     max_row = rows_of_max_in_col(max_col);
-    % Assing value of mi and lambda in respect to the main A matrix
+    % Assign value of mu and lambda in respect to the main A matrix
     [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]));
@@ -36,8 +36,8 @@ for k = 1 : n-1
     end
 end
 
-U = triu(A);
-L = tril(A, -1) + eye(n);
 I = eye(n);
+U = triu(A);
+L = tril(A, -1) + I;
 P = I(p, :);
 Q = I(:, q);

Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m

@@ -1,17 +1,18 @@
-% Algorithm 3: Forward Substitution (Alg. 3.1.1)
 function b = Alg3_forward_substitution(L, b)
+% Algorithm 3: Forward Substitution (Golub, Loan, Alg. 3.1.1)
 
-[n, m] = size(L);
+    [m, n] = size(L);
-if n ~= m
+    if m ~= n
-    error('Matrix is not squared!')
+        error('Matrix is not squared!')
-end
+    end
-
+
-if length(b) ~= n
+    if length(b) ~= m
-    error('Vector b has wrong length!')
+        error('Vector b has wrong length!')
-end
+    end
-
+
-b(1, :) = b(1, :)/L(1,1);
+    b(1, :) = b(1, :)/L(1,1);
-for i = 2:n
+    for i = 2:m
-    b(i, :) = (b(i, :) - L(i, 1:i-1)*b(1:i-1, :))/L(i, i);
+        b(i, :) = (b(i, :) - L(i, 1:i-1)*b(1:i-1, :))/L(i, i);
+    end
 end
 

Direct Methods for Solving Linear Systems/Alg4_back_substitution.m

@@ -1,5 +1,5 @@
 function b = Alg4_back_substitution(U,b)
-% ADDME Argorithm 4: Back Substitution (Golub, Loan, Alg. 3.1.2)
+% Argorithm 4: Back Substitution (Golub, Loan, Alg. 3.1.2)
 % Returns vetor b with solution to he Ux = b.
 
 [m, n] = size(U);

Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m

@@ -2,19 +2,17 @@ function A = Alg5_gauss_jordan_elimination(A)
 % Algorithm 5: Gauss-Jordan Elimination
 % Argument A is an augmented matrix
 
-% M – rows, N – columns
+[m, n] = size(A);
 
-[M, N] = size(A);
+for k = 1 : m
-
-for m = 1 : M
     
-    row = A(m, :);
+    row = A(k, :);
-    row = row/row(m);
+    row = row/row(k);
-    A(m, :) = row;
+    A(k, :) = row;
     
-    for n = 1 : M
+    for i = 1 : m
-        if n ~= m
+        if i ~= k
-            A(n, :) = A(n, :)-(A(n, m))*row;
+            A(i, :) = A(i, :)-(A(i, k))*row;
         end
     end
 end

Direct Methods for Solving Linear Systems/Alg6_RREF.m

@@ -1,9 +1,7 @@
 function A = Alg6_RREF(A)
-% ADDME Algorithm 6: Reduced Row Echelon Form (RREF)
+% Algorithm 6: Reduced Row Echelon Form (RREF)
 % A = Alg6_RREF(A) returns RREF of matrix A. 
 
-% M – rows, N – columns
-
 [M, N] = size(A);
 
 n = 0;

Direct Methods for Solving Linear Systems/Alg7.m

@@ -1,7 +1,8 @@
 function [L, U] = Alg7(A)
-% ADDME LU Factorization without pivoting
+% LU Factorization without pivoting
-% [L, U] = Alg7(A) decomposes matrix A into U – upper triangular matrix and
+% [L, U] = Alg7(A) decomposes matrix A using Crout's algorithm into
-% L – lower unit triangular matrix such, that A = LU.
+% U – upper triangular matrix and L – unit lower triangular matrix
+% such that A = LU.
 
 [m, n] = size(A);
 
@@ -14,21 +15,21 @@ function [L, U] = Alg7(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)
+L = zeros(m, n);
-U = eye(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);
+        U(k, j) = A(k, j) - dot(L(k, 1:k-1), U(1:k-1, j));
     end
+    
+    for i = k:n
+        L(i, k) = (A(i, k) - dot(L(i,1:k-1), U(1:k-1, k))) / U(k,k);
+    end
+    
+    
 end
 
 end

Direct Methods for Solving Linear Systems/Alg8.m

@@ -0,0 +1,54 @@
+function [L, U, P] = Alg8(A)
+% LU Factorization with partial pivoting
+% [L, U, P] = Alg8(A) decomposes matrix A using Crout's algorithm into
+% U – upper triangular matrix and L – unit lower triangular matrix
+% using partial pivoting, interchanging A's rows, such that AP = LU.
+
+[m, n] = size(A);
+
+L = zeros(m, n);
+U = eye(m, n);
+P = eye(m, n);
+
+for k = 1:m
+    if A(k, k) == 0
+        non_zero_col = find(A(k, :), 1);
+        A(:, [k non_zero_col]) = deal(A(:, [non_zero_col k]));
+        P(:, [k non_zero_col]) = deal(P(:, [non_zero_col k]));
+    end
+end
+
+% TODO: Permutation A*P does not work
+
+
+% 
+% 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;
+%    end
+%    A(i:end, m+1:end);
+%    if A(i, m+1) == 0
+%        non_zero_row = find(A(i:end,m+1), 1);
+%        if isempty(non_zero_row)
+%            continue
+%        end
+%        A([i, i+non_zero_row-1], :) = deal(A([i+non_zero_row-1, i], :));
+%    end
+% end
+
+for k = 1:n
+    
+    for j = k : n
+        U(k, j) = A(k, j) - dot(L(k, 1:k-1), U(1:k-1, j));
+    end
+    
+    for i = k:n
+        L(i, k) = (A(i, k) - dot(L(i,1:k-1), U(1:k-1, k))) / U(k,k);
+    end
+    
+    
+end
+
+end