{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Partial Pivoting\n",
    "\n",
    "**References:**\n",
    "\n",
    "- Section 2.4.1 *Partial Pivoting* of [Sauer](../frontmatter/references.html#Sauer)\n",
    "\n",
    "- Section 6.2 *Pivoting Stratgies* of [Burden&Faires](../frontmatter/references.html#Burden-Faires)\n",
    "\n",
    "- Section 7.1 of [Chenney&Kincaid](../frontmatter/references.html#Chenney-Kincaid)\n",
    "\n",
    "**Note:**\n",
    "Some references describe the method of *scaled* partial pivoting, but here we present instead a version without the \"scaling\", because not only is it simpler, but modern research shows that it is esentially always as good, once the problem is set up in a \"sane\" way."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "The basic row reduction method can fail due to divisoion by zero (and to have very large rouding errors when a denominator is extremely close to zero.\n",
    "A more robust modification is to swap the order of the equations to avaid these problems: *partial pivotng*.\n",
    "Here we look at a particularly robust version of this strategy, *Maximal Element Partial Pivoting*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What can go wrong with Naive Gaussian Elimination?\n",
    "\n",
    "We have noted two problems with the naive algorithm for Gaussian elimination: total failure due the division be zero, and loss of precision due to dividing by very small values — or more preciselt calculations the lead to intermediate values far larger than the final results.\n",
    "The culprits in all cases are the same: the denominators are first the *pivot elements* $a_{k,k}^{(k-1)}$ in evaluation of $l_{i,k}$ during row reduction and then the $u_{k,k}$ in back substitution.\n",
    "Further, those $a_{k,k}^{(k-1)}$ are the final updated values at indices $(k,k)$, so are the same as $u_{k,k}$.\n",
    "Thus it is exactly these *main diagonal elements* that we must deal with."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The basic fix: Partial Pivoting\n",
    "\n",
    "The basic strategy is that at step $k$, we can swap equation $k$ with any equation $i$, $i > k$.\n",
    "Note that this involves swapping those rows of array `A` and also those elements of the array `b` for the right-hand side: $b_k \\leftrightarrow b_i$.\n",
    "\n",
    "This approach of swapping equations (swapping rows in arrays `A` and`b`) is called **pivoting**, or more specifically **partial pivoting**, to distinguish from the more elaborate strategy where to columns of `A` are also reordered (which is equivalent to reordeting the unknowns in the equations).\n",
    "The row that is swapped with row $k$ is sometimes called the **pivot row**, and the new denominator is the corresponding **pivot element**.\n",
    "\n",
    "This approach is robust so long as one is using exact arithmetic: it works for any well-posed system because so long as the $Ax = b$ has a unique solution — so that the original matrix $A$ is non-singular — at least one of the $a_{i,k}^{(k-1)}, i \\geq k$ will be non-zero, and thus the swap will give a new element in position $(k,k)$ that is non-zero.\n",
    "(I will stop caring about superscripts to distinguish updates, but if you wish to, the elements of the new row $k$ could be called either $a_{k,j}^{(k)}$ or even $u_{k,j}$, since those values are in their final state.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Handling rounding error: Maximal Element Partial Pivoting\n",
    "\n",
    "The final refinement is to seek the smallest possible magnitudes for intermediate values, and thus the smallest absolute errors in them, by making the multipliers $l_{i,k}$ small, in turn by making the denominator $a_{k,k}^{(k-1)} = u_{k,k}$ as large as possible in magnitude:\n",
    "\n",
    "*At step $k$, choose the pivot row $p_k \\geq k$ so that $|a_{p_k,k}^{(k-1)}| \\geq |a_{i,k}^{(k-1)}|$ for all $i \\geq k$. If there is more that one such element of largest magnitude, use the lowest value:\n",
    "in particular, if $p_k = k$ works, use it and do not swap!*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Julia note 1.** Swapping values, including rows.\n",
    "\n",
    "In Julia (as in Python) the value of two variables `a` and `b` can be swapped via tuples with `(a, b) = (b, a)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 1)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = 1\n",
    "b = 2\n",
    "(a, b) = (b, a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combined with array slicing, this can also be used to swap rows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " 11  12  13\n",
       " 21  22  23\n",
       " 31  32  33"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [11 12 13 ; 21 22 23 ; 31 32 33]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " 21  22  23\n",
       " 11  12  13\n",
       " 31  32  33"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A[1,:], A[2,:]) = (A[2,:], A[1,:])\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "### Exercise 1\n",
    "\n",
    "Explain why we cannot just swap the relevant elements of rows $k$ and $p$ with:\n",
    "\n",
    "    for j in 1:n\n",
    "        A[k,j] = A[p,j]\n",
    "        A[p,j] = A[k,j]\n",
    "    end\n",
    "\n",
    "or in vectorized form:\n",
    "\n",
    "    A[k,:] = A[p,:]\n",
    "    A[p,:] = A[k,:]\n",
    "\n",
    "Describe what happens instead."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some demonstrations\n",
    "\n",
    "No row reduction is done here, so entire rows are swapped rather than just the elements from column $k$ onward:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, get the matrix pretty-printer seen in {doc}`linear-equations-1-row-reduction-julia`;\n",
    "this time from {doc}`NumericalMethods`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"NumericalMethods.jl\")\n",
    "using .NumericalMethods: printmatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initially A is\n",
      "[ 1 -6 2 \n",
      "  3 5 -6 \n",
      "  4 2 7 ]\n"
     ]
    }
   ],
   "source": [
    "A = [1 -6 2 ; 3 5 -6 ; 4 2 7]\n",
    "n = 3\n",
    "println(\"Initially A is\")\n",
    "printmatrix(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After swapping rows 1 <-> 3 using slicing and a temporary row, A is\n",
      "[ 4 2 7 \n",
      "  3 5 -6 \n",
      "  1 -6 2 ]\n"
     ]
    }
   ],
   "source": [
    "k = 1\n",
    "p = 3\n",
    "temp = copy(A[k,:])\n",
    "A[k,:] = A[p,:]\n",
    "A[p,:] = temp\n",
    "println(\"After swapping rows 1 <-> 3 using slicing and a temporary row, A is\")\n",
    "printmatrix(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After swapping rows 2 <-> 3 using a loop and tuples of elements (no temp) A is:\n",
      "[ 4 2 7 \n",
      "  1 -6 2 \n",
      "  3 5 -6 ]\n"
     ]
    }
   ],
   "source": [
    "k = 2\n",
    "p = 3\n",
    "for j in 1:n\n",
    "    ( A[k,j] , A[p,j] ) = ( A[p,j] , A[k,j] )\n",
    "end\n",
    "println(\"After swapping rows 2 <-> 3 using a loop and tuples of elements (no temp) A is:\")\n",
    "printmatrix(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After swapping rows 1 <-> 2 using tuples of slices (no loop or temp) A is:\n",
      "[ 1 -6 2 \n",
      "  4 2 7 \n",
      "  3 5 -6 ]\n"
     ]
    }
   ],
   "source": [
    "k = 1\n",
    "p = 2\n",
    "( A[k,:] , A[p,:] ) = ( copy(A[p,:]) , copy(A[k,:]) )\n",
    "println(\"After swapping rows 1 <-> 2 using tuples of slices (no loop or temp) A is:\")\n",
    "printmatrix(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "This work is licensed under [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0/)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.1",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}