From d6b6b159af3734853fd624a28c031f9d1f7aa641 Mon Sep 17 00:00:00 2001
From: EdwardEisenhauer <sergiusz.warga@pwr.edu.pl>
Date: Sat, 13 Mar 2021 19:02:40 +0100
Subject: [PATCH] Crout's algorithm works just fine

---
 .../Alg1_outer_product_gaussian_elimination.m | 10 +++----
 .../Alg7.m                                    | 27 ++++++++++++++++---
 2 files changed, 29 insertions(+), 8 deletions(-)

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