Intro to Julia
Julia Intro for QNumerics Summer School by Raye Kimmerer.
Other resources
What scientists must know about hardware to write fast code
Linear Algebra
using LinearAlgebra
using SparseArrays
using InteractiveUtils
A = rand(100,100)
100×100 Matrix{Float64}:
0.191983 0.647775 0.455669 … 0.302315 0.879271 0.983743
0.448385 0.260415 0.839166 0.701187 0.723328 0.62188
0.86304 0.758405 0.603672 0.734188 0.859247 0.162333
0.717269 0.813982 0.846238 0.838115 0.488723 0.597041
0.959853 0.00885612 0.170822 0.732229 0.914596 0.717885
0.877949 0.698864 0.284917 … 0.483404 0.929086 0.542634
0.474096 0.41639 0.884238 0.00551019 0.475207 0.210616
0.7661 0.741337 0.850102 0.996041 0.8356 0.139016
0.0770887 0.0634909 0.512287 0.847197 0.121 0.799592
0.966052 0.430345 0.343592 0.485494 0.847403 0.478891
⋮ ⋱
0.80454 0.661706 0.199265 0.478533 0.60718 0.166643
0.531067 0.579756 0.409595 0.591683 0.781103 0.0103975
0.552929 0.160814 0.274547 0.893512 0.00906589 0.363667
0.488364 0.943916 0.497546 0.825736 0.262646 0.425224
0.753456 0.585014 0.404382 … 0.350104 0.0668837 0.16917
0.192935 0.101178 0.816083 0.185139 0.202154 0.899578
0.195405 0.236401 0.288868 0.16529 0.916725 0.364441
0.770633 0.488813 0.878098 0.344763 0.0433383 0.274334
0.845638 0.355938 0.541016 0.658037 0.174233 0.13629
Transpose and other operations are non-allocating, and instead create a wrapper of the original array. This prevents unecessary allocation and allows more opportunity to dispatch on array types.
B = transpose(A)
100×100 transpose(::Matrix{Float64}) with eltype Float64:
0.191983 0.448385 0.86304 0.717269 … 0.195405 0.770633 0.845638
0.647775 0.260415 0.758405 0.813982 0.236401 0.488813 0.355938
0.455669 0.839166 0.603672 0.846238 0.288868 0.878098 0.541016
0.642147 0.411275 0.00272084 0.592007 0.460712 0.21975 0.923579
0.488443 0.836379 0.0735283 0.275104 0.769075 0.25776 0.791245
0.0566473 0.424226 0.571391 0.198242 … 0.672614 0.338262 0.933929
0.183487 0.118888 0.867415 0.648042 0.76261 0.701895 0.667913
0.964718 0.129729 0.792813 0.63873 0.161975 0.251405 0.42952
0.721241 0.236871 0.836687 0.120308 0.915513 0.316075 0.0303537
0.677705 0.475734 0.894099 0.959848 0.252561 0.360503 0.514171
⋮ ⋱
0.394756 0.531737 0.193921 0.902524 0.168066 0.818934 0.796105
0.701097 0.509285 0.12677 0.554471 0.971499 0.458677 0.816168
0.485895 0.479988 0.870884 0.792286 0.138737 0.293608 0.976333
0.821456 0.321325 0.00167535 0.3388 0.197603 0.67803 0.0770422
0.476918 0.723504 0.796713 0.167669 … 0.693842 0.203568 0.354676
0.201813 0.93247 0.702721 0.43493 0.886152 0.223881 0.03226
0.302315 0.701187 0.734188 0.838115 0.16529 0.344763 0.658037
0.879271 0.723328 0.859247 0.488723 0.916725 0.0433383 0.174233
0.983743 0.62188 0.162333 0.597041 0.364441 0.274334 0.13629
C = copy(B)
100×100 Matrix{Float64}:
0.191983 0.448385 0.86304 0.717269 … 0.195405 0.770633 0.845638
0.647775 0.260415 0.758405 0.813982 0.236401 0.488813 0.355938
0.455669 0.839166 0.603672 0.846238 0.288868 0.878098 0.541016
0.642147 0.411275 0.00272084 0.592007 0.460712 0.21975 0.923579
0.488443 0.836379 0.0735283 0.275104 0.769075 0.25776 0.791245
0.0566473 0.424226 0.571391 0.198242 … 0.672614 0.338262 0.933929
0.183487 0.118888 0.867415 0.648042 0.76261 0.701895 0.667913
0.964718 0.129729 0.792813 0.63873 0.161975 0.251405 0.42952
0.721241 0.236871 0.836687 0.120308 0.915513 0.316075 0.0303537
0.677705 0.475734 0.894099 0.959848 0.252561 0.360503 0.514171
⋮ ⋱
0.394756 0.531737 0.193921 0.902524 0.168066 0.818934 0.796105
0.701097 0.509285 0.12677 0.554471 0.971499 0.458677 0.816168
0.485895 0.479988 0.870884 0.792286 0.138737 0.293608 0.976333
0.821456 0.321325 0.00167535 0.3388 0.197603 0.67803 0.0770422
0.476918 0.723504 0.796713 0.167669 … 0.693842 0.203568 0.354676
0.201813 0.93247 0.702721 0.43493 0.886152 0.223881 0.03226
0.302315 0.701187 0.734188 0.838115 0.16529 0.344763 0.658037
0.879271 0.723328 0.859247 0.488723 0.916725 0.0433383 0.174233
0.983743 0.62188 0.162333 0.597041 0.364441 0.274334 0.13629
When we copy the array it allocates and is no longer a wrapper
Diagonal(A)
100×100 LinearAlgebra.Diagonal{Float64, Vector{Float64}}:
0.191983 ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅
⋅ 0.260415 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 0.603672 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.592007 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋮ ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 0.16529 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 0.0433383 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.13629
subtypes(LinearAlgebra.Factorization)
21-element Vector{Any}:
LinearAlgebra.AdjointFactorization
LinearAlgebra.BunchKaufman
LinearAlgebra.Cholesky
LinearAlgebra.CholeskyPivoted
LinearAlgebra.Eigen
LinearAlgebra.GeneralizedEigen
LinearAlgebra.GeneralizedSVD
LinearAlgebra.GeneralizedSchur
LinearAlgebra.Hessenberg
LinearAlgebra.LDLt
⋮
LinearAlgebra.QR
LinearAlgebra.QRCompactWY
LinearAlgebra.QRPivoted
LinearAlgebra.SVD
LinearAlgebra.Schur
LinearAlgebra.TransposeFactorization
SparseArrays.CHOLMOD.Factor
SparseArrays.SPQR.QRSparse
SparseArrays.UMFPACK.UmfpackLU
issymmetric(A*A')
true
How do we communicate this to the dispatch system?
D = Symmetric(A*A')
100×100 LinearAlgebra.Symmetric{Float64, Matrix{Float64}}:
26.2458 21.9174 21.5573 24.6514 … 20.7541 18.731 22.7841 20.2111
21.9174 35.1197 26.1359 29.417 24.5176 24.1426 26.6386 25.327
21.5573 26.1359 34.1386 27.923 24.3698 22.3534 25.7664 24.7429
24.6514 29.417 27.923 38.6597 26.3025 23.0557 29.8224 28.325
19.9921 23.8199 23.0674 24.0203 22.0853 20.6375 24.1552 23.4745
21.4973 24.9992 25.2468 27.0158 … 23.7177 22.1222 25.3736 24.9971
20.9693 26.9475 25.0604 29.2302 24.4049 21.9628 27.4436 26.8895
24.2297 27.5378 26.9258 30.4282 26.2334 24.7361 29.3078 28.8234
21.0203 27.013 24.2184 27.217 24.9104 22.4303 25.8028 27.0007
21.7824 25.1906 25.0665 26.9468 21.7958 21.8307 23.1402 23.3865
⋮ ⋱
20.9425 25.3974 24.6736 27.771 23.0963 21.4178 26.8237 25.0653
20.9833 25.5236 26.0744 27.0185 22.3034 25.0623 27.5233 25.6094
19.3539 23.4718 22.1093 25.2943 23.1247 21.3525 23.594 23.6654
22.3317 24.7157 25.5742 27.1824 21.7149 21.3288 23.8024 24.0661
20.6212 23.8606 25.6432 27.0251 … 21.4184 20.5398 24.5434 24.8289
20.7541 24.5176 24.3698 26.3025 30.3977 22.0569 25.6438 22.6868
18.731 24.1426 22.3534 23.0557 22.0569 28.2522 21.5008 22.6817
22.7841 26.6386 25.7664 29.8224 25.6438 21.5008 34.56 26.0204
20.2111 25.327 24.7429 28.325 22.6868 22.6817 26.0204 33.7825
Now we can dispatch on the fact that this matrix is symmetric, allowing us to select more efficient algorithms for this case.
Symmetric(A)
100×100 LinearAlgebra.Symmetric{Float64, Matrix{Float64}}:
0.191983 0.647775 0.455669 … 0.302315 0.879271 0.983743
0.647775 0.260415 0.839166 0.701187 0.723328 0.62188
0.455669 0.839166 0.603672 0.734188 0.859247 0.162333
0.642147 0.411275 0.00272084 0.838115 0.488723 0.597041
0.488443 0.836379 0.0735283 0.732229 0.914596 0.717885
0.0566473 0.424226 0.571391 … 0.483404 0.929086 0.542634
0.183487 0.118888 0.867415 0.00551019 0.475207 0.210616
0.964718 0.129729 0.792813 0.996041 0.8356 0.139016
0.721241 0.236871 0.836687 0.847197 0.121 0.799592
0.677705 0.475734 0.894099 0.485494 0.847403 0.478891
⋮ ⋱
0.394756 0.531737 0.193921 0.478533 0.60718 0.166643
0.701097 0.509285 0.12677 0.591683 0.781103 0.0103975
0.485895 0.479988 0.870884 0.893512 0.00906589 0.363667
0.821456 0.321325 0.00167535 0.825736 0.262646 0.425224
0.476918 0.723504 0.796713 … 0.350104 0.0668837 0.16917
0.201813 0.93247 0.702721 0.185139 0.202154 0.899578
0.302315 0.701187 0.734188 0.16529 0.916725 0.364441
0.879271 0.723328 0.859247 0.916725 0.0433383 0.274334
0.983743 0.62188 0.162333 0.364441 0.274334 0.13629
This is 'casting' to a symmetric matrix, it will happily create a symmetric matrix out of a non-symmetric one.
Sparse Arrays
A = sprand(100,100,0.1)
100×100 SparseArrays.SparseMatrixCSC{Float64, Int64} with 989 stored entries:
⎡⠀⡄⠀⠢⢈⢠⠀⠁⠄⠠⢀⠠⠀⢄⠁⡁⠀⠐⡈⠀⡀⡁⠀⠤⡀⠐⡡⠢⠈⠄⠀⠆⠂⠀⠠⠀⠐⠠⠀⠐⎤
⎢⠀⠙⠀⠂⠠⠠⠀⢂⠔⠃⠐⠀⠂⠀⠩⢥⠀⠃⠄⠀⠠⢀⠐⠀⠂⠀⠨⠀⠨⠡⠀⠀⠠⠄⡹⠀⠀⠈⠰⠀⎥
⎢⠐⠄⠐⠔⠂⠠⠀⠅⠄⠅⠠⠀⠄⡊⠤⠀⠠⠀⠦⠐⢐⠀⠀⠀⠉⠉⢐⢀⠁⠈⢀⢀⠀⠀⠀⠄⠠⠠⠀⠀⎥
⎢⠡⠀⠠⠚⡄⢂⠀⠁⠂⠀⠠⠏⠈⡀⠀⠠⠀⠁⠀⠄⡐⠆⠸⣌⠀⠰⢤⠀⡃⠀⢠⠈⠀⠂⠠⠁⠀⠀⠐⢀⎥
⎢⠀⢀⠠⡀⠀⠀⠀⠀⠀⠁⠠⠀⠄⠅⠶⠃⠀⠀⠀⠀⢂⠈⢢⡢⠀⠀⠁⠀⢠⠤⠠⠈⢠⠤⠩⠀⠨⠠⣤⠀⎥
⎢⠐⠠⠀⡀⠀⠠⠇⠀⠐⡠⠁⠀⡀⡄⠀⠜⠐⠀⢠⠀⣀⠄⢤⠈⠀⠴⠰⠰⡠⠔⠀⠀⡀⡀⠠⠀⠈⡠⠄⠀⎥
⎢⢄⠊⢨⠄⢠⠘⠀⠭⠂⡊⡀⠀⢀⠁⡀⠀⠄⠧⠀⠆⠠⡃⠴⠀⡠⢀⠒⠊⠀⠤⢠⠠⠔⠌⠐⠀⠄⡡⠃⠄⎥
⎢⠄⠠⠀⠀⠀⠡⠀⢀⠀⠄⡀⠔⠨⢄⠀⡌⢄⠀⠐⡡⡀⠨⠠⠀⠨⠆⠰⠠⠠⢄⠐⠀⠈⠄⠀⡀⠈⠀⡠⠂⎥
⎢⢀⡸⠌⠂⠂⠦⠠⢄⠀⢕⢀⡁⠀⠀⠀⡄⠐⢆⠀⠂⠀⠄⠀⠀⠄⢑⠈⠂⠀⠀⢆⠠⠄⢀⢀⠀⠠⡀⠢⠀⎥
⎢⠀⠀⠠⠐⠀⠥⠠⠄⡁⡁⠀⠀⠀⠂⠀⡂⠀⠸⠀⠀⠢⠸⢠⠂⠀⡤⢠⠈⠠⠠⠠⠀⢀⠄⠀⢰⠈⡀⠀⠀⎥
⎢⠐⠠⠄⡠⡈⠂⢀⠀⠄⠈⠀⡠⠀⠀⠐⠂⠀⠀⠀⠘⠐⠀⠆⠈⡠⠂⠨⠀⡃⠱⠈⢠⠀⠀⠘⠁⠂⢰⠈⢄⎥
⎢⡫⠁⠀⠧⠀⠅⡉⠂⠀⠆⠀⠠⠠⠀⢀⠈⠀⠠⠐⠆⢃⢀⠨⠂⠠⠀⢐⠐⢂⠄⠐⠐⠡⢀⠀⠀⠂⢄⠀⠄⎥
⎢⢠⠄⠀⡀⠄⠘⡀⠔⠀⠒⢁⠂⠀⡇⠀⠀⡀⡀⠀⠂⠀⡀⠀⠂⠚⢀⠈⡈⠐⡐⠒⠈⠀⡂⢂⠀⠀⠰⡅⠀⎥
⎢⠀⡃⠀⠂⠄⠆⠀⠂⠂⠈⡕⠀⡂⠃⠈⠀⠀⠂⠂⣃⢐⠀⢊⠀⠂⠐⡠⠠⠒⠀⠐⢁⠂⠀⠬⠀⠐⠄⠐⠀⎥
⎢⢠⠁⠠⠒⢁⠐⠀⢄⡐⠂⠀⠒⡚⠀⠌⠧⠂⠀⠀⡀⠈⠀⠄⠐⢐⠀⡰⠀⢤⠀⠐⡐⠠⠠⠂⠠⠠⠐⠀⠄⎥
⎢⠀⠑⠊⠁⠑⠑⢀⡀⠂⣁⠀⠐⠀⠂⠀⠄⠸⠨⠀⠀⣈⡑⠠⠈⠢⠀⠰⠀⢀⡄⣐⡀⠒⠄⠤⠐⠀⠀⡑⣀⎥
⎢⠀⠂⠠⠀⠀⢂⠀⠇⠨⠔⠌⠆⠈⠄⠠⡀⠂⡁⠠⠙⢰⠂⠀⠠⠀⠊⢄⠔⠀⠠⡀⠀⠒⠃⠂⢀⠂⠂⢉⠐⎥
⎢⠁⠀⠀⡀⢀⠨⠐⡱⡀⡔⠤⠌⠀⠀⠀⠉⠐⠀⠀⠀⠈⡀⢀⠀⠐⠰⢐⠠⢐⢑⠘⠀⠀⠀⠠⠂⠂⠐⢔⢐⎥
⎢⠠⠢⠀⠣⢐⡀⠀⡀⠚⡅⠀⠆⢁⠁⢐⠀⢀⠓⠀⡐⡈⠐⡀⠀⠐⡀⠀⡂⢀⠐⠐⡀⡠⠐⠀⠀⠀⠠⠀⠆⎥
⎣⠐⢀⠐⠠⠈⠄⡀⠁⠐⠀⠔⠀⠆⢂⡀⠅⠅⠄⠀⡀⡎⠀⡸⠀⠀⠀⠠⠑⡳⠀⠀⠀⠀⠀⠖⠀⠃⠀⢂⠐⎦