Added Algorithm 7 and refactored prevoius Algorithms

Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m

@@ -1,14 +1,15 @@
-% Algorithm 1: Outer Product Gaussian Elimination (Alg. 3.2.1)
 function A = Alg1_outer_product_gaussian_elimination(A)
+% ADDME Algorithm 1: Outer Product Gaussian Elimination (Golub, Loan, 3.2.1)
+
 
 [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
  
  for k = 1 : n-1
      rows = k + 1 : n;

Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m

@@ -1,4 +1,6 @@
 function [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A)
+% ADDME 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

Direct Methods for Solving Linear Systems/Alg4_back_substitution.m

@@ -1,24 +1,30 @@
-% Argorithm 4: Back Substitution (Alg. 3.1.2)
 function b = Alg4_back_substitution(U,b)
+% ADDME Argorithm 4: Back Substitution (Golub, Loan, Alg. 3.1.2)
+% Returns vetor b with solution to he Ux = b.
 
-[n, m] = size(U);
+[m, n] = size(U);
-if n ~= m
+
+if U ~= triu(U)
+    error('Matrix is not upper triangular!')
+end
+
+if m ~= n
     error('Matrix is not squared!')
 end
 
-if length(b) ~= n
+if length(b) ~= m
     error('Vector b has wrong length!')
 end
 
-if det(U) == 0
+% if det(U) < 0.001
-    error('Matrix is not nonsingular!')
+%     error('Matrix is not nonsingular!')
-end
+% end
 
-% b(n, :) so that matrices are also accepted
+% b(m, :) so that matrices are also accepted
 
-b(n, :) = b(n, :)/U(n, n);
+b(m, :) = b(m, :)/U(m, m);
-for i = n-1:-1:1
+for i = m-1:-1:1
-    b(i, :) = (b(i, :) - U(i, i+1 : n)*b(i+1 : n, :))/U(i, i);
+    b(i, :) = (b(i, :) - U(i, i+1 : m)*b(i+1 : m, :))/U(i, i);
 end
 
 end

Direct Methods for Solving Linear Systems/Alg6_RREF.m

@@ -1,5 +1,6 @@
 function A = Alg6_RREF(A)
-%  Algorithm 6: Reduced Row Echelon Form (RREF)
+% ADDME Algorithm 6: Reduced Row Echelon Form (RREF)
+% A = Alg6_RREF(A) returns RREF of matrix A. 
 
 % M – rows, N – columns
 
@@ -8,17 +9,15 @@ function A = Alg6_RREF(A)
 n = 0;
 
 for m = 1 : M
-    n = n + 1
+    n = n + 1;
     if n > N
         break
     end
-    A
     % We want the left-most coefficient to be 1 (pivot)
     row = A(m, :);
     if row(m) == 0
-        n = n + 1
+        n = n + 1;
     end
-    [m ,n]
     row = row/row(n);
     A(m, :) = row;
     
@@ -28,8 +27,6 @@ for m = 1 : M
         end
     end
     
-    A
-    
     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

Direct Methods for Solving Linear Systems/Alg7.m

@@ -0,0 +1,13 @@
+function [L, U] = Alg7(A)
+% ADDME LU Factorization without pivoting
+% [L, U] = Alg7(A) decomposes matrix A into U – upper triangular matrix and
+% L – lower unit triangular matrix such, that A = LU.
+
+[m, n] = size(A);
+
+A = Alg1_outer_product_gaussian_elimination(A);
+
+U = triu(A);
+L = tril(A, -1) + eye(m);
+
+end

Direct Methods for Solving Linear Systems/main.m

@@ -1,4 +1,5 @@
 clear all;
+clc;
 
 A = [2, -1, 0, 0;
     -1, 2, -1, 0;
@@ -9,8 +10,7 @@ b = [0;0;0;5];
 
 B = Alg1_outer_product_gaussian_elimination(A);
 U = triu(B);
-L = tril(B, -1);
+x = Alg4_back_substitution(U, b);
-x = back_substitution(U, b);
 
 %% Problem 1
 clear all;
@@ -44,9 +44,9 @@ b = [1;2;1];
 
 b = P*b;
 % Ly = b and Ux = y
-y = forward_substitution(L, b);
+y = Alg3_forward_substitution(L, b);
 
-x = Q*back_substitution(U, y)
+x = Q*Alg4_back_substitution(U, y)
 
 % L*U