{
 "metadata": {
  "name": "",
  "signature": "sha256:f9f0225755125b52906cbeac69554c162ab9fb69a6f78503294c83f65d1daf92"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Gaussian Elimination"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 77
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This function performs LU decomposition of given matrix A."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def LU(A):\n",
      "    n = A.shape[0]\n",
      "    L = np.zeros((n,n))\n",
      "    U = np.array(A)\n",
      "    for i in range(n-1):\n",
      "        L[i,i] = 1.0\n",
      "        L[i+1:n,i] = U[i+1:n,i]/U[i,i]\n",
      "        for j in range(i+1,n):\n",
      "            U[j,i+1:n] = U[j,i+1:n] - L[j,i] * U[i,i+1:n]\n",
      "        U[i+1:n,i] = 0.0\n",
      "    L[n-1,n-1] = 1.0\n",
      "    return L,U"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 86
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This performs performs solution of\n",
      "$$\n",
      "LUx=b\n",
      "$$\n",
      "in two steps. In first step, solve\n",
      "$$\n",
      "Ly = b\n",
      "$$\n",
      "In second step, solve\n",
      "$$\n",
      "Ux = y\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def LUSolve(L,U,b):\n",
      "    n = L.shape[0]\n",
      "    # solve Ly = b\n",
      "    y = 0*b\n",
      "    for i in range(n):\n",
      "        y[i] = (b[i] - L[i,0:i].dot(y[0:i]))/L[i,i]\n",
      "    # solve Ux = y\n",
      "    x = 0*b\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": 86
    },
    {
     "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 = LU(A)\n",
      "print A\n",
      "print L\n",
      "print U\n",
      "print A - L.dot(U)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 0.23000755  0.32161034  0.37253327]\n",
        " [ 0.87099231  0.62303385  0.06389315]\n",
        " [ 0.24376526  0.92755117  0.03188404]]\n",
        "[[ 1.          0.          0.        ]\n",
        " [ 3.78679872  1.          0.        ]\n",
        " [ 1.05981414 -0.98632268  1.        ]]\n",
        "[[ 0.23000755  0.32161034  0.37253327]\n",
        " [ 0.         -0.59483979 -1.34681537]\n",
        " [ 0.          0.         -1.69132653]]\n",
        "[[ 0.  0.  0.]\n",
        " [ 0.  0.  0.]\n",
        " [ 0.  0.  0.]]\n"
       ]
      }
     ],
     "prompt_number": 87
    },
    {
     "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,b)\n",
      "print A.dot(x) - b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[  0.00000000e+00   0.00000000e+00   2.22044605e-16]\n"
       ]
      }
     ],
     "prompt_number": 88
    }
   ],
   "metadata": {}
  }
 ]
}