# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.3.3
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# #### Julia Tutorial: Diagonal vs diagm
#
# This tutorial discusses the three ways to make diagonal matrices in Julia:  
# `Diagonal` and `diagm` and `spdiag`
#
# Of these, the `Diagonal` command is by far the most efficient.  
# Always use `Diagonal` instead of `diagm` for a pure diagonal matrix!
#
# Jeff Fessler, University of Michigan  
# 2019-01-15, Julia 1.0.3  
# 2020-08-05, Julia 1.5.0  
# 2021-08-23, Julia 1.6.2

using LinearAlgebra # this package comes with Julia, no need to install it

# #### The `Diagonal` function  
# The `Diagonal` function in the `LinearAlgebra` package
# creates a diagonal "matrix" from a vector of its diagonal elements.  
#
# I put "matrix" in quotes because the data type of the output of `Diagonal`
# is not an ordinary 2D Array but rather a special `Diagonal` type
# that behaves like a matrix but is stored efficiently.

# Example of creating a diagonal matrix
# (internally Julia stores only the nonzero values)
D = Diagonal(10:10:50)

# We can access individual elements of `D` just like one would with an ordinary matrix.
D[3,3], D[3,1], D[1,2]

# If we access multiple elements of `D` then we get sparse array
D[1:2,1:4]

# #### The `diagm` function  
# The `diagm` function in the `LinearAlgebra` package
# creates a diagonal matrix from a vector of its diagonal elements,
# or a banded matrix from multiple band vectors.
#
# Here, matrix means a standard 2D `Array` in Julia.

M = diagm(0 => 10:10:50)

# #### The `spdiagm` function  
# The `spdiagm` function in the `SparseArrays` package
# creates a diagonal "matrix" from a vector of its diagonal elements,
# or a banded matrix from multiple band vectors.
#
# Here, "matrix" is quoted again because it is not a standard 2D Array in Julia,
# but rather a special `SparseMatrix` type that efficiently stores
# only the nonzero elements and their array indexes.
# However, for the special case of a diagonal matrix
# the sparse data type is not as efficient.

using SparseArrays: spdiagm
S = spdiagm(0 => 10:10:50)

# #### Comparing matrix multiplication speeds
#
# Next we compare the execution times of multiplying
# each of the above 3 diagonal matrix data types
# by a vector.
# The differences are dramatic!

# +
using BenchmarkTools
N = 2000 # medium sized problem to illustrate the differences

D = Diagonal(1:N)
M = diagm(0 => 1:N)
S = spdiagm(0 => 1:N)

x = rand(N) # random test vector
y = zeros(N) # space for output

@btime y = D * x; # fastest
@btime y = S * x; # 3× slower
@btime y = M * x; # >100× slower!
# -


# The preceding test has to reallocate y everytime, which is inefficient.
# Here is a faster version that overwrites `y` each time:

@btime mul!(y, D, x); # fastest
@btime mul!(y, S, x); # 3× slower
@btime mul!(y, M, x); # >100× slower!


# The notation mul!(y, A, x) does not look much like the math y=A*x
# The following package solves this issue
using InplaceOps
@! y = D * x; # typical use
@btime @! $y = D * x; # very fast
@btime @! $y = S * x; # slower
@btime @! $y = M * x; # very slow



# #### So what use is `diagm` then?
#
# For making banded matrices with multiple diagonals, as follows.
# NOT for making an ordinary diagonal matrix!

diagm(1 => 1:5, -2 => 3:7)


# ### Rectangular Diagonal matrices  
# The `Diagonal` function in the `LinearAlgebra` package
# can also create a diagonal matrix
# from the diagonal elements of a (possibly rectangular) matrix.

# Construct a diagonal matrix out of the elements of a (possibly rectangular) matrix
# When the matrix is rectangular, the size of the output diagonal matrix
# is the minimum of the dimensions of the dimensions of the input matrix
B = [1 2; 3 4; 5 6] # 3 × 2
display(B)
DB = Diagonal(B) # 2 × 2

C = [1 2 3; 4 5 6] # 2 × 3
display(C)
DC = Diagonal(C) # 2 × 2

# One must again be careful with the difference
# between a 1D vector and a 2D array that
# happens to have one dimension equal to 1,
# as follows.

Diagonal(ones(3)) # a 1D vector works as expected

colmat = ones(3,1) # This is a 2D array of size "3 × 1"
display(colmat)
Diagonal(colmat) # so this is 1 × 1

rowmat = ones(1,3) # This is a 2D array of size 1 × 3
display(rowmat)
Diagonal(rowmat) # this is 1 × 1