{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7bf7ea8e",
   "metadata": {},
   "source": [
    "Quiz 5 MCS 471, on Friday 23 September 2022, due 10:50am."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac0ac1c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "Consider the linear system $A {\\bf x} = {\\bf b}$ with\n",
    "\n",
    "$$\n",
    "   A = \n",
    "   \\left[\n",
    "     \\begin{array}{cc}\n",
    "        10^{-10} & 1 \\\\\n",
    "        2 & 10^{-10}\n",
    "     \\end{array}\n",
    "   \\right]\n",
    "   \\quad\n",
    "   {\\bf b} = \\left[ \\begin{array}{c} 1 \\\\ 1 \\end{array} \\right].\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0149217e",
   "metadata": {},
   "source": [
    "1. Compute the LU factorization with the ``lu`` of the ``LinearAlgebra`` module and use the output of ``lu`` to compute $\\bf x$.\n",
    "\n",
    "2. Compute the LU factorization *without* pivoting and use this to solve the system.\n",
    "\n",
    "3. Explain the difference between the two solutions, once obtained with and once without pivoting."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00bc10a9",
   "metadata": {},
   "source": [
    "# Answer to question 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "714e6c87",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0b8d657c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " 1.0e-10  1.0\n",
       " 2.0      1.0e-10"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = reshape([1.0e-10, 2.0, 1.0, 1.0e-10], (2,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a2453069",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LU{Float64, Matrix{Float64}, Vector{Int64}}\n",
       "L factor:\n",
       "2×2 Matrix{Float64}:\n",
       " 1.0      0.0\n",
       " 5.0e-11  1.0\n",
       "U factor:\n",
       "2×2 Matrix{Float64}:\n",
       " 2.0  1.0e-10\n",
       " 0.0  1.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L, U, P = lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b88b2edb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 1.0\n",
       " 1.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = [1.0, 1.0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e674646f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 1.0\n",
       " 0.99999999995"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = L\\b[P]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f143415c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 0.49999999995\n",
       " 0.99999999995"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = U\\y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7b78efc8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(b - A*x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e665d822",
   "metadata": {},
   "source": [
    "# Answer to Question 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7be7ddbf",
   "metadata": {},
   "source": [
    "To solve this without pivoting, we apply the function ``lufac`` of Lecture 8."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e7ca9a1c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lufac"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Returns the matrices L and U in an LU factorization of A.\n",
    "\"\"\"\n",
    "function lufac(mat::Array{Float64}) \n",
    "    nbrows, nbcols = size(mat)\n",
    "    low = Matrix{Float64}(I,nbrows, nbcols)\n",
    "    upp = zeros(nbrows, nbcols)\n",
    "    for j=1:nbcols\n",
    "        for i=1:j\n",
    "            upp[i,j] = mat[i,j]\n",
    "            for k=1:i-1\n",
    "                upp[i,j] = upp[i,j] - low[i,k]*upp[k,j]\n",
    "            end\n",
    "        end\n",
    "        for i=j+1:nbrows\n",
    "            low[i, j] = mat[i,j]\n",
    "            for k=1:j-1\n",
    "                low[i, j] = low[i, j] - low[i, k]*upp[k, j]\n",
    "            end\n",
    "            low[i, j] = low[i, j]/upp[j, j]\n",
    "        end\n",
    "    end\n",
    "    return (low, upp)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "95da135b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([1.0 0.0; 2.0e10 1.0], [1.0e-10 1.0; 0.0 -2.0e10])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L2, U2 = lufac(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "581398c8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       "  1.0\n",
       " -1.9999999999e10"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y2 = L2\\b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "dcf891ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 0.5000000413701855\n",
       " 0.99999999995"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x2 = U2\\y2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "34035b2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.284037100736441e-8"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(b - A*x2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eea0dd08",
   "metadata": {},
   "source": [
    "# Anwer to question 3."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2f473b7",
   "metadata": {},
   "source": [
    "The first solution computed with pivoting is ``x``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0856280a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 0.49999999995\n",
       " 0.99999999995"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5244e3a5",
   "metadata": {},
   "source": [
    "Without pivoting the solution is in ``x2``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "801ee8b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 0.5000000413701855\n",
       " 0.99999999995"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "54409292",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.1420185503682205e-8"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(x - x2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06132f5f",
   "metadata": {},
   "source": [
    "We see that the difference is ``4.0e-8`` which is also the computed residual, or the backward error of ``x2``.  The error of ``x2`` is due to the large numbers in the output of the LU factorization computed without row pivoting."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}