{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "97d1d53f",
   "metadata": {},
   "source": [
    "In lecture 8 of MCS 471, we consider the LU factorization of a matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2b02703",
   "metadata": {},
   "source": [
    "The module ``LinearAlgebra`` is builtin with Julia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5eb47dcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bfac1ca",
   "metadata": {},
   "source": [
    "For formatting matrices, we add the following definition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1cf66238",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Printf\n",
    "Base.show(io::IO, f::Float64) = @printf(io, \"%.3e\", f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fea04d98",
   "metadata": {},
   "source": [
    "# 1. LU Factorization Without Pivoting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02abb915",
   "metadata": {},
   "source": [
    "The function below illustrates the formulas for LU factorization without pivoting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "55bc9ab5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lufac"
      ]
     },
     "execution_count": 3,
     "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": "markdown",
   "id": "0eb2d958",
   "metadata": {},
   "source": [
    "We test the function ``lufac`` on a random 3-by-3 matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e1351018",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 7.105e-01  4.465e-02  6.740e-01\n",
       " 6.425e-01  7.774e-01  8.954e-01\n",
       " 6.570e-01  3.858e-01  3.651e-02"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(3, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "cbb4c26f",
   "metadata": {},
   "outputs": [],
   "source": [
    "L, U = lufac(A);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d50b4219",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 1.000e+00  0.000e+00  0.000e+00\n",
       " 9.043e-01  1.000e+00  0.000e+00\n",
       " 9.248e-01  4.674e-01  1.000e+00"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7be7fd54",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 7.105e-01  4.465e-02   6.740e-01\n",
       " 0.000e+00  7.370e-01   2.859e-01\n",
       " 0.000e+00  0.000e+00  -7.205e-01"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "553307a6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 7.105e-01  4.465e-02  6.740e-01\n",
       " 6.425e-01  7.774e-01  8.954e-01\n",
       " 6.570e-01  3.858e-01  3.651e-02"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L*U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c4ba7480",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 0.000e+00  0.000e+00  0.000e+00\n",
       " 0.000e+00  0.000e+00  0.000e+00\n",
       " 0.000e+00  0.000e+00  0.000e+00"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = A - L*U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "53ab42ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.000e+00"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(D)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34e5d480",
   "metadata": {},
   "source": [
    "# 2. Solving Linear Systems with Blackslash"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48330545",
   "metadata": {},
   "source": [
    "The fastest way to solve a linear system is with the backslash operator."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbafb10b",
   "metadata": {},
   "source": [
    "We make a random linear system that has a vector of ones as solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "fcc38783",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 2.631e-02  1.310e-01  1.799e-01\n",
       " 5.635e-02  2.665e-01  7.534e-01\n",
       " 9.833e-01  7.675e-02  6.986e-01"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(3, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "424a1991",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 1.000e+00\n",
       " 1.000e+00\n",
       " 1.000e+00"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = ones(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ab1dd510",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 3.372e-01\n",
       " 1.076e+00\n",
       " 1.759e+00"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = A*x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2ca7e5b",
   "metadata": {},
   "source": [
    "The application of the backslash operator, as below, solves the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a6bcd903",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 1.000e+00\n",
       " 1.000e+00\n",
       " 1.000e+00"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = A\\b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7b806223",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.220e-16"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(x - y)"
   ]
  }
 ],
 "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
}