{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Gaussian elimination II"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Importing functions we wrote in the last class : packages."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In python, one can just store a bunch of function definitions in a file. Then the file can be imported like any other package. We have stored our functions in `gauss_elim_old_fn.py`. We import it as"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from gauss_elim_old_fn import *"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[1.0, 2.0, 3.0], [2.0, 1.0, 2.0], [3.0, 2.0, 1.0]]\n",
        "         1          2          3\n",
        "         2          1          2\n",
        "         3          2          1\n",
        " 1  2  3\n",
        " 2  1  2\n",
        " 3  2  1\n",
        "--------------------------\n",
        " 4  8 12\n",
        " 2  1  2\n",
        " 3  2  1\n",
        "-------------------------\n",
        " 4  8 12\n",
        " 2  1  2\n",
        " 3  2  1\n",
        " 4  8 12\n",
        " 2  1  2\n",
        " 3  2  1\n",
        "--------------------\n",
        " 0  6  8\n",
        " 2  1  2\n",
        " 3  2  1\n",
        "---------------------\n",
        " 0  6  8\n",
        " 2  1  2\n",
        " 3  2  1\n",
        "(1, 1)\n",
        "The matrix is the zero matrix.\n",
        "(-1, -1)\n",
        "The matrix is the zero matrix.\n",
        "(-1, -1)\n",
        "Please input a natural number.\n",
        "None\n",
        "1\n",
        "1\n",
        "2\n",
        "3\n",
        "5\n",
        "8\n",
        "13\n",
        "21\n",
        "34\n",
        "55\n",
        "Please input an positive integer.\n",
        "None\n",
        "Please input a natural number.\n",
        "None\n",
        "1\n",
        "1\n",
        "2\n",
        "3\n",
        "5\n",
        "8\n",
        "13\n",
        "21\n",
        "34\n",
        "55\n",
        "Please input a natural number.\n",
        "None\n",
        "9.33262154439e+157 Expected a non-negative integer.\n",
        "None Expected a non-negative integer.\n",
        "None 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000\n",
        "(2, 1)\n",
        "  1   2   3   4\n",
        "  5   6   7   8\n",
        "  9  10  11  12\n",
        "\n",
        "  1   2   3   4\n",
        "  5   6   7   8\n",
        "  9  10  11  12\n",
        "  1   2   3   4\n",
        "  5   6   7   8\n",
        "  9  10  11  12\n",
        "---------------\n",
        "  1   2   3   4\n",
        "  5   6   7  -2\n",
        "  9  10  11  -3\n",
        "   0.00000    0.50000    1.00000    0.00000\n",
        "   2.00000   11.00000    0.00000    0.00000\n",
        "   1.00000    0.00000    0.00000    9.00000\n",
        "--------------------------------------------------\n",
        "The matrix is empty.\n",
        "Nothing to copy: S is empty\n",
        "   1.00000    5.50000    0.00000    0.00000\n",
        "   0.00000    1.00000    2.00000    0.00000\n",
        "   0.00000    0.00000    1.00000    0.81818\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def gauss_elim_prelim(matrix) :\n",
      "    row = 0\n",
      "    row_offset = 0\n",
      "    col_offset = 0\n",
      "    \n",
      "    (r, c) = first_non_zero(matrix)\n",
      "    if r == -1 :\n",
      "        print \"Nothing to do: The matrix is zero.\"\n",
      "    else :\n",
      "        matrix = elem1(matrix, 0, r)\n",
      "        matrix = sweep(matrix, (0, c))\n",
      "        print_matrix(matrix, \"%6.2f\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gauss_elim_prelim([[0,0,2,2],\\\n",
      "            [1,4,4,4],\\\n",
      "            [0,1,3,3]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "  1.00   4.00   4.00   4.00\n",
        "  0.00   0.00   2.00   2.00\n",
        "  0.00   1.00   3.00   3.00\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def gauss_elim(matrix) :\n",
      "    row = 0\n",
      "    row_offset = 0\n",
      "    col_offset = 0\n",
      "    no_of_rows = len(matrix)\n",
      "    if no_of_rows == 0 :\n",
      "        print \"The matrix is empty.\"\n",
      "    else :\n",
      "        no_of_cols = len(matrix[0])\n",
      "        if no_of_cols == 0 :\n",
      "            print \"The matrix is empty.\"\n",
      "        else :\n",
      "            while row < no_of_rows :\n",
      "                submatrix = []\n",
      "                # Copy the submatrix to the right, bottom of (row_offset, col_offset)\n",
      "                # to submatrix\n",
      "                for i in range(no_of_rows - row_offset) :\n",
      "                    row_submatrix = []\n",
      "                    for j in range(no_of_cols - col_offset) :\n",
      "                        row_submatrix.append(matrix[i + row_offset][j + col_offset])\n",
      "                    submatrix.append(row_submatrix)\n",
      "                    #print_matrix(submatrix, \"%6.2f\")\n",
      "                # print row, row_offset, col_offset, submatrix, matrix\n",
      "                # print \"-\"*20\n",
      "                                     \n",
      "                (r, c) = first_non_zero(submatrix)\n",
      "                if r != -1 :\n",
      "                    # print \"Nothing to do: The (sub)matrix is zero.\"\n",
      "                    # else :\n",
      "                    r += row_offset\n",
      "                    c += col_offset\n",
      "                    matrix = elem1(matrix, row_offset, r)\n",
      "                    matrix = sweep(matrix, (row_offset, c))\n",
      "                    row_offset += 1\n",
      "                    col_offset = c+1\n",
      "                row += 1\n",
      "    return matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "B = read_matrix(\"files/B.txt\")\n",
      "print_matrix(B, \"%6.2f\")\n",
      "print \"-\"*40\n",
      "print_matrix(gauss_elim(B), \"%6.2f\")\n",
      "print \"=\"*40\n",
      "C = [[0,0,0,0], [0,1,1,0], [0,0,1,0]]\n",
      "print_matrix(C, \"%3d\")\n",
      "print \"-\"*40\n",
      "print_matrix(gauss_elim(C), \"%3d\")\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "  0.00   2.00   3.00   2.00\n",
        "  1.00   4.00   5.00   5.00\n",
        "  3.00   0.00   2.00   6.00\n",
        "----------------------------------------\n",
        "  1.00   0.00   0.00   1.60\n",
        "  0.00   1.00   0.00   0.10\n",
        "  0.00   0.00   1.00   0.60\n",
        "========================================\n",
        "  0   0   0   0\n",
        "  0   1   1   0\n",
        "  0   0   1   0\n",
        "----------------------------------------\n",
        "The matrix is the zero matrix.\n",
        "  0   1   0   0\n",
        "  0   0   1   0\n",
        "  0   0   0   0\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Solving a system of linear equations"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First suppose we are given an equation of the type $Ax = b$ where $A$ is an $m \\times n$ matrix and $x$ is an $n$-vector and $b$ is an $m$-vector. We first form the matrix $A:b$ which is obtained by sticking $b$ after $A$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def join(A, B) :\n",
      "    \"\"\"Returns A:B\"\"\"\n",
      "    no_rows_A = len(A)\n",
      "    no_rows_B = len(B)\n",
      "    \n",
      "    if no_rows_A != no_rows_B :\n",
      "        print \"Matrices should have the same number of rows to be joined.\"\n",
      "        retM = None\n",
      "    else :\n",
      "        retM = []\n",
      "        if no_rows_A == 0 :\n",
      "            retM = A\n",
      "        else :\n",
      "            no_cols_A = len(A[0])\n",
      "            no_cols_B = len(B[0])\n",
      "            for i in range(no_rows_A) :\n",
      "                Mrow = []\n",
      "                for j in range(no_cols_A) :\n",
      "                    Mrow.append(A[i][j])\n",
      "                for k in range(no_cols_B) :\n",
      "                    Mrow.append(B[i][k])\n",
      "                retM.append(Mrow)\n",
      "        return retM"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "join([[4],[5]], [[4,5],[9,8]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "[[4, 4, 5], [5, 9, 8]]"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### To solve $Ax = b$\n",
      "Gaussian elimination of $A:b$, if $A$ is *invertible* is just finding a matrix $P$ such that $P(A:b)$ is of the form $I_3:c$. Then $P$ will be the inverse of $A$ and $c = Pb$ will be the solution."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Let A = \n",
      "A = [[0, 2, 3], [1, 0, 1], [0, 0, 1]]\n",
      "b = [[5], [2], [1]]\n",
      "# A:b = \n",
      "Ab = join(A, b)\n",
      "print_matrix(gauss_elim(Ab), \"%5.2g\")\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "    1     0     0     1\n",
        "    0     1     0     1\n",
        "    0     0     1     1\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Thus $x = \\begin{pmatrix}1\\\\1\\\\1\\end{pmatrix}$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To invert a matrix, one just appends identity on the right and does Gaussian elimination."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def idmat(n) :\n",
      "    \"\"\"returns the identity matrix of size n\"\"\"\n",
      "    rmat = []\n",
      "    for i in range(n) :\n",
      "        row = []\n",
      "        for j in range(n) :\n",
      "            if i == j :\n",
      "                row.append(1)\n",
      "            else :\n",
      "                row.append(0)\n",
      "        rmat.append(row)\n",
      "    return rmat"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print_matrix(idmat(4), \"%2d\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1  0  0  0\n",
        " 0  1  0  0\n",
        " 0  0  1  0\n",
        " 0  0  0  1\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def invert(A) :\n",
      "    size = len(A)\n",
      "    bigmat = join(A,idmat(size))\n",
      "    biginv = gauss_elim(bigmat)\n",
      "    inv = []\n",
      "    for i in range(size) :\n",
      "        row = []\n",
      "        for j in range(size) :\n",
      "            row.append(biginv[i][j + size])\n",
      "        inv.append(row)\n",
      "    return inv"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "UT1 = [[1,1,1],[0,1,1],[0,0,1]]\n",
      "print \"Inverse of\"\n",
      "print_matrix(UT1, \"%3d\")\n",
      "print \"is\"\n",
      "print_matrix(invert(UT1), \"%3d\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Inverse of\n",
        "  1   1   1\n",
        "  0   1   1\n",
        "  0   0   1\n",
        "is\n",
        "  1  -1   0\n",
        "  0   1  -1\n",
        "  0   0   1\n"
       ]
      }
     ],
     "prompt_number": 12
    }
   ],
   "metadata": {}
  }
 ]
}