05.5.1 Multiplying upper triangular matrices - Answer.ipynb

{
 "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 <code> U </code> and <code> R </code>. The upper triangular part of matrix <code> C </code> 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": [
      "<h2>The algorithm  </h2>\n",
      "\n",
      "<image src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/5_5_1_10_Answer.png\" alt=\"Upper triangular matrix-matrix multiplication\" width=\"80%\">\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> The routine <code> Trtrmm_uu_unb_var1( U, R, C ) </code> </h2>\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 <code> U </code>, <code> R </code>, and <code> C </code> are not to be \"touched\".\n",
      "    \n",
      "The specific laff function you will want to use is some subset of\n",
      "<ul>\n",
      "<li> <code> laff.gemv( trans, alpha, A, x, beta, y ) </code> which computes $ y := \\alpha A x + \\beta y $ or $ y := \\alpha A^T x + \\beta y $ depending on whether <code> trans = 'No transpose' </code> or <code> trans = 'Transpose' </code> </li>\n",
      "    <li> <code> laff.ger( alpha, x, y, A ) </code> which computes the rank-1 update (adds a multiple of an outer product to a matrix)\n",
      "$ A := \\alpha x y^T + A $. </li>\n",
      "    <li> <code> laff.axpy( alpha, x, y ) </code></li>\n",
      "    <li>    <code> laff.dots( x, y, alpha ) </code></li>\n",
      "</ul>\n",
      "\n",
      "<font color=red> Hint:</font>  If you multiply with $ U_{00} $, you will want to use <code> np.triu( U00 ) </code> to make sure you don't compute with the nonzeroes below the diagonal.\n",
      "\n",
      "Use the <a href=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/index.html\"> Spark webpage</a> 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 <a href=\"http://edx-org-utaustinx.s3.amazonaws.com/UT501x/PictureFlame/PictureFLAME.html\"> PictureFLAME </a>, 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": {}
  }
 ]
}