{ "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", "
laff.trsv( uplo, trans, diag, A, b )
Solves $Ax = b$ where $x$ and $b$ are vectors.\n",
"laff.trsm( uplo, trans, diag, A, B )
Solves $AX = B$ where $X$ and $B$ are matrices.\n",
"laff.gemm( alpha, A, B, beta, C )
$C := \\alpha A B + \\beta C$\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": [
"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",
" 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": [
" b - x
should yield a zero vector whose two-norm, $||b -x||_2$, is close to zero..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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",
"