{
 "metadata": {
  "name": "",
  "signature": "sha256:812d54e00eb07d70eb4212f216f75c30c0992c010fee7905c0ec79fdd91bff37"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Gaussian Elimination with row pivoting"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This function performs LU decomposition of given matrix A with pivoting to obtain\n",
      "$$\n",
      "PA = LU\n",
      "$$\n",
      "where $P$ is a permutation matrix. Here we dont store $P$ as a matrix, but only as a vector."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def LU(A):\n",
      "    n = A.shape[0]\n",
      "    L = np.identity(n)\n",
      "    P = np.arange(n,dtype=int) # Permutation matrix\n",
      "    U = np.array(A)\n",
      "    for k in range(n-1):\n",
      "        i = np.argmax(abs(U[k:n,k]))\n",
      "        i = i+k\n",
      "        U[[k,i],k:n] = U[[i,k],k:n] # swap row i and k\n",
      "        L[[k,i],0:k] = L[[i,k],0:k] # swap row i and k\n",
      "        P[[k,i]] = P[[i,k]]         # swap row i and k\n",
      "        for j in range(k+1,n):\n",
      "            L[j,k] = U[j,k]/U[k,k]\n",
      "            U[j,k:n] = U[j,k:n] - L[j,k]*U[k,k:n]\n",
      "    return L,U,P"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This performs solution of\n",
      "$$\n",
      "LUx = Pb\n",
      "$$\n",
      "in two steps. In first step, solve\n",
      "$$\n",
      "Ly = Pb\n",
      "$$\n",
      "In second step, solve\n",
      "$$\n",
      "Ux = y\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def LUSolve(L,U,P,b):\n",
      "    n = L.shape[0]\n",
      "    # solve Ly = Pb\n",
      "    pb = b[P]\n",
      "    y = 0.0*pb\n",
      "    for i in range(n):\n",
      "        y[i] = (pb[i] - L[i,0:i].dot(y[0:i]))/L[i,i]\n",
      "    #solve Ux = y\n",
      "    x = 0.0*pb\n",
      "    for i in range(n-1,-1,-1):\n",
      "        x[i] = (y[i] - U[i,i+1:n].dot(x[i+1:n]))/U[i,i]\n",
      "    return x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we test the above function for LU decomposition."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 3\n",
      "A = np.random.rand(n,n)\n",
      "L,U,P = LU(A)\n",
      "print A\n",
      "print L\n",
      "print U\n",
      "print P\n",
      "print P.dot(A) - L.dot(U)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 0.41300456  0.43021602  0.38608235]\n",
        " [ 0.41796129  0.86040814  0.0505991 ]\n",
        " [ 0.09585373  0.27291233  0.11822924]]\n",
        "[[ 1.          0.          0.        ]\n",
        " [ 0.98814069  1.          0.        ]\n",
        " [ 0.22933638 -0.17997989  1.        ]]\n",
        "[[ 0.41796129  0.86040814  0.0505991 ]\n",
        " [ 0.         -0.41998827  0.33608332]\n",
        " [ 0.          0.          0.16711326]]\n",
        "[1 0 2]\n",
        "[[ 0.18675073  0.11563255  0.57194172]\n",
        " [ 0.19170746  0.54582466  0.23645847]\n",
        " [ 0.50885829  0.70312835  0.50431159]]\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Solve the linear system."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "b = np.random.rand(n)\n",
      "x = LUSolve(L,U,P,b)\n",
      "print A.dot(x) - b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[  0.00000000e+00  -1.11022302e-16  -1.11022302e-16]\n"
       ]
      }
     ],
     "prompt_number": 6
    }
   ],
   "metadata": {}
  }
 ]
}