{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example of an ill-conditioned system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "cond (generic function with 15 methods)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import LinearAlgebra\n",
    "import Random\n",
    "Random.seed!(0)\n",
    "\n",
    "# Shorthand notation\n",
    "norm = LinearAlgebra.norm\n",
    "cond = LinearAlgebra.cond"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "# # Example of an ill-conditioned system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative error on b: 2.011335237074438503729410008359924320584757100186531629388143594799196078353068e-17\n",
      "Relative error on A: 1.106086006075478567832268914387623212025069211192183133525679964913483566248121e-17\n"
     ]
    }
   ],
   "source": [
    "# Parameters of linear system\n",
    "A = [1 1; 1 (1-1e-12)]\n",
    "b = [0; 1e-12]\n",
    "\n",
    "# Use much more accurate BigFloat format\n",
    "exact_A = [1 1; 1 (1- BigFloat(\"1e-12\"))]\n",
    "exact_b = [0, BigFloat(\"1e-12\")]\n",
    "\n",
    "# Relative error on data\n",
    "relative_error_b = norm(b - exact_b) / norm(b)\n",
    "relative_error_A = norm(A - exact_A) / norm(A)\n",
    "println(\"Relative error on b: $relative_error_b\")\n",
    "println(\"Relative error on A: $relative_error_A\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative error on x: 2.2121720121379184e-5\n"
     ]
    }
   ],
   "source": [
    "# Exact and approximate solutions\n",
    "x_exact = [1; -1]\n",
    "x_approx = A\\b\n",
    "\n",
    "relative_error_x = norm(x_exact - x_approx)/norm(x_approx)\n",
    "println(\"Relative error on x: $relative_error_x\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.000660598318207e12"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Condition number of A\n",
    "cond_A = cond(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The relative error on the solution is much larger than the relative error on the data. The ratio between the two is of the same order of magnitude as the condition number of the matrix $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LU decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: 1.8246083669574607e-12, time: 0.456947637\n",
      "Error: 9.975383169794571e-13, time: 0.005081298\n",
      "Error: 3.4400222673534816e-11, time: 0.100179736\n",
      "Error: 1.2980101641877466e-10, time: 0.737644245\n",
      "Error: 1.9962768190723483e-9, time: 2.961592128\n"
     ]
    }
   ],
   "source": [
    "function lu(A)\n",
    "    n = size(A)[1]\n",
    "    L = [i == j ? 1.0 : 0.0 for i in 1:n, j in 1:n]\n",
    "    U = copy(A)\n",
    "    for i in 1:n-1\n",
    "        for r in i+1:n\n",
    "            U[i, i] == 0 && error(\"Pivotal entry is zero!\")\n",
    "            ratio = U[r, i] / U[i, i]\n",
    "            L[r, i] = ratio\n",
    "            U[r, i:end] -= U[i, i:end] * ratio\n",
    "        end\n",
    "    end\n",
    "    return L, U\n",
    "end\n",
    "\n",
    "ns = 2 .^ collect(5:9)\n",
    "times = zeros(length(ns))\n",
    "for (i, n) in enumerate(ns)\n",
    "    A = randn(n, n)\n",
    "    result = @timed lu(A)\n",
    "    L, U = result.value\n",
    "    times[i] = result.time\n",
    "    error = norm(L*U - A)\n",
    "    println(\"Error: $(LinearAlgebra.norm(L*U - A)), time: $(result.time)\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LU decomposition with pivoting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "-- n = 32 --\n",
      "Condition numbers without pivoting: κ(L) = 178637.01929136374, κ(U) = 161595.53057380614\n",
      "Condition numbers with    pivoting: κ(L) = 28.691391551934966, κ(U) = 193.19052697115316\n",
      "\n",
      "-- n = 64 --\n",
      "Condition numbers without pivoting: κ(L) = 243901.7856336236, κ(U) = 206804.18071457956\n",
      "Condition numbers with    pivoting: κ(L) = 83.9559265157693, κ(U) = 402.1768188798624\n",
      "\n",
      "-- n = 128 --\n",
      "Condition numbers without pivoting: κ(L) = 5.856999082920966e6, κ(U) = 6.169501793769366e6\n",
      "Condition numbers with    pivoting: κ(L) = 207.60793737470857, κ(U) = 599.3815911475028\n",
      "\n",
      "-- n = 256 --\n",
      "Condition numbers without pivoting: κ(L) = 4.2755305949334687e8, κ(U) = 4.729207635483682e8\n",
      "Condition numbers with    pivoting: κ(L) = 447.0937209341298, κ(U) = 1009.1733665132496\n",
      "\n",
      "-- n = 512 --\n",
      "Condition numbers without pivoting: κ(L) = 1.4221750935427153e8, κ(U) = 1.3737670790079305e8\n",
      "Condition numbers with    pivoting: κ(L) = 989.5470052618462, κ(U) = 3785.056520009892\n"
     ]
    }
   ],
   "source": [
    "function lu_pivot(A)\n",
    "    n = size(A)[1]\n",
    "    L, U = zeros(n, 0), copy(A)\n",
    "    P = [i == j ? 1.0 : 0.0 for i in 1:n, j in 1:n]\n",
    "    function swap_rows!(i, j, matrices...)\n",
    "        for M in matrices\n",
    "            M_row_i = M[i, :]\n",
    "            M[i, :] = M[j, :]\n",
    "            M[j, :] = M_row_i\n",
    "        end\n",
    "    end\n",
    "    for i in 1:n-1\n",
    "\n",
    "        # Pivoting\n",
    "        index_row_pivot = i - 1 + argmax(abs.(U[i:end, i]))\n",
    "        swap_rows!(i, index_row_pivot, U, L, P)\n",
    "\n",
    "        # Usual Gaussian transformation\n",
    "        c = [zeros(i-1); 1.0; zeros(n-i)]\n",
    "        for r in i+1:n\n",
    "            ratio = U[r, i] / U[i, i]\n",
    "            c[r] = ratio\n",
    "            U[r, i:end] -= U[i, i:end] * ratio\n",
    "        end\n",
    "        L = [L c]\n",
    "    end\n",
    "    L = [L [zeros(n-1); 1.0]]\n",
    "    return P, L, U\n",
    "end\n",
    "\n",
    "ns = 2 .^ collect(5:9)\n",
    "times = zeros(length(ns))\n",
    "for (i, n) in enumerate(ns)\n",
    "    A = randn(n, n)\n",
    "    L, U = lu(A)\n",
    "    P, Lpivot, Upivot = lu_pivot(A)\n",
    "    println(\"\\n-- n = $n --\")\n",
    "    println(\"Condition numbers without pivoting: κ(L) = $(cond(L)), κ(U) = $(cond(U))\")\n",
    "    println(\"Condition numbers with    pivoting: κ(L) = $(cond(Lpivot)), κ(U) = $(cond(Upivot))\")\n",
    "end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reducing the bandwidth via permutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8×8 SparseArrays.SparseMatrixCSC{Float64,Int64} with 24 stored entries:\n",
       "  [1, 1]  =  10.0\n",
       "  [2, 1]  =  1.0\n",
       "  [8, 1]  =  1.0\n",
       "  [1, 2]  =  1.0\n",
       "  [2, 2]  =  10.0\n",
       "  [3, 2]  =  1.0\n",
       "  [2, 3]  =  1.0\n",
       "  [3, 3]  =  10.0\n",
       "  [4, 3]  =  1.0\n",
       "  [3, 4]  =  1.0\n",
       "  [4, 4]  =  10.0\n",
       "  [5, 4]  =  1.0\n",
       "  [4, 5]  =  1.0\n",
       "  [5, 5]  =  10.0\n",
       "  [6, 5]  =  1.0\n",
       "  [5, 6]  =  1.0\n",
       "  [6, 6]  =  10.0\n",
       "  [7, 6]  =  1.0\n",
       "  [6, 7]  =  1.0\n",
       "  [7, 7]  =  10.0\n",
       "  [8, 7]  =  1.0\n",
       "  [1, 8]  =  1.0\n",
       "  [7, 8]  =  1.0\n",
       "  [8, 8]  =  10.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "8×8 SparseArrays.SparseMatrixCSC{Float64,Int64} with 24 stored entries:\n",
       "  [1, 1]  =  10.0\n",
       "  [2, 1]  =  1.0\n",
       "  [3, 1]  =  1.0\n",
       "  [1, 2]  =  1.0\n",
       "  [2, 2]  =  10.0\n",
       "  [4, 2]  =  1.0\n",
       "  [1, 3]  =  1.0\n",
       "  [3, 3]  =  10.0\n",
       "  [5, 3]  =  1.0\n",
       "  [2, 4]  =  1.0\n",
       "  [4, 4]  =  10.0\n",
       "  [6, 4]  =  1.0\n",
       "  [3, 5]  =  1.0\n",
       "  [5, 5]  =  10.0\n",
       "  [7, 5]  =  1.0\n",
       "  [4, 6]  =  1.0\n",
       "  [6, 6]  =  10.0\n",
       "  [8, 6]  =  1.0\n",
       "  [5, 7]  =  1.0\n",
       "  [7, 7]  =  10.0\n",
       "  [8, 7]  =  1.0\n",
       "  [6, 8]  =  1.0\n",
       "  [7, 8]  =  1.0\n",
       "  [8, 8]  =  10.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import SparseArrays\n",
    "A = LinearAlgebra.diagm(0 => 10*ones(8), 1 => ones(7), -1 => ones(7), 7 => ones(1), -7 => ones(1))\n",
    "A = SparseArrays.SparseMatrixCSC(A)\n",
    "display(A)\n",
    "\n",
    "# Construct permutation matrix obtained by Cuthill-McKee\n",
    "σ = [1, 2, 4, 6, 8, 7, 5, 3]\n",
    "Pσ = SparseArrays.spzeros(8, 8)\n",
    "for j in 1:8\n",
    "    Pσ[σ[j], j] = 1\n",
    "end\n",
    "\n",
    "# Matrix after permutation\n",
    "new_A = Pσ*A*Pσ'\n",
    "display(new_A)"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all"
  },
  "kernelspec": {
   "display_name": "Julia 1.5.3",
   "language": "julia",
   "name": "julia-1.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}