04.3.3.4 Symmetric Matrix Vector Multiply Routines Challenge Question.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Symmetric Matrix Vector Multiply Routines (stored in lower triangle)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook walks you through how to implement $ y := A x + y $ where $ A $ is symmetric.  The challenge is to only access the matrix by columns."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Getting started"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We will use some functions that are part of our laff library (of which this function will become a part) as well as some routines from the FLAME API (Application Programming Interface) that allows us to write code that closely resembles how we typeset algorithms using the FLAME notation.  These functions are imported with the \"import laff as laff\" and \"import flame\" statements."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Algorithm for the case where the matrix is stored in upper triangular part of $ A $"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The routine <br> <code> Symv_u_unb_var3( A, x, y ) </code>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This routine, given symmetric matrix $ A \\in \\mathbb{R}^{n \\times n} $, $ x \\in \\mathbb{R}^n $, and $ y \\in \\mathbb{R}^n $, computes $ y := A x + y $.  The \"_u_\" in the name of the routine indicates that $ A $ is stored in the upper triangular part of the matrix.\n",
      "\n",
      "The specific laff functions you may want to consider using are\n",
      "<ul>\n",
      "<li> <code> laff.dots( x, y, alpha ) </code> which computes $ \\alpha := x^T y + \\alpha $.  </li>\n",
      "<li> <code> laff.axpy( alpha, x, y ) </code> which computes $ y := \\alpha x + y $.  </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": [
      "#Insert code here..."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Testing"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's quickly test the routine by creating a 4 x 4 matrix and related vectors, performing the computation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import random\n",
      "from numpy import matrix\n",
      "\n",
      "A = matrix( random.rand( 4,4 ) )\n",
      "x = matrix( random.rand( 4,1 ) )\n",
      "y = matrix( random.rand( 4,1 ) )\n",
      "yold = matrix( random.rand( 4,1 ) )\n",
      "\n",
      "print( 'A before =' )\n",
      "print( A )\n",
      "\n",
      "print( 'x before =' )\n",
      "print( x )\n",
      "\n",
      "print( 'y before =' )\n",
      "print( y )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "# Notice that A is not symmetric.  We now \"symmetrize it\"\n",
      "Asymm = np.triu( A ) + np.transpose( np.triu( A, 1 ) )\n",
      "\n",
      "laff.copy( y, yold )   # save the original vector y\n",
      "\n",
      "Symv_u_unb_var3( A, x, y )\n",
      "\n",
      "print( 'y after =' )\n",
      "print( y )\n",
      "\n",
      "print( 'y - ( Asymm * x + yold ) = ' )\n",
      "print( y - ( Asymm * x + yold ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Bingo, it seems to work!  (Notice that we are doing floating point computations, which means that due to rounding you may not get an exact \"0\".)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Watch your code in action!"
     ]
    },
    {
     "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."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Algorithm for the case where the matrix is stored in lower triangular part of $ A $"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The routine <br> <code> Symv_l_unb_var3( A, x, y ) </code>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This routine, given symmetric matrix $ A \\in \\mathbb{R}^{n \\times n} $, $ x \\in \\mathbb{R}^n $, and $ y \\in \\mathbb{R}^n $, computes $ y := A x + y $.  The \"_l_\" in the name of the routine indicates that $ A $ is stored in the lower triangular part of the matrix.\n",
      "\n",
      "The specific laff functions you may want to consider using are\n",
      "<ul>\n",
      "<li> <code> laff.dots( x, y, alpha ) </code> which computes $ \\alpha := x^T y + \\alpha $.  </li>\n",
      "<li> <code> laff.axpy( alpha, x, y ) </code> which computes $ y := \\alpha x + y $.  </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": [
      "#Insert code here..."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Testing"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's quickly test the routine by creating a 4 x 4 matrix and related vectors, performing the computation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import random\n",
      "from numpy import matrix\n",
      "\n",
      "A = matrix( random.rand( 4,4 ) )\n",
      "x = matrix( random.rand( 4,1 ) )\n",
      "y = matrix( random.rand( 4,1 ) )\n",
      "yold = matrix( random.rand( 4,1 ) )\n",
      "\n",
      "print( 'A before =' )\n",
      "print( A )\n",
      "\n",
      "print( 'x before =' )\n",
      "print( x )\n",
      "\n",
      "print( 'y before =' )\n",
      "print( y )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "# Notice that A is not symmetric.  We now \"symmetrize it\"\n",
      "Asymm = np.tril( A ) + np.transpose( np.tril( A, -1 ) )\n",
      "\n",
      "laff.copy( y, yold )   # save the original vector y\n",
      "\n",
      "Symv_l_unb_var3( A, x, y )\n",
      "\n",
      "print( 'y after =' )\n",
      "print( y )\n",
      "\n",
      "print( 'y - ( Asymm * x + yold ) = ' )\n",
      "print( y - ( Asymm * x + yold ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Bingo, it seems to work!  (Notice that we are doing floating point computations, which means that due to rounding you may not get an exact \"0\".)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Watch your code in action!"
     ]
    },
    {
     "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."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}