{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Julia Tutorial: Diagonal vs diagm\n",
    "\n",
    "This tutorial discusses the three ways to make diagonal matrices in Julia:  \n",
    "`Diagonal` and `diagm` and `spdiag`\n",
    "\n",
    "Of these, the `Diagonal` command is by far the most efficient.  \n",
    "Always use `Diagonal` instead of `diagm` for a pure diagonal matrix!\n",
    "\n",
    "Jeff Fessler, University of Michigan  \n",
    "2019-01-15, Julia 1.0.3  \n",
    "2020-08-05, Julia 1.5.0  \n",
    "2021-08-23, Julia 1.6.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra # this package comes with Julia, no need to install it"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The `Diagonal` function  \n",
    "The `Diagonal` function in the `LinearAlgebra` package\n",
    "creates a diagonal \"matrix\" from a vector of its diagonal elements.  \n",
    "\n",
    "I put \"matrix\" in quotes because the data type of the output of `Diagonal`\n",
    "is not an ordinary 2D Array but rather a special `Diagonal` type\n",
    "that behaves like a matrix but is stored efficiently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example of creating a diagonal matrix\n",
    "# (internally Julia stores only the nonzero values)\n",
    "D = Diagonal(10:10:50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can access individual elements of `D` just like one would with an ordinary matrix.\n",
    "D[3,3], D[3,1], D[1,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If we access multiple elements of `D` then we get sparse array\n",
    "D[1:2,1:4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The `diagm` function  \n",
    "The `diagm` function in the `LinearAlgebra` package\n",
    "creates a diagonal matrix from a vector of its diagonal elements,\n",
    "or a banded matrix from multiple band vectors.\n",
    "\n",
    "Here, matrix means a standard 2D `Array` in Julia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = diagm(0 => 10:10:50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The `spdiagm` function  \n",
    "The `spdiagm` function in the `SparseArrays` package\n",
    "creates a diagonal \"matrix\" from a vector of its diagonal elements,\n",
    "or a banded matrix from multiple band vectors.\n",
    "\n",
    "Here, \"matrix\" is quoted again because it is not a standard 2D Array in Julia,\n",
    "but rather a special `SparseMatrix` type that efficiently stores\n",
    "only the nonzero elements and their array indexes.\n",
    "However, for the special case of a diagonal matrix\n",
    "the sparse data type is not as efficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using SparseArrays: spdiagm\n",
    "S = spdiagm(0 => 10:10:50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparing matrix multiplication speeds\n",
    "\n",
    "Next we compare the execution times of multiplying\n",
    "each of the above 3 diagonal matrix data types\n",
    "by a vector.\n",
    "The differences are dramatic!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "using BenchmarkTools\n",
    "N = 2000 # medium sized problem to illustrate the differences\n",
    "\n",
    "D = Diagonal(1:N)\n",
    "M = diagm(0 => 1:N)\n",
    "S = spdiagm(0 => 1:N)\n",
    "\n",
    "x = rand(N) # random test vector\n",
    "y = zeros(N) # space for output\n",
    "\n",
    "@btime y = D * x; # fastest\n",
    "@btime y = S * x; # 3× slower\n",
    "@btime y = M * x; # >100× slower!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The preceding test has to reallocate y everytime, which is inefficient.\n",
    "Here is a faster version that overwrites `y` each time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "@btime mul!(y, D, x); # fastest\n",
    "@btime mul!(y, S, x); # 3× slower\n",
    "@btime mul!(y, M, x); # >100× slower!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The notation mul!(y, A, x) does not look much like the math y=A*x\n",
    "# The following package solves this issue\n",
    "using InplaceOps\n",
    "@! y = D * x; # typical use\n",
    "@btime @! $y = D * x; # very fast\n",
    "@btime @! $y = S * x; # slower\n",
    "@btime @! $y = M * x; # very slow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### So what use is `diagm` then?\n",
    "\n",
    "For making banded matrices with multiple diagonals, as follows.\n",
    "NOT for making an ordinary diagonal matrix!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "diagm(1 => 1:5, -2 => 3:7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rectangular Diagonal matrices  \n",
    "The `Diagonal` function in the `LinearAlgebra` package\n",
    "can also create a diagonal matrix\n",
    "from the diagonal elements of a (possibly rectangular) matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct a diagonal matrix out of the elements of a (possibly rectangular) matrix\n",
    "# When the matrix is rectangular, the size of the output diagonal matrix\n",
    "# is the minimum of the dimensions of the dimensions of the input matrix\n",
    "B = [1 2; 3 4; 5 6] # 3 × 2\n",
    "display(B)\n",
    "DB = Diagonal(B) # 2 × 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = [1 2 3; 4 5 6] # 2 × 3\n",
    "display(C)\n",
    "DC = Diagonal(C) # 2 × 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One must again be careful with the difference\n",
    "between a 1D vector and a 2D array that\n",
    "happens to have one dimension equal to 1,\n",
    "as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Diagonal(ones(3)) # a 1D vector works as expected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colmat = ones(3,1) # This is a 2D array of size \"3 × 1\"\n",
    "display(colmat)\n",
    "Diagonal(colmat) # so this is 1 × 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rowmat = ones(1,3) # This is a 2D array of size 1 × 3\n",
    "display(rowmat)\n",
    "Diagonal(rowmat) # this is 1 × 1"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.2",
   "language": "julia",
   "name": "julia-1.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}