08.2.5 Alternative Gauss Jordon Algorithm - Answer.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "8.2.5 Alternative Gauss Jordon Algorithm"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this notebook, you will implement the alternative Gauss Jordan algorithm that overwrites $ A $ in one sweep with the identity matrix and $ B $ with the inverse of the original matrix $ A $.\n",
      "\n",
      "<font color=red> Be sure to make a copy!!!! </font>\n",
      "\n",
      "<h2> First, let's create a matrix $ A $ and set $ B $ to the identity. </h2>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import laff\n",
      "import flame\n",
      "\n",
      "L = np.matrix( ' 1, 0, 0. 0;\\\n",
      "                -2, 1, 0, 0;\\\n",
      "                 1,-3, 1, 0;\\\n",
      "                 2, 3,-1, 1' )\n",
      "\n",
      "U = np.matrix( ' 2,-1, 3,-2;\\\n",
      "                 0,-2, 1,-1;\\\n",
      "                 0, 0, 1, 2;\\\n",
      "                 0, 0, 0, 3' )\n",
      "\n",
      "A = L * U\n",
      "Aold = np.matrix( np.copy( A ) ) \n",
      "B = np.matrix( np.eye( 4 ) )\n",
      "\n",
      "print( 'A = ' )\n",
      "print( A )\n",
      "\n",
      "print( 'B = ' )\n",
      "print( B )\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> Implement the alternative Gauss-Jordan algorithm from 8.2.5 </h2>\n",
      "\n",
      "Here is the algorithm:\n",
      "\n",
      "<img src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/8_2_5_Algorithm.png\" alt=\"Alternative Gauss-Jordan algorithm\" width=100%>\n",
      "    \n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working, and then the one where you are working. </font>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Create the routine\n",
      "<code> GJ_Inverse_alt </code>\n",
      "with the <a href=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/index.html\"> Spark webpage</a> for the algorithm"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import flame\n",
      "\n",
      "def GJ_Inverse_alt(A, B):\n",
      "\n",
      "    ATL, ATR, \\\n",
      "    ABL, ABR  = flame.part_2x2(A, \\\n",
      "                               0, 0, 'TL')\n",
      "\n",
      "    BTL, BTR, \\\n",
      "    BBL, BBR  = flame.part_2x2(B, \\\n",
      "                               0, 0, 'TL')\n",
      "\n",
      "    while ATL.shape[0] < A.shape[0]:\n",
      "\n",
      "        A00,  a01,     A02,  \\\n",
      "        a10t, alpha11, a12t, \\\n",
      "        A20,  a21,     A22   = flame.repart_2x2_to_3x3(ATL, ATR, \\\n",
      "                                                       ABL, ABR, \\\n",
      "                                                       1, 1, 'BR')\n",
      "\n",
      "        B00,  b01,    B02,  \\\n",
      "        b10t, beta11, b12t, \\\n",
      "        B20,  b21,    B22   = flame.repart_2x2_to_3x3(BTL, BTR, \\\n",
      "                                                      BBL, BBR, \\\n",
      "                                                      1, 1, 'BR')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        # a01 := a01 / alpha11\n",
      "        # a21 := a21 / alpha11\n",
      "        laff.invscal( alpha11, a01 )\n",
      "        laff.invscal( alpha11, a21 )\n",
      "        \n",
      "        # A02 := A02 - a01 * a12t\n",
      "        laff.ger( -1.0, a01, a12t, A02 )\n",
      "        # A22 := A22 - a21 * a12t\n",
      "        laff.ger( -1.0, a21, a12t, A22 )\n",
      "        \n",
      "        # B00 := B00 - a01 * b01t\n",
      "        laff.ger( -1.0, a01, b10t, B00 )\n",
      "        # B20 := B20 - a21 * b01t\n",
      "        laff.ger( -1.0, a21, b10t, B20 )\n",
      "        \n",
      "        # b01 := - a01 (= - u01 in the discussion)\n",
      "        laff.copy( a01, b01 )\n",
      "        laff.scal( -1.0, b01 )\n",
      "        # b21 := - a21 (= - l21 in the discussion)\n",
      "        laff.copy( a21, b21 )\n",
      "        laff.scal( -1.0, b21 )\n",
      "        \n",
      "        # a12t:= a21t / alpha11  \n",
      "        laff.invscal( alpha11, a12t )\n",
      "        # b10t:= b10t / alpha11 \n",
      "        laff.invscal( alpha11, b10t )\n",
      "        \n",
      "        # beta11 := 1.0 / alpha11\n",
      "        laff.invscal( alpha11, beta11 )\n",
      "       \n",
      "        # a01 = 0  (zero vector)\n",
      "        # alpha11 = 1\n",
      "        # a21 = 0  (zero vector)\n",
      "        laff.zerov( a01 )\n",
      "        laff.onev( alpha11 )\n",
      "        laff.zerov( a21 )\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        ATL, ATR, \\\n",
      "        ABL, ABR  = flame.cont_with_3x3_to_2x2(A00,  a01,     A02,  \\\n",
      "                                               a10t, alpha11, a12t, \\\n",
      "                                               A20,  a21,     A22,  \\\n",
      "                                               'TL')\n",
      "\n",
      "        BTL, BTR, \\\n",
      "        BBL, BBR  = flame.cont_with_3x3_to_2x2(B00,  b01,    B02,  \\\n",
      "                                               b10t, beta11, b12t, \\\n",
      "                                               B20,  b21,    B22,  \\\n",
      "                                               'TL')\n",
      "\n",
      "    flame.merge_2x2(ATL, ATR, \\\n",
      "                    ABL, ABR, A)\n",
      "\n",
      "    flame.merge_2x2(BTL, BTR, \\\n",
      "                    BBL, BBR, B)\n",
      "\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h3> Test the routine </h3>\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working, and then the one where you are working. </font>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "GJ_Inverse_alt( A, B )\n",
      "\n",
      "print( A )\n",
      "print( B )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Matrix $ A $ should now be an identity matrix and $ B $ should no longer be an identity matrix.\n",
      "\n",
      "Check if $ B $ now equals (approximately) the inverse of the original matrix $ A $:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print( Aold * B )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}