l2/z3

l2/3.jl

@@ -0,0 +1,43 @@
+#!/usr/bin/env julia
+
+# Jacek Poziemski 272389
+
+using LinearAlgebra
+
+include("hilb.jl")
+include("matcond.jl")
+
+struct Errors
+    gauss::Float64
+    inverse::Float64
+end
+
+function solve(A::Matrix{Float64}, n::Int)::Errors
+    x::Vector{Float64} = ones(Float64, n)
+    b::Vector{Float64} = A * x
+    gauss::Float64 = norm((A \ b) - x) / norm(x)
+    inverse::Float64 = norm((inv(A) * b) - x) / norm(x)
+
+    return Errors(gauss, inverse)
+end
+
+println("hilbert matrices")
+println()
+
+for n in 1:30
+    A::Matrix{Float64} = hilb(n)
+    errors = solve(A, n)
+    println("n: $n; cond: $(cond(A)); rank: $(rank(A)); gauss error: $(errors.gauss); inverse error: $(errors.inverse)")
+end
+
+println()
+println("random matrices")
+println()
+
+for n in [5, 10, 20]
+    for c in [0, 1, 3, 7, 12, 16]
+        A::Matrix{Float64} = matcond(n, 10.0^c)
+        errors = solve(A, n)
+        println("n: $n; c: $c; cond: $(cond(A)); rank: $(rank(A)); gauss error: $(errors.gauss); inverse error: $(errors.inverse)")
+    end
+end

l2/hilb.jl

@@ -0,0 +1,16 @@
+# https://cs.pwr.edu.pl/zielinski/lectures/scna/hilb.jl
+function hilb(n::Int)
+# Function generates the Hilbert matrix  A of size n,
+#  A (i, j) = 1 / (i + j - 1)
+# Inputs:
+#	n: size of matrix A, n>=1
+#
+#
+# Usage: hilb(10)
+#
+# Pawel Zielinski
+        if n < 1
+         error("size n should be >= 1")
+        end
+        return [1 / (i + j - 1) for i in 1:n, j in 1:n]
+end

l2/matcond.jl

@@ -0,0 +1,22 @@
+# https://cs.pwr.edu.pl/zielinski/lectures/scna/matcond.jl
+using LinearAlgebra
+
+function matcond(n::Int, c::Float64)
+# Function generates a random square matrix A of size n with
+# a given condition number c.
+# Inputs:
+#	n: size of matrix A, n>1
+#	c: condition of matrix A, c>= 1.0
+#
+# Usage: matcond(10, 100.0)
+#
+# Pawel Zielinski
+        if n < 2
+         error("size n should be > 1")
+        end
+        if c< 1.0
+         error("condition number  c of a matrix  should be >= 1.0")
+        end
+        (U,S,V)=svd(rand(n,n))
+        return U*diagm(0 =>[LinRange(1.0,c,n);])*V'
+end