{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Solving $ A x = b $ via LU factorization and triangular solves" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook, you will implement an 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 $.\n", "\n", " Be sure to make a copy!!!! \n", "\n", "
laff.dots( x, y, alpha )
$\\alpha := x^T y + \\alpha$\n",
"laff.invscal( alpha, x )
$x := x / \\alpha$ (See note below)\n",
"laff.axpy( alpha, x, y )
$y := \\alpha x + y$\n",
"laff.ger( alpha, x, y, A )
$A := \\alpha x y^T + A$\n",
"laff.invscal( alpha, x)
was added 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.\n",
"\n",
" LU_unb_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_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",
"\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" Ltrsv_unb_var1( L, b )
\n",
"with the Spark webpage for the algorithm given above."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import flame\n",
"\n",
"def Ltrsv_unb_var1(L, b):\n",
"\n",
" LTL, LTR, \\\n",
" LBL, LBR = flame.part_2x2(L, \\\n",
" 0, 0, 'TL')\n",
"\n",
" bT, \\\n",
" bB = flame.part_2x1(b, \\\n",
" 0, 'TOP')\n",
"\n",
" while LTL.shape[0] < L.shape[0]:\n",
"\n",
" L00, l01, L02, \\\n",
" l10t, lambda11, l12t, \\\n",
" L20, l21, L22 = flame.repart_2x2_to_3x3(LTL, LTR, \\\n",
" LBL, LBR, \\\n",
" 1, 1, 'BR')\n",
"\n",
" b0, \\\n",
" beta1, \\\n",
" b2 = flame.repart_2x1_to_3x1(bT, \\\n",
" bB, \\\n",
" 1, 'BOTTOM')\n",
"\n",
" #------------------------------------------------------------#\n",
"\n",
" laff.axpy( -beta1, l21, b2 ) # b2 := b2 - beta1*l21\n",
"\n",
" #------------------------------------------------------------#\n",
"\n",
" LTL, LTR, \\\n",
" LBL, LBR = flame.cont_with_3x3_to_2x2(L00, l01, L02, \\\n",
" l10t, lambda11, l12t, \\\n",
" L20, l21, L22, \\\n",
" 'TL')\n",
"\n",
" bT, \\\n",
" bB = flame.cont_with_3x1_to_2x1(b0, \\\n",
" beta1, \\\n",
" b2, \\\n",
" 'TOP')\n",
"\n",
" flame.merge_2x1(bT, \\\n",
" bB, b)\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"LU_unb_var5
, and use it to solve $ L z = b $, where $ L $ is unit lower triangular with its strictly lower triangular part stored in the strictly lower triangular part of $ A $.\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": [
"print( A )\n",
"print( b )\n",
"Ltrsv_unb_var1( A, b )\n",
"\n",
"print( 'updated b' )\n",
"print( b )"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To be able to perform the back substitution before we implement the upper triangular solve routine described in 6.3.3, we do it the hard way here. This allows us to see if the LU factorization followed by the solve with the unit lower triangular system gives the correct intermediate result.\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": [
"x[ 3 ] = b[ 3 ] / A[ 3,3 ]\n",
"x[ 2 ] = ( b[ 2 ] - A[ 2,3 ] * x[ 3 ] ) / A[ 2,2 ]\n",
"x[ 1 ] = ( b[ 1 ] - A[ 1,2 ] * x[ 2 ] - A[ 1,3 ] * x[ 3 ] ) / A[ 1,1 ]\n",
"x[ 0 ] = ( b[ 0 ] - A[ 0,1 ] * x[ 1 ] - A[ 0,2 ] * x[ 2 ]- A[ 0,3 ] * x[ 3 ] ) / A[ 0,0 ]"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print( 'x = ' )\n",
"print( x )\n",
"\n",
"print( 'x - xold' )\n",
"print( x - xold )"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" x - xold
should yield a zero vector"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Utrsv_unb_var1( U, b )
\n",
"with the Spark webpage for the algorithm given above.\n",
"\n",
" Hint: Implement $ \\beta_1 := \\beta - u_{12}^T b_2 $ as laff.dots( -u12t, b2, beta1 )
"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import flame\n",
"import laff as laff\n",
"\n",
"def Utrsv_unb_var1(U, b):\n",
"\n",
" UTL, UTR, \\\n",
" UBL, UBR = flame.part_2x2(U, \\\n",
" 0, 0, 'BR')\n",
"\n",
" bT, \\\n",
" bB = flame.part_2x1(b, \\\n",
" 0, 'BOTTOM')\n",
"\n",
" while UBR.shape[0] < U.shape[0]:\n",
"\n",
" U00, u01, U02, \\\n",
" u10t, upsilon11, u12t, \\\n",
" U20, u21, U22 = flame.repart_2x2_to_3x3(UTL, UTR, \\\n",
" UBL, UBR, \\\n",
" 1, 1, 'TL')\n",
"\n",
" b0, \\\n",
" beta1, \\\n",
" b2 = flame.repart_2x1_to_3x1(bT, \\\n",
" bB, \\\n",
" 1, 'TOP')\n",
"\n",
" #------------------------------------------------------------#\n",
" \n",
" laff.dots( -u12t, b2, beta1 ) # beta1 := beta1 - u21 * b2\n",
" laff.invscal( upsilon11, beta1 ) # beta1 := beta1/upsilon11\n",
"\n",
" #------------------------------------------------------------#\n",
"\n",
" UTL, UTR, \\\n",
" UBL, UBR = flame.cont_with_3x3_to_2x2(U00, u01, U02, \\\n",
" u10t, upsilon11, u12t, \\\n",
" U20, u21, U22, \\\n",
" 'BR')\n",
"\n",
" bT, \\\n",
" bB = flame.cont_with_3x1_to_2x1(b0, \\\n",
" beta1, \\\n",
" b2, \\\n",
" 'BOTTOM')\n",
"\n",
" flame.merge_2x1(bT, \\\n",
" bB, b)\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"LU_unb_var5
and Ltrsv_unb_var1
and use it to solve $ U x = b $, where $ U $ is upper triangular and stored in the upper triangular part of $ A $.\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",
"LU_unb_var5( A )\n",
"Ltrsv_unb_var1( A, b )\n",
"Utrsv_unb_var1( A, b )\n",
"\n",
"print( 'updated b' )\n",
"print( b )\n",
"print( 'original x' )\n",
"print( x )\n",
"print( 'b - x' )\n",
"print( b - x )\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In theory, b - x
should yield a zero vector..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" 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",
" LU_unb_var5( A )\n",
" Ltrsv_unb_var1( A, b )\n",
" Utrsv_unb_var1( A, b )\n",
" \n",
" "
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" b - x
should yield a zero vector..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solve
and use it with Ltrsv_unb_var1
and Utrsv_unb_var1
:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# insert appropriate calls here.\n",
"\n",
"Ltrsv_unb_var1( A, y )\n",
"Utrsv_unb_var1( A, y )\n",
"\n",
"print( 'y = ' )\n",
"print( y )\n",
"\n",
"print( 'x2 - y =' )\n",
"print( x2 - y )"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"x2 - y should evaluate to the zero vector.\n",
"\n",
"\n",
"