{ "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", " Be sure to make a copy!!!! \n", "\n", "

Preliminaries

\n", "\n", "Here is a list of laff routines that you might want to use in this notebook:\n", "\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", " Make sure you call the routine LU_unb_var5 \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": [ "

Now, let's create a matrix $ A $ and right-hand side $ b $

" ] }, { "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": [ "

Implement the blocked LU factorization routine from 6.4.1

\n", "\n", "Here is the algorithm:\n", "\n", "\"Blocked\n", " \n", " 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. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the routine\n", " LU_blk_var5( A ) \n", "with the Spark webpage 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": [ "

Test the routine

\n", "\n", "Note that the code you generated using Spark has two input parameters, A and nb_alg. This nb_alg 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 nb.\n", "\n", "
\n", "\n", " 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. " ] }, { "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", " 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. \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": [ "

Implement the routine Solve from 6.3.4 using the blocked LU instead of the regular LU

\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", "\n", "\n", " 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. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the routine\n", " Solve( A, b ) " ] }, { "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": [ "

Test Solve

\n", "\n", " 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. " ] }, { "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, b - x should yield a zero vector whose two-norm, $||b -x||_2$, is close to zero..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

What if a new right-hand side comes along?

\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", " 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. \n", "\n", "Notice that you can take the matrix $A $ that was modified by Solve and use it with the appropriate calls to laff.trsv:" ] }, { "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", "

Important: you should not refactor $ A $!!!!

" ] } ], "metadata": {} } ] }