{
"metadata": {
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 1,
"metadata": {},
"source": [
"Multiplying upper triangular matrices"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook helps you implement the operation $ C := U R $ where $ C, U, R \\in \\mathbb{R}^{n \\times n} $, and $ U $ and $ R $ are upper triangular. $ U $ and $ R $ are stored in the upper triangular part of numpy matrices U and R . The upper triangular part of matrix C is to be overwritten with the resulting upper triangular matrix. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, we create some matrices."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import numpy as np\n",
"\n",
"n = 5\n",
"\n",
"C = np.matrix( np.random.random( (n, n) ) )\n",
"print( 'C = ' )\n",
"print( C )\n",
"\n",
"Cold = np.matrix( np.zeros( (n,n ) ) )\n",
"Cold = np.matrix( np.copy( C ) ) # an alternative way of doing a \"hard\" copy, in this case of a matrix\n",
" \n",
"U = np.matrix( np.random.random( (n, n) ) )\n",
"print( 'U = ' )\n",
"print( U )\n",
"\n",
"R = np.matrix( np.random.random( (n, n) ) )\n",
"print( 'R = ' )\n",
"print( R )"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"The algorithm
\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" The routine Trtrmm_uu_unb_var1( U, R, C )
\n",
"\n",
"This routine computes $ C := U R + C $. The \"\\_uu\\_\" means that $ U $ and $ R $ are upper triangular matrices (which means $ C $ is too). However, the lower triangular parts of numpy matrices U , R , and C are not to be \"touched\".\n",
" \n",
"The specific laff function you will want to use is some subset of\n",
"\n",
"-
laff.gemv( trans, alpha, A, x, beta, y ) which computes $ y := \\alpha A x + \\beta y $ or $ y := \\alpha A^T x + \\beta y $ depending on whether trans = 'No transpose' or trans = 'Transpose' \n",
" -
laff.ger( alpha, x, y, A ) which computes the rank-1 update (adds a multiple of an outer product to a matrix)\n",
"$ A := \\alpha x y^T + A $. \n",
" -
laff.axpy( alpha, x, y ) \n",
" -
laff.dots( x, y, alpha ) \n",
"
\n",
"\n",
" Hint: If you multiply with $ U_{00} $, you will want to use np.triu( U00 ) to make sure you don't compute with the nonzeroes below the diagonal.\n",
"\n",
"Use the Spark webpage to generate a code skeleton. (Make sure you adjust the name of the routine.)"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import laff as laff\n",
"import flame\n",
"\n",
"def Trtrmm_uu_unb_var1(U, R, C):\n",
"\n",
" UTL, UTR, \\\n",
" UBL, UBR = flame.part_2x2(U, \\\n",
" 0, 0, 'TL')\n",
"\n",
" RTL, RTR, \\\n",
" RBL, RBR = flame.part_2x2(R, \\\n",
" 0, 0, 'TL')\n",
"\n",
" CTL, CTR, \\\n",
" CBL, CBR = flame.part_2x2(C, \\\n",
" 0, 0, 'TL')\n",
"\n",
" while UTL.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, 'BR')\n",
"\n",
" R00, r01, R02, \\\n",
" r10t, rho11, r12t, \\\n",
" R20, r21, R22 = flame.repart_2x2_to_3x3(RTL, RTR, \\\n",
" RBL, RBR, \\\n",
" 1, 1, 'BR')\n",
"\n",
" C00, c01, C02, \\\n",
" c10t, gamma11, c12t, \\\n",
" C20, c21, C22 = flame.repart_2x2_to_3x3(CTL, CTR, \\\n",
" CBL, CBR, \\\n",
" 1, 1, 'BR')\n",
"\n",
" #------------------------------------------------------------#\n",
"\n",
" # imporant: You do not need to recompute C00 = U00 * R00!!!!\n",
" laff.gemv( 'No transpose', 1.0, np.triu( U00 ), r01, 0.0, c01 )\n",
" laff.axpy( rho11, u01, c01 )\n",
" laff.dot( rho11, upsilon11, gamma11 )\n",
"\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",
" 'TL')\n",
"\n",
" RTL, RTR, \\\n",
" RBL, RBR = flame.cont_with_3x3_to_2x2(R00, r01, R02, \\\n",
" r10t, rho11, r12t, \\\n",
" R20, r21, R22, \\\n",
" 'TL')\n",
"\n",
" CTL, CTR, \\\n",
" CBL, CBR = flame.cont_with_3x3_to_2x2(C00, c01, C02, \\\n",
" c10t, gamma11, c12t, \\\n",
" C20, c21, C22, \\\n",
" 'TL')\n",
"\n",
" flame.merge_2x2(CTL, CTR, \\\n",
" CBL, CBR, C)\n",
"\n",
"\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"C = np.matrix( np.copy( Cold ) ) # restore C \n",
"\n",
"Trtrmm_uu_unb_var1( U, R, C )\n",
"\n",
"# compute it using numpy *. This is a little complex, since we want to make sure we\n",
"# don't change the lower triangular part of C\n",
"# Cref = np.tril( Cold, -1 ) + np.triu( U ) * np.triu( R )\n",
"Cref = np.triu( U ) * np.triu( R ) + np.tril( Cold, -1 )\n",
"print( 'C - Cref' )\n",
"print( C - Cref )"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In theory, the result matrix should be (approximately) zero."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Watch the algorithm at work!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Copy and paste the code into PictureFLAME , a webpage where you can watch your routine in action. Just cut and paste into the box. \n",
"\n",
"Disclaimer: we implemented a VERY simple interpreter. If you do something wrong, we cannot guarantee the results. But if you do it right, you are in for a treat.\n",
"\n",
"If you want to reset the problem, just click in the box into which you pasted the code and hit \"next\" again."
]
}
],
"metadata": {}
}
]
}