05.3.4 Matrix-matrix multiplication via rank-1 updates - Answer.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Matrix-matrix multiplication via rank-1 updates"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We continue to look at how the FLAMEPy API can be used to implement different matrix-matrix multiplication algorithms.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First, we create some matrices."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "m = 4\n",
      "n = 3\n",
      "k = 5\n",
      "\n",
      "C = np.matrix( np.random.random( (m, n) ) )\n",
      "print( 'C = ' )\n",
      "print( C )\n",
      "\n",
      "Cold = np.matrix( np.zeros( (m,n ) ) )\n",
      "Cold = np.matrix( np.copy( C ) )           # an alternative way of doing a \"hard\" copy, in this case of a matrix\n",
      "    \n",
      "A = np.matrix( np.random.random( (m, k) ) )\n",
      "print( 'A = ' )\n",
      "print( A )\n",
      "\n",
      "B = np.matrix( np.random.random( (k, n) ) )\n",
      "print( 'B = ' )\n",
      "print( B )"
     ],
     "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/Gemm_nn_unb_var3.png\" alt=\"Matrix-matrix multiplication via rank-1 updates picture\" width=\"80%\">\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> The routine <code> Gemm_nn_unb_var3( A, B, C ) </code> </h2>\n",
      "\n",
      "This routine computes $ C := A B + C $ via rank-1 updates.  The \"\\_nn\\_\" means that this is the \"No transpose, No transpose\" case of matrix multiplication.  \n",
      "The reason for this is that the operations $ C := A^T B + C $ (\"\\_tn\\_\" or \"Transpose, No transpose\"), $ C := A B^T + C $ (\"\\_nt\\_\" or \"No transpose, Transpose\"), and $ C := A^T B^T + C $ (\"\\_tt\\_\" or \"Transpose, Transpose\") are also encountered.  \n",
      "    \n",
      "The specific laff function we will use is\n",
      "<ul>\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 $. \n",
      "    </li>\n",
      "</ul>\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 flame\n",
      "import laff as laff\n",
      "\n",
      "def Gemm_nn_unb_var3(A, B, C):\n",
      "\n",
      "    AL, AR = flame.part_1x2(A, \\\n",
      "                            0, 'LEFT')\n",
      "\n",
      "    BT, \\\n",
      "    BB  = flame.part_2x1(B, \\\n",
      "                         0, 'TOP')\n",
      "\n",
      "    while AL.shape[1] < A.shape[1]:\n",
      "\n",
      "        A0, a1, A2 = flame.repart_1x2_to_1x3(AL, AR, \\\n",
      "                                             1, 'RIGHT')\n",
      "\n",
      "        B0,  \\\n",
      "        b1t, \\\n",
      "        B2   = flame.repart_2x1_to_3x1(BT, \\\n",
      "                                       BB, \\\n",
      "                                       1, 'BOTTOM')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        laff.ger( 1.0, a1, b1t, C )\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        AL, AR = flame.cont_with_1x3_to_1x2(A0, a1, A2, \\\n",
      "                                            'LEFT')\n",
      "\n",
      "        BT, \\\n",
      "        BB  = flame.cont_with_3x1_to_2x1(B0,  \\\n",
      "                                         b1t, \\\n",
      "                                         B2,  \\\n",
      "                                         'TOP')\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C = np.matrix( np.copy( Cold ) )               # restore C \n",
      "\n",
      "Gemm_nn_unb_var3( A, B, C )\n",
      "\n",
      "print( 'C - ( Cold + A * B )' )\n",
      "print( C - ( Cold + A * B ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Bingo! It works!"
     ]
    },
    {
     "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": {}
  }
 ]
}