Optimize partial pivoting LU factorization implementation

- Simplify row swapping using slice assignment
- Improve pivot selection with more concise findmax syntax
- Update test case with a smaller matrix to highlight numerical stability

Author: yangjunjie0320

src/partial-pivoting-lu-factorization.jl

@@ -27,21 +27,12 @@ function kernel(prob::PartialPivotingLUFactorizationProblem)
     p = collect(1:n)
 
     for k in 1:n-1
-        v, x = findmax(abs.(u[k:n, k]))
+        v, x = findmax(i->abs(u[i, k]), k:n)
         x += k - 1
 
         if x != k
-            # swap row k and row x of matrix u
-            for i in 1:n
-                u[k, i], u[x, i] = u[x, i], u[k, i]
-            end
-
-            # swap row k and row x of matrix l
-            for i in 1:n
-                l[k, i], l[x, i] = l[x, i], l[k, i]
-            end
-
-            # swap row k and row x of matrix p
+            u[k, :], u[x, :] = u[x, :], u[k, :]
+            l[k, :], l[x, :] = l[x, :], l[k, :]
             p[k], p[x] = p[x], p[k]
         end
 

src/qr-factorization.jl

@@ -0,0 +1,42 @@
+struct QRFactorizationSolution <: SolutionMixin
+    q::AbstractMatrix
+    r::AbstractMatrix
+    function QRFactorizationSolution(q, r)
+        sol = new(q, r)
+        @assert isortho(q) "matrix must be orthogonal"
+        @assert istriu(r) "matrix must be upper triangular"
+        return sol
+    end
+end
+
+struct QRFactorizationProblem <: ProblemMixin
+    a::AbstractMatrix
+end
+
+QRFactorization = QRFactorizationProblem
+
+function householder_reflection!(q, r)
+    m, n = size(r)
+    
+    if m == 1
+        return q, r
+    else
+        h = Matrix{Float64}(I, m, m)
+        for i in 1:m-1
+            h[i+1:m, i] = -2 * r[i+1:m, i] / r[i, i]
+        end
+        q = h * q
+        r = h * r
+        householder_reflection!(q, r)
+    end
+end
+
+function kernel(prob::QRFactorizationProblem)
+    r = copy(prob.a)
+    n = size(r, 1)
+    q = Matrix{Float64}(I, n, n)    
+
+
+end
+
+export PartialPivotingLUFactorization

test/partial-pivoting-lu-factorization.jl

@@ -1,7 +1,7 @@
 @testset "Partial pivoting LU factorization" begin
-    n = 4
-    tol = 1e-10
-    a = rand(n, n)
+    tol = 1e-6
+    a = [1e-8 1; 1 1]
+    n = size(a, 1)
 
     p = PartialPivotingLUFactorization(a, tol)
     s = kernel(p)