08.4.0 - (Optional Enrichment) Blocked LU Factorization - Answer.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Solving $ A x = b $ via Blocked LU Factorization"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this notebook, you will implement a blocked LU factorization, solve a system with a unit lower triangular matrix, and solve a system with an upper triangular matrix.  This notebook culminates in a routine that combines these three steps into a routine that solves $ A x = b $ in a computationally efficient way.\n",
      "\n",
      "<font color=red> Be sure to make a copy!!!! </font>\n",
      "\n",
      "<h2>Preliminaries</h2>\n",
      "\n",
      "Here is a list of laff routines that you might want to use in this notebook:\n",
      "<ul>\n",
      "<li> <code>laff.trsv( uplo, trans, diag, A, b )</code> Solves $Ax = b$ where $x$ and $b$ are vectors.\n",
      "<li> <code>laff.trsm( uplo, trans, diag, A, B )</code> Solves $AX = B$ where $X$ and $B$ are matrices.\n",
      "<li> <code>laff.gemm( alpha, A, B, beta, C )</code> $C := \\alpha A B + \\beta C$\n",
      "</ul>\n",
      "\n",
      "These three methods were added to the laff library after the course started. If you're having problems with git and are manually downloading the notebooks, check the wiki page for notebooks for help on getting updated laff routines as well.\n",
      "\n",
      "\n",
      "And last but not least, __*copy and paste your method from 6.3 Solving A x b via LU factorization and triangular solves into the box below.*__ We'll be using it during this notebook. Recall that it overwrites $A$ with $L$ in the strictly lower triangular part and $U$ in the upper triangular part.\n",
      "<font color=red> Make sure you call the routine LU_unb_var5 </font>\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import flame\n",
      "import laff\n",
      "\n",
      "def LU_unb_var5(A):\n",
      "\n",
      "    ATL, ATR, \\\n",
      "    ABL, ABR  = flame.part_2x2(A, \\\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",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        laff.invscal( alpha11, a21 )        #  a21 := a21 / alpha11\n",
      "        laff.ger( -1.0, a21, a12t, A22 )    #  A22 := A22 - a21 * a12t\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",
      "    flame.merge_2x2(ATL, ATR, \\\n",
      "                    ABL, ABR, A)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> Now, let's create a matrix $ A $ and right-hand side $ b $</h2>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import laff\n",
      "import flame\n",
      "\n",
      "# Form the matrices\n",
      "m = 200\n",
      "\n",
      "L = np.matrix( np.tril(np.random.rand(m,m), -1 ) + np.eye(m) )\n",
      "U = np.matrix( np.triu(np.random.rand(m,m) ) + np.eye(m) * m )\n",
      "\n",
      "A = L * U\n",
      "\n",
      "# Create a large, random solution vector x\n",
      "x = np.matrix( np.random.rand(m,1) )\n",
      "\n",
      "# Store the original value of x\n",
      "xold = np.matrix( np.copy( x ) )\n",
      "\n",
      "# Create a solution vector b so that A x = b\n",
      "b = A * x\n",
      "\n",
      "\n",
      "# Later, we are also going to solve A x = b2.  Here we create that b2:\n",
      "x2 = np.matrix( np.random.rand(m,1) )\n",
      "b2 = A * x2\n",
      "\n",
      "# Set the block size, we'll use it later. \n",
      "# Since we're just messing around with blocked algorithms,\n",
      "# we set this totally arbitrarily\n",
      "nb = 20"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> Implement the blocked LU factorization routine from 6.4.1 </h2>\n",
      "\n",
      "Here is the algorithm:\n",
      "\n",
      "<img src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/LU_blk_var5.png\" alt=\"Blocked LU factorization algorithm\" width=50%>\n",
      "    \n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Create the routine\n",
      "<code> LU_blk_var5( A ) </code>\n",
      "with the <a href=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/index.html\">Spark webpage</a> for the algorithm given above."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import flame\n",
      "import laff\n",
      "\n",
      "def LU_blk_var5(A, nb_alg):\n",
      "\n",
      "    ATL, ATR, \\\n",
      "    ABL, ABR  = flame.part_2x2(A, \\\n",
      "                               0, 0, 'TL')\n",
      "\n",
      "    while ATL.shape[0] < A.shape[0]:\n",
      "\n",
      "        block_size = min(ABR.shape[0], nb_alg)\n",
      "\n",
      "        A00, A01, A02, \\\n",
      "        A10, A11, A12, \\\n",
      "        A20, A21, A22  = flame.repart_2x2_to_3x3(ATL, ATR, \\\n",
      "                                                 ABL, ABR, \\\n",
      "                                                 block_size, block_size, 'BR')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        LU_unb_var5( A11 )\n",
      "        laff.trsm( 'Lower triangular', 'No transpose', 'Unit diagonal', A11, A12 )\n",
      "        laff.trsm( 'Upper triangular', 'Transpose', 'Nonunit diagonal', A11, A21 )\n",
      "        laff.gemm( -1.0, A21, A12, 1.0, A22 )\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        ATL, ATR, \\\n",
      "        ABL, ABR  = flame.cont_with_3x3_to_2x2(A00, A01, A02, \\\n",
      "                                               A10, A11, A12, \\\n",
      "                                               A20, A21, A22, \\\n",
      "                                               'TL')\n",
      "\n",
      "    flame.merge_2x2(ATL, ATR, \\\n",
      "                    ABL, ABR, A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h3> Test the routine </h3>\n",
      "\n",
      "Note that the code you generated using Spark has two input parameters, <code>A</code> and <code>nb_alg</code>. This <code>nb_alg</code> is the block size that you want to use to do your blocked LU decomposition, we'll set it arbitrarily to 20 for now and store it in the variable <code>nb</code>.\n",
      "\n",
      "<br>\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# recreate matrix A\n",
      "A = L * U\n",
      "\n",
      "# recreate the right-hand side\n",
      "b = A * xold\n",
      "\n",
      "# apply blocked LU to matrix A\n",
      "# remember nb holds our block size\n",
      "LU_blk_var5( A, nb )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Compare the overwritten matrix, $ A $, to the original matrices, $ L $ and $ U $.  The upper triangular part of $ A $ should equal $ U $ and the strictly lower triangular part of $ A $ should equal the strictly lower triangular part of $ L $. If this is the case, the maximum value in the matrix $A - L - U$ should be close to zero.\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compare A to the original L and U matrices\n",
      "print( 'Maximum value of (A - L - U) after factorization' )\n",
      "print( np.max( np.abs( A - np.tril(L,-1) - U ) ) ) #The \"-1\" ignores the diagonal"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> Implement the routine Solve from 6.3.4 using the blocked LU instead of the regular LU </h2>\n",
      "\n",
      "(if you have not yet visited Unit 6.3.4, do so now!)\n",
      "\n",
      "This time, we do NOT use Spark!  What we need to do is write a routine that, when given a matrix $ A $ and right-hand side vector $ b $, solves $ A x = b $, overwriting $ A $ with the LU factorization and overwriting $ b $ with the solution vector $ x $:\n",
      "\n",
      "<ul>\n",
      "<li>\n",
      "$ A \\rightarrow L U $, overwriting $ A $ with $ L $ and $ U $. Use the Blocked version.\n",
      "</li>\n",
      "<li>\n",
      "Solve $ L z = b $, overwriting $ b $ with $ z $.\n",
      "</li>\n",
      "<li>\n",
      "Solve $ U x = z $, where $ z $ is stored in vector $ b $ and $ x $ overwrites $ b $.\n",
      "</li>\n",
      "</ul>\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Create the routine\n",
      "<code> Solve( A, b ) </code>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Solve( A, b ):\n",
      "    \n",
      "    # insert appropriate calls to routines you have written here!\n",
      "    # remember the variable nb holds our block size\n",
      "    LU_blk_var5( A, 1 )\n",
      "    laff.trsv( 'Lower triangular', 'No transpose', 'Unit diagonal', A, b )\n",
      "    laff.trsv( 'Upper triangular', 'No transpose', 'Nonunit diagonal', A, b )\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h3> Test Solve </h3>\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# just to be sure, let's start over.  We'll recreate A, x, and b, run all the routines, and\n",
      "# then compare the updated b to the original vector x.\n",
      "\n",
      "A = L * U\n",
      "b = A * x\n",
      "\n",
      "Solve( A, b )\n",
      "\n",
      "print( '2-Norm of Updated b - original x' )\n",
      "print( np.linalg.norm(b - x) )\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In theory, <code> b - x </code> should yield a zero vector whose two-norm, $||b -x||_2$, is close to zero..."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h3> What if a new right-hand side comes along? </h3>\n",
      "\n",
      "What if we are presented with a new right-hand side, call it $ b_2 $, with which we want to solve $ A x = b_2 $, overwriting $ b_2 $ with the solution?  (We created such a $ b_2 $ at the top of this notebook.)\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>\n",
      "\n",
      "Notice that you can take the matrix $A $ that was modified by <code>Solve</code> and use it with the appropriate calls to <code>laff.trsv</code>:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# insert appropriate calls here.\n",
      "\n",
      "laff.trsv( 'Lower triangular', 'No transpose', 'Unit diagonal', A, b2 )\n",
      "laff.trsv( 'Upper triangular', 'No transpose', 'Nonunit diagonal', A, b2 )\n",
      "\n",
      "print( '2-Norm of updated b2 - original x2' )\n",
      "print( np.linalg.norm(b2 - x2) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$||x_2 - b_2||_2$ should be close to zero...\n",
      "\n",
      "\n",
      "<h2> <font color=red> Important: you should not refactor $ A $!!!! <font> </h2>"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}