05.3.1 Lots of loops.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Lots of loops"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook illustrates the different ways in which loops for matrix-matrix multiplication can be ordered.  Let's start by creating 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 basic algorithm</h2\n",
      "\n",
      "<p>\n",
      "Given $ A \\in \\mathbb{R}^{m \\times k} $, $ B \\in \\mathbb{R}^{k \\times n} $, and $ C \\in \\mathbb{R}^{m \\times n} $, we will consider $ C := A B + C $.\n",
      "</p>\n",
      "    \n",
      "<p>\n",
      "    Now, recall that the $ i,j $ element of $ A B $ is computed as the dot product of \n",
      "the $ i $th row of $ A $ with the $ j $th column of $ B $:\n",
      "</p>\n",
      "\n",
      "<p>\n",
      "    $\\sum_{p=0}^{k-1} \\alpha_{i,j} \\beta_{i,j}$\n",
      "</p>\n",
      "\n",
      "<p>\n",
      "    and here, by adding to $ C $ we get\n",
      "</p>\n",
      "\n",
      "<p>\n",
      "$ \\gamma_{i,j} = \\sum_{p=0}^{k-1} \\alpha_{i,j} \\beta_{i,j} + \\gamma_{i,j}.$\n",
      "</p>\n",
      "\n",
      "<p>\n",
      "    Now, we have to loop over all elements of $ C $.  The code, without the FLAMEpy API, becomes\n",
      "</p>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def MMmult_lots_of_loops( A, B, C ):\n",
      "\n",
      "    m, n = np.shape( C )\n",
      "    m, k = np.shape( A )\n",
      "    \n",
      "    # i,j,p\n",
      "    for i in range( m ):                     \n",
      "        for j in range( n ):                    \n",
      "            for p in range( k ):                    \n",
      "                C[ i,j ] = A[ i,p ] * B[ p, j ] + C[ i,j ]\n",
      "    \n",
      "    # Comment out the loops above and do the other 5 implementations here one at a time..."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C = np.matrix( np.copy( Cold ) )             # restore C\n",
      "\n",
      "MMmult_lots_of_loops( 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": [
      "Now, go back and systematically move the loops around, so that in the end you try out all six orders of the loops: three choices for the first, outermost, loop; two choices for the secod loop; one choice for the third loop, for a total of $ 3! $ (3 factorial) choices. Check that you get the right answer, regardless. \n",
      "\n",
      "(We suggest you just change the box in which the routine is defined and comment out variations that you've already tested.  Be careful with indentation.)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Why $ C := A B + C $ rather than $ C := A B $?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice that we could have written a routine to compute $ C := A B $ instead, given below."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def MMmult_C_eq_AB( A, B, C ):\n",
      "\n",
      "    m, n = np.shape( C )\n",
      "    m, k = np.shape( A )\n",
      "    \n",
      "    for i in range( m ):                     \n",
      "        for j in range( n ):   \n",
      "            C[ i,j ] = 0.0\n",
      "            for p in range( k ):                    \n",
      "                C[ i,j ] = A[ i,p ] * B[ p, j ] + C[ i,j ]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C = np.matrix( np.copy( Cold ) )             # restore C\n",
      "\n",
      "MMmult_C_eq_AB( A, B, C )\n",
      "\n",
      "print( 'C - ( A * B )' )\n",
      "print( C - ( A * B ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, start changing the order of the loops.  You notice it is not quite as simple.  But, if you have a routine for computing $ C := A B + C $, you can always initialize $ C = 0 $ (the zero matrix) and then use it to call $ C := A B $:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C = np.matrix( np.zeros( np.shape( C ) ) )          # initialize C = 0 \n",
      "\n",
      "MMmult_lots_of_loops( A, B, C )\n",
      "\n",
      "print( 'C - ( A * B )' )\n",
      "print( C - ( A * B ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}