From c597641f65e01d0e43d44407909990735c968a1a Mon Sep 17 00:00:00 2001
From: =?utf8?q?Nils=20Forss=C3=A9n?= <forssennils@gmail.com>
Date: Thu, 22 Sep 2022 15:04:00 +0200
Subject: [PATCH] update

---
 gaussSeidel.asv | 35 -----------------------------------
 lufact.asv      | 48 ++++++++++++++++++++++++++++++++++++++++++++++++
 lufact.m        |  9 ++++++++-
 3 files changed, 56 insertions(+), 36 deletions(-)
 delete mode 100644 gaussSeidel.asv
 create mode 100644 lufact.asv

diff --git a/gaussSeidel.asv b/gaussSeidel.asv
deleted file mode 100644
index c3502ea..0000000
--- a/gaussSeidel.asv
+++ /dev/null
@@ -1,35 +0,0 @@
-function [x,k] = gaussSeidel(A,b,x0,tol)
-%%
-%  Gauss-Seidel iterative method to approximate the solution of a
-%  linear system Ax=b up to a user defined tolerance
-%
-%  INPUT: 
-%    A   - n by n square, non-singular matrix
-%    b   - n by 1 right hand side vector
-%    x0  - n by 1 vector containing that initial guess for the iteration
-%    tol - user set tolerance for the stopping condition in the iteration 
-%
-%  OUTPUT:
-%    x - n by 1 vector containing the iterative solution
-%    k - number of iterations
-%%
-%  get the system size
-    n = length(A);
-    max_its = 250;
-    x = zeros(n, 1);
-    for k = 1:max_its
-        for i = 1:n
-            sigma = 0;
-            for j = 1:n
-                if j ~= i
-                    sigma = sigma + (A(i,j)*x(j));
-            x(i) = (b(i) - sigma)/a(i,i);
-        if norm(b - A*x, 2) <<
-
-%%
-%  Gauss-Seidel iteration which overwrites the current approximate solution
-%  with the new approximate solution (pseudocode available in the lecture
-%  notes on page 46)
-
-%
-end
\ No newline at end of file
diff --git a/lufact.asv b/lufact.asv
new file mode 100644
index 0000000..17f9c12
--- /dev/null
+++ b/lufact.asv
@@ -0,0 +1,48 @@
+function [L,U,P] = lufact(A)
+%%
+%  Computes the LU factorization of the matrix A. The factorization is done
+%  in-place and then the L and U matrices are extracted at the end for
+%  output. The factorization is PA = LU
+%
+%  INPUT: 
+%    A - n by n square, non-singular matrix
+%
+%  OUTPUT:
+%    L - lower triangular matrix with ones along the main diagonal
+%    U - upper triangular matrix
+%    P - permutation matrix for pivoting
+%%
+%  get the size of the system
+   n  = length(A);
+%%
+%  initialize the pivoting matrix
+   P = eye(n);
+%%
+%  Vectorized in-place LU factorization (with row pivoting) that keeps
+%  track of the total permutations by scrambleing the matrix P
+   for j = 1:n-1
+
+        [v i] = max(abs(A(:, j)));
+        A([k i],:) = A([i k],:);
+        P([k i],:) = P([i k],:);
+%%
+% Find the index of the largest pivot element in the current column
+
+
+%%
+% Swap the rows within the in-place array as well as the permutation matrix P
+    
+
+
+%%
+% Perform the in-place elimination and save the new column of L
+      i = j+1:n; % indices for the "active" matrix portion
+      A(i,j) = A(i,j)/A(j,j);
+      A(i,i) = A(i,i) - A(i,j)*A(j,i);
+      %A, return
+   end
+%%
+% Extract L and U from the in-place form
+   U = triu(A);
+   L = eye(n) + tril(A,-1);
+end
diff --git a/lufact.m b/lufact.m
index 68dcdce..5ca3886 100644
--- a/lufact.m
+++ b/lufact.m
@@ -21,11 +21,18 @@ function [L,U,P] = lufact(A)
 %  Vectorized in-place LU factorization (with row pivoting) that keeps
 %  track of the total permutations by scrambleing the matrix P
    for j = 1:n-1
+        A
+        [v i] = max(abs(A(j:n, j)));
+        A([j i+j-1],:) = A([i+j-1 j],:);
+        P([j i+j-1],:) = P([i+j-1 j],:);
+        A
 %%
 % Find the index of the largest pivot element in the current column
 
+
 %%
 % Swap the rows within the in-place array as well as the permutation matrix P
+    
 
 
 %%
@@ -33,7 +40,7 @@ function [L,U,P] = lufact(A)
       i = j+1:n; % indices for the "active" matrix portion
       A(i,j) = A(i,j)/A(j,j);
       A(i,i) = A(i,i) - A(i,j)*A(j,i);
-%      A, return
+      %A, return
    end
 %%
 % Extract L and U from the in-place form
-- 
2.30.2