01.5.4 Implementing an axpy routine - Answer.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Implementing an axpy routine"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Preliminaries"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<p><font color=red> Again, copy this notebook so that you don't corrupt the original!  Then you can \"play\" with the copy of the notebook all you want! </font> </p>\n",
      "\n",
      "<p>\n",
      "<b>\n",
      "NOTE: A common problem that students have with IPython notebooks is not understanding that when the code in the gray boxes (cells) is executed, it assigns variables that persist the whole time that the notebook is open. Further, some cells rely on variables assigned by earlier cells.  If you execute these cells out of order, or if you execute the same cell twice, then you may end up changing the value of the variables.  To correct this, click on \"Cell\" at the top and execute \"run all above\" or \"run all\".  You can also reset all cells by clicking \"Cell -> All Output -> Clear\"\n",
      "</b>\n",
      "</p>\n",
      "\n",
      "<p>In this notebook, you are asked to write the loop that implements an axpy operation.</p>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's start by importing numpy and creating vectors $ x = \\left( \\begin{array}{r} 1 \\\\ 2 \\\\ 3 \\end{array} \\right) $, $ y = \\left( \\begin{array}{r} -1 \\\\ 0 \\\\ -2 \\end{array} \\right) $ and a scalar $ \\alpha = 2.5 $.  \n",
      "\n",
      "Execute the code in the box by clicking in the box and then on \"Cell -> Run\".  Alternative, click on the box and push \"Shift\" and \"Return\" (or \"Enter\") together."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np    # This imports a package called \"numpy\" that will make working with matrices \n",
      "                      # simpler\n",
      "\n",
      "# create two two-dimensional matrices of only one column each.  \n",
      "x = np.matrix( '1.;2.;3.' )\n",
      "print( 'x = ' )\n",
      "print( x )\n",
      "\n",
      "y = np.matrix( '-1.;0.;-2.' )\n",
      "print( 'y = ' )\n",
      "print( y )\n",
      "\n",
      "alpha = 2.5\n",
      "print( 'alpha = ' )\n",
      "print( alpha )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "With Python, you can simply compute $ y := \\alpha x + y $.  We will first make a copy of $ y $ so we can easily reset $ y $ to its original contents in the future.  For this, we will use the laff.copy routine."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import laff\n",
      "\n",
      "yold = np.matrix( '0;0;0' )\n",
      "\n",
      "laff.copy( y, yold )\n",
      "\n",
      "print( 'yold = ' )\n",
      "print( y )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, let's update $ y := \\alpha x + y $:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y = yold\n",
      "\n",
      "print( 'y before axpy:')\n",
      "\n",
      "print( y )\n",
      "\n",
      "y = alpha * x + y\n",
      "\n",
      "print( 'y after axpy: ' )\n",
      "print( y )\n",
      "\n",
      "print( 'compare new y to  alpha * x + yold:' )\n",
      "print( y - ( alpha * x + yold ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Computing an axpy with a loop"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, we want you to write a loop that updates the three entries in $ y $ with $ \\alpha x + y $.  \n",
      "Recall that the loop for copying three entries of a vector $ x $ into a vector $ y $ was given by \n",
      "\n",
      "<code>\n",
      "        for i in range( 3 ):\n",
      "           y[ i, 0 ] = x[ i, 0 ] \n",
      "</code>\n",
      "\n",
      "Don't forget the indentation before <code> y[ i, 0 ] </code> since that indicates what statements are part of the \"loop body\" (in other words, what statements are executed as part of the \"for loop\")."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "laff.copy( yold, y )\n",
      "\n",
      "print( 'y before axpy:')\n",
      "\n",
      "print( y )\n",
      "\n",
      "for i in range( 3 ):\n",
      "    y[ i, 0 ] = alpha * x[ i, 0 ] + y[ i, 0 ]\n",
      "\n",
      "print( 'y after axpy: ' )\n",
      "print( y )\n",
      "\n",
      "print( 'compare new y to  alpha * x + yold:' )\n",
      "print( y - ( alpha * x + yold ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The result should be:\n",
      "\n",
      "<code>\n",
      "y before axpy:\n",
      "[[-1.]\n",
      " [ 0.]\n",
      " [-2.]]\n",
      "y after axpy: \n",
      "[[ 1.5]\n",
      " [ 5. ]\n",
      " [ 5.5]]\n",
      "compare new y to  alpha * x + yold:\n",
      "[[ 0.]\n",
      " [ 0.]\n",
      " [ 0.]]\n",
      "</code>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Axpy as a simple routine"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Writing the \"for loop\" every time you want to compute $ y := \\alpha * x + y $ is labor intensive and unnecessary.  Obviously, you can do \"<code> y = alpha * x + y </code>\", but the point of this exercise is for you to write your own routine.  For this reason, you are going to write a routine, axpy( alpha, x, y ).\n",
      "\n",
      "Recall what it means to update.  $ y $ becomes $ \\alpha x + y $ when $ x $ and $ y $ are both of size $ m $:\n",
      "$$\n",
      "\\left( \\begin{array}{c}\n",
      "\\psi_0 \\\\\n",
      "\\psi_1 \\\\\n",
      "\\vdots \\\\\n",
      "\\psi_{m-1}\n",
      "\\end{array}\n",
      "\\right)\n",
      ":=\n",
      "\\left( \\begin{array}{c}\n",
      "\\alpha \\chi_0 \\\\\n",
      "\\alpha \\chi_1 \\\\\n",
      "\\vdots \\\\\n",
      "\\alpha \\chi_{m-1}\n",
      "\\end{array}\n",
      "\\right)\n",
      "+\n",
      "\\left( \\begin{array}{c}\n",
      "\\psi_0 \\\\\n",
      "\\psi_1 \\\\\n",
      "\\vdots \\\\\n",
      "\\psi_{m-1}\n",
      "\\end{array}\n",
      "\\right)\n",
      "=\n",
      "\\left( \\begin{array}{c}\n",
      "\\alpha \\chi_0 + \\psi_0 \\\\\n",
      "\\alpha \\chi_1 + \\psi_1 \\\\\n",
      "\\vdots \\\\\n",
      "\\alpha \\chi_{m-1} + \\psi_{m-1}\n",
      "\\end{array}\n",
      "\\right)\n",
      "$$\n",
      "\n",
      "As an algorithm this can be written as\n",
      "<p>\n",
      "<b> for </b> $ i = 0, \\ldots , m-1 $  <br>\n",
      " $ ~~~ \\psi_i := \\alpha \\chi_i + \\psi_i $ <br>\n",
      "<b> endfor </b>\n",
      "</p>\n",
      "\n",
      "<p>\n",
      "Complete the following routine to implement this:\n",
      "</p>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def axpy( alpha, x, y ):\n",
      "\n",
      "    m, n = np.shape( x )\n",
      "    \n",
      "    for i in range( m ):\n",
      "        y[ i, 0 ] = alpha * x[ i, 0 ] + y[ i, 0 ]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Be sure the run the box above, or this notebook won't know about the routine!!!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, if you execute"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "laff.copy( yold, y )\n",
      "\n",
      "print( 'y before axpy:')\n",
      "\n",
      "print( y )\n",
      "\n",
      "axpy( alpha, x, y )\n",
      "\n",
      "print( 'y after axpy: ' )\n",
      "print( y )\n",
      "\n",
      "print( 'compare new y to  alpha * x + yold:' )\n",
      "print( y - ( alpha * x + yold ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The result should be:\n",
      "\n",
      "<code>\n",
      "y before axpy:\n",
      "[[-1.]\n",
      " [ 0.]\n",
      " [-2.]]\n",
      "y after axpy: \n",
      "[[ 1.5]\n",
      " [ 5. ]\n",
      " [ 5.5]]\n",
      "compare new y to  alpha * x + yold:\n",
      "[[ 0.]\n",
      " [ 0.]\n",
      " [ 0.]]\n",
      "</code>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "A complete axpy function as part of the LAFF library"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As we proceed with routines that become progressively more advanced, we are going to need a general axpy routine where $ x $ and $ y $ can be row and/or column vectors.  \n",
      "\n",
      "This routine is part of the 'laff' library.  If you do\n",
      "\n",
      "<code>\n",
      "import laff\n",
      "</code>\n",
      "\n",
      "then <code> laff.axpy( alpha, x, y ) </code> will perform the desired axpy, when <code> x </code> and <code> y </code> are column and/or a row vectors.  If you really want to see what this routine looks like, then ask for it on the discussion forum and we'll point you to where it can be found."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import laff\n",
      "\n",
      "laff.copy( yold, y )\n",
      "\n",
      "print( 'y before axpy:')\n",
      "\n",
      "print( y )\n",
      "\n",
      "laff.axpy( alpha, x, y )\n",
      "\n",
      "print( 'y after axpy: ' )\n",
      "print( y )\n",
      "\n",
      "print( 'compare new y to  alpha * x + yold:' )\n",
      "print( y - ( alpha * x + yold ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Need a challenge?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In \"1.5.2 Implementing a copy routine\", we gave a complete routine that implements a copy from a row/column vector to a row/column vector, checked whether the parameters were legal, and had comments in it.  If you feel up to the challenge, below write a similar routine <code> axpy </code> that works for row and column vectors, checks the parameters, and has comments.  Be sure to test your implementation.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def axpy( alpha, x, y ):\n",
      "    ### You fill in the rest!"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}