04.3.2.7 Lower Triangular Matrix Vector Multiply Routines (overwriting x) - Answer.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Lower Triangular Matrix Vector Multiply Routines"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook walks you through how to implement $ x := L x $ where $ L $ is lower triangular.  Vector $ y $ is not to be touched (and, indeed, not even passed into the routines).  This is a little trickier than you might think.  Order matters, since you don't want to corrupt values of $ x $ that are still needed.  Indeed, you may want to do a few small problems by hand if you don't get the right answer.  Also, PictureFLAME may help you see what is going on."
     ]
    },
    {
     "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 that takes dot products"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The routine <br> <code> Trmv_ln_unb_var1( L, x ) </code>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This routine, given lower triangular $ L \\in \\mathbb{R}^{n \\times n} $ and $ x \\in \\mathbb{R}^n $ computes $ x := L x $.  The \"_un_\" in the name of the routine indicates this is the \"lower, no transpose\" matrix-vector multiplication.  \n",
      "\n",
      "The specific laff functions we will use are \n",
      "<ul>\n",
      "<li> <code> laff.dots( x, y, alpha ) </code> which computes $ \\alpha := x^T y + \\alpha $.  </li>\n",
      "<li> <code> laff.scal( alpha, x ) </code> which computes $ x := \\alpha x $.  </li> (You will want to use this one to update $ \\chi_1 $ if you want to use PictureFrame.)\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": "markdown",
     "metadata": {},
     "source": [
      "<h1> <font color=red> NOTE: The algorithms for this operation march from bottom-right to top-left!!! </font> </h1>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import flame\n",
      "import laff as laff\n",
      "\n",
      "def Trmv_ln_unb_var1(L, x):\n",
      "\n",
      "    LTL, LTR, \\\n",
      "    LBL, LBR  = flame.part_2x2(L, \\\n",
      "                               0, 0, 'BR')\n",
      "\n",
      "    xT, \\\n",
      "    xB  = flame.part_2x1(x, \\\n",
      "                         0, 'BOTTOM')\n",
      "\n",
      "    while LBR.shape[0] < L.shape[0]:\n",
      "\n",
      "        L00,  l01,      L02,  \\\n",
      "        l10t, lambda11, l12t, \\\n",
      "        L20,  l21,      L22   = flame.repart_2x2_to_3x3(LTL, LTR, \\\n",
      "                                                        LBL, LBR, \\\n",
      "                                                        1, 1, 'TL')\n",
      "\n",
      "        x0,   \\\n",
      "        chi1, \\\n",
      "        x2    = flame.repart_2x1_to_3x1(xT, \\\n",
      "                                        xB, \\\n",
      "                                        1, 'TOP')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        laff.scal( lambda11, chi1 )\n",
      "        laff.dots( l10t, x0, chi1 )\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        LTL, LTR, \\\n",
      "        LBL, LBR  = flame.cont_with_3x3_to_2x2(L00,  l01,      L02,  \\\n",
      "                                               l10t, lambda11, l12t, \\\n",
      "                                               L20,  l21,      L22,  \\\n",
      "                                               'BR')\n",
      "\n",
      "        xT, \\\n",
      "        xB  = flame.cont_with_3x1_to_2x1(x0,   \\\n",
      "                                         chi1, \\\n",
      "                                         x2,   \\\n",
      "                                         'BOTTOM')\n",
      "\n",
      "    flame.merge_2x1(xT, \\\n",
      "                    xB, x)\n",
      "\n"
     ],
     "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",
      "L = matrix( random.rand( 4,4 ) )\n",
      "x = matrix( random.rand( 4,1 ) )\n",
      "xold = matrix( random.rand( 4,1 ) )\n",
      "\n",
      "# Notice that L is not lower triangular.  We will only use the lower triangular part.\n",
      "\n",
      "print( 'L before =' )\n",
      "print( L )\n",
      "\n",
      "print( 'x before =' )\n",
      "print( x )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "laff.copy( x, xold )   # save the original vector x\n",
      "\n",
      "Trmv_ln_unb_var1( L, x )\n",
      "\n",
      "print( 'x after =' )\n",
      "print( x )\n",
      "\n",
      "print( np.tril( L ) * xold )\n",
      "\n",
      "print( 'x - ( np.tril( L ) * xold ) = ' ) #np.tril makes a matrix lower triangular\n",
      "print( x - ( np.tril( L ) * xold ) )"
     ],
     "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 that uses axpys"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The routine <br> <code> Trmv_ln_unb_var2( L, x ) </code>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This routine, given lower triangular $ L \\in \\mathbb{R}^{n \\times n} $ and $ x \\in \\mathbb{R}^n $ computes $ x := L x $.  The \"_ln_\" in the name of the routine indicates this is the \"lower triangular, no transpose\" matrix-vector multiplication.  \n",
      "\n",
      "The specific laff functions we will use are \n",
      "<ul>\n",
      "<li> <code> laff.axpy( alpha, x, y ) </code> which computes $ y := \\alpha x +  y  $.  </li>\n",
      "<li> <code> laff.scal( alpha, x ) </code> which computes $ x := \\alpha x $.  </li> (You will want to use this one to update $ \\chi_1 $ if you want to use PictureFrame.)\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 Trmv_ln_unb_var2(L, x):\n",
      "\n",
      "    LTL, LTR, \\\n",
      "    LBL, LBR  = flame.part_2x2(L, \\\n",
      "                               0, 0, 'BR')\n",
      "\n",
      "    xT, \\\n",
      "    xB  = flame.part_2x1(x, \\\n",
      "                         0, 'BOTTOM')\n",
      "\n",
      "    while LBR.shape[0] < L.shape[0]:\n",
      "\n",
      "        L00,  l01,      L02,  \\\n",
      "        l10t, lambda11, l12t, \\\n",
      "        L20,  l21,      L22   = flame.repart_2x2_to_3x3(LTL, LTR, \\\n",
      "                                                        LBL, LBR, \\\n",
      "                                                        1, 1, 'TL')\n",
      "\n",
      "        x0,   \\\n",
      "        chi1, \\\n",
      "        x2    = flame.repart_2x1_to_3x1(xT, \\\n",
      "                                        xB, \\\n",
      "                                        1, 'TOP')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        laff.axpy( chi1, l21, x2 )\n",
      "        laff.scal( lambda11, chi1 )\n",
      "        \n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        LTL, LTR, \\\n",
      "        LBL, LBR  = flame.cont_with_3x3_to_2x2(L00,  l01,      L02,  \\\n",
      "                                               l10t, lambda11, l12t, \\\n",
      "                                               L20,  l21,      L22,  \\\n",
      "                                               'BR')\n",
      "\n",
      "        xT, \\\n",
      "        xB  = flame.cont_with_3x1_to_2x1(x0,   \\\n",
      "                                         chi1, \\\n",
      "                                         x2,   \\\n",
      "                                         'BOTTOM')\n",
      "\n",
      "    flame.merge_2x1(xT, \\\n",
      "                    xB, x)\n",
      "\n"
     ],
     "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",
      "L = matrix( random.rand( 4,4 ) )\n",
      "x = matrix( random.rand( 4,1 ) )\n",
      "xold = matrix( random.rand( 4,1 ) )\n",
      "\n",
      "# L is not lower triangular.  We will only use the lower triangular part.\n",
      "\n",
      "print( 'L before =' )\n",
      "print( L )\n",
      "\n",
      "print( 'x before =' )\n",
      "print( x )\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "laff.copy( x, xold )   # save the original vector y\n",
      "\n",
      "Trmv_ln_unb_var2( L, x )\n",
      "\n",
      "print( 'x after =' )\n",
      "print( x )\n",
      "\n",
      "print( 'x - ( np.tril( L ) * xold ) = ' ) #np.tril makes the matrix lower triangular\n",
      "print( x - ( np.tril( L ) * xold ) )"
     ],
     "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": {}
  }
 ]
}