05.3.2 Matrix-matrix multiplication by columns - Answer.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Matrix-matrix multiplication by columns"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We now 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_var1.png\" alt=\"Matrix-matrix multiplication by columns picture\" width=\"80%\">\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> The routine <code> Gemm_nn_unb_var1( A, B, C ) </code> </h2>\n",
      "\n",
      "This routine computes $ C := A B + C $ by columns.  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.gemv( trans, alpha, A, x, beta, y ) </code> which computes \n",
      "$ y := \\alpha A x + \\beta y $ or $ y := \\alpha A^T x + \\beta y $, depending on \n",
      "        parameter <code> trans</code>.  In particular, \n",
      "        <ul>\n",
      "        <li>\n",
      "        <code> laff.gemv( 'No transpose', alpha, A, x, beta, y ) </code> computes $ y := \\alpha A x + \\beta y $.\n",
      "            </li>\n",
      "        <li>\n",
      "        <code> laff.gemv( 'Transpose', alpha, A, x, beta, y ) </code> computes $ y := \\alpha A^T x + \\beta y $.\n",
      "            </li>\n",
      "            </ul>\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_var1(A, B, C):\n",
      "\n",
      "    BL, BR = flame.part_1x2(B, \\\n",
      "                            0, 'LEFT')\n",
      "\n",
      "    CL, CR = flame.part_1x2(C, \\\n",
      "                            0, 'LEFT')\n",
      "\n",
      "    while BL.shape[1] < B.shape[1]:\n",
      "\n",
      "        B0, b1, B2 = flame.repart_1x2_to_1x3(BL, BR, \\\n",
      "                                             1, 'RIGHT')\n",
      "\n",
      "        C0, c1, C2 = flame.repart_1x2_to_1x3(CL, CR, \\\n",
      "                                             1, 'RIGHT')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        laff.gemv( 'No transpose', 1.0, A, b1, 1.0, c1 )\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        BL, BR = flame.cont_with_1x3_to_1x2(B0, b1, B2, \\\n",
      "                                            'LEFT')\n",
      "\n",
      "        CL, CR = flame.cont_with_1x3_to_1x2(C0, c1, C2, \\\n",
      "                                            'LEFT')\n",
      "\n",
      "    flame.merge_1x2(CL, CR, C)\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_var1( 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": {}
  }
 ]
}