{
"cells": [
{
"cell_type": "markdown",
"id": "4475f0bc",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2021-10-15 Sparse Linear Algebra"
]
},
{
"cell_type": "markdown",
"id": "fa664424",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"## Last time\n",
"\n",
"* Note on p-Laplacian activity\n",
"* Weak forms and jump conditions\n",
"* 2D wave equations via gas dynamics\n",
"* Appropriate time integrators\n",
"* Hamiltonian structure and energy drift"
]
},
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"## Today\n",
"\n",
"* Recap and questions\n",
"* Sparse direct solvers\n",
" * matrix orderings\n",
" * impact on formulation\n",
" * cost scaling\n",
"* Why iterative solvers\n",
"* Discussion with Ioana Fleming"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "bf84b17a",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"laplacian_matrix (generic function with 1 method)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"using LinearAlgebra\n",
"using SparseArrays\n",
"\n",
"default(linewidth=4)\n",
"\n",
"function plot_stability(Rz, title; xlims=(-2, 2), ylims=(-2, 2))\n",
" x = LinRange(xlims[1], xlims[2], 100)\n",
" y = LinRange(ylims[1], ylims[2], 100)\n",
" heatmap(x, y, (x, y) -> abs(Rz(x + 1im*y)), c=:bwr, clims=(0, 2), aspect_ratio=:equal, title=title)\n",
"end\n",
"\n",
"struct RKTable\n",
" A::Matrix\n",
" b::Vector\n",
" c::Vector\n",
" function RKTable(A, b)\n",
" s = length(b)\n",
" A = reshape(A, s, s)\n",
" c = vec(sum(A, dims=2))\n",
" new(A, b, c)\n",
" end\n",
"end\n",
"\n",
"function rk_stability(z, rk)\n",
" s = length(rk.b)\n",
" 1 + z * rk.b' * ((I - z*rk.A) \\ ones(s))\n",
"end\n",
"\n",
"rk4 = RKTable([0 0 0 0; .5 0 0 0; 0 .5 0 0; 0 0 1 0], [1, 2, 2, 1] / 6)\n",
"\n",
"function ode_rk_explicit(f, u0; tfinal=1., h=0.1, table=rk4)\n",
" u = copy(u0)\n",
" t = 0.\n",
" n, s = length(u), length(table.c)\n",
" fY = zeros(n, s)\n",
" thist = [t]\n",
" uhist = [u0]\n",
" while t < tfinal\n",
" tnext = min(t+h, tfinal)\n",
" h = tnext - t\n",
" for i in 1:s\n",
" ti = t + h * table.c[i]\n",
" Yi = u + h * sum(fY[:,1:i-1] * table.A[i,1:i-1], dims=2)\n",
" fY[:,i] = f(ti, Yi)\n",
" end\n",
" u += h * fY * table.b\n",
" t = tnext\n",
" push!(thist, t)\n",
" push!(uhist, u)\n",
" end\n",
" thist, hcat(uhist...)\n",
"end\n",
"\n",
"function laplacian_matrix(n)\n",
" \"Laplacian in n×n periodic domain\"\n",
" h = 2 / n\n",
" rows = Vector{Int64}()\n",
" cols = Vector{Int64}()\n",
" vals = Vector{Float64}()\n",
" wrap(i) = (i + n - 1) % n + 1\n",
" idx(i, j) = (wrap(i)-1)*n + wrap(j)\n",
" stencil_diffuse = [-1, -1, 4, -1, -1] / h^2\n",
" for i in 1:n\n",
" for j in 1:n\n",
" append!(rows, repeat([idx(i,j)], 5))\n",
" append!(cols, [idx(i-1,j), idx(i,j-1), idx(i,j), idx(i+1,j), idx(i,j+1)])\n",
" append!(vals, stencil_diffuse)\n",
" end\n",
" end\n",
" sparse(rows, cols, vals)\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "73ec6381",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Two forms of acoustic wave equation"
]
},
{
"cell_type": "markdown",
"id": "12de4385",
"metadata": {
"cell_style": "split"
},
"source": [
"Divide the momentum equation through by background density and dropping the tildes yields the standard form.\n",
"$$\\begin{pmatrix} \\rho \\\\ \\mathbf u \\end{pmatrix}_t + \\nabla\\cdot \\begin{bmatrix}\n",
"\\bar\\rho \\mathbf u \\\\ \\rho \\frac{a^2}{\\bar\\rho} I \\end{bmatrix} = 0 .$$\n",
"\n",
"Examine second equation\n",
"$$ \\frac{a^2}{\\bar\\rho} \\nabla\\cdot\\big[ \\rho I \\big] = \\frac{a^2}{\\bar\\rho} \\nabla \\rho $$\n",
"and thus\n",
"$$\\begin{pmatrix} \\rho \\\\ \\mathbf u \\end{pmatrix}_t +\n",
"\\begin{bmatrix} & \\bar\\rho \\nabla\\cdot \\\\\n",
"\\frac{a^2}{\\bar\\rho} \\nabla & \\\\\n",
"\\end{bmatrix}\n",
"\\begin{pmatrix} \\rho \\\\ \\mathbf u \\end{pmatrix} = 0.$$\n"
]
},
{
"cell_type": "markdown",
"id": "30cb6485",
"metadata": {
"cell_style": "split"
},
"source": [
"Let's differentiate the first equation,\n",
"$$ \\rho_{tt} + \\bar\\rho\\nabla\\cdot(\\mathbf u_t) = 0$$\n",
"and substitute in the second equation\n",
"$$ \\rho_{tt} = a^2 \\nabla\\cdot(\\nabla \\rho)$$\n",
"\n",
"* Note: we had to assume these derivatives exist!\n",
"\n",
"We can reduce this to a first order system as\n",
"\n",
"$$\\begin{pmatrix} \\rho \\\\ \\dot \\rho \\end{pmatrix}_t + \\begin{bmatrix} & -I \\\\ -a^2 \\nabla\\cdot\\nabla & \\end{bmatrix}\n",
"\\begin{pmatrix} \\rho \\\\ \\dot\\rho \\end{pmatrix} = 0$$"
]
},
{
"cell_type": "markdown",
"id": "732c2f8b",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Question\n",
"* How is the problem size different?\n",
"* What might we be concerned about in choosing the second formulation?"
]
},
{
"cell_type": "markdown",
"id": "531b0118",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Wave operator\n",
"\n",
"$$\\begin{pmatrix} \\rho \\\\ \\dot \\rho \\end{pmatrix}_t = \\begin{bmatrix} & I \\\\ a^2 \\nabla\\cdot\\nabla & \\end{bmatrix}\n",
"\\begin{pmatrix} \\rho \\\\ \\dot\\rho \\end{pmatrix}$$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "08a74470",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"8×8 SparseMatrixCSC{Float64, Int64} with 16 stored entries:\n",
" ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ \n",
" ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0\n",
" -4.0 2.0 2.0 ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" 2.0 -4.0 ⋅ 2.0 ⋅ ⋅ ⋅ ⋅ \n",
" 2.0 ⋅ -4.0 2.0 ⋅ ⋅ ⋅ ⋅ \n",
" ⋅ 2.0 2.0 -4.0 ⋅ ⋅ ⋅ ⋅ "
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function wave_matrix(n; a=1)\n",
" Z = spzeros(n^2, n^2)\n",
" L = laplacian_matrix(n)\n",
" [Z I; -a^2*L Z]\n",
"end\n",
"wave_matrix(2)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "dd821792",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = wave_matrix(8; a=2) * .1\n",
"ev = eigvals(Matrix(A))\n",
"plot_stability(z -> rk_stability(z, rk4), \"RK4\", xlims=(-4, 4), ylims=(-4, 4))\n",
"scatter!(real(ev), imag(ev), color=:black)"
]
},
{
"cell_type": "markdown",
"id": "1aa3d38b",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Question: would forward Euler work?"
]
},
{
"cell_type": "markdown",
"id": "6395bc23",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Example 2D wave solver with RK4"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "3e5585db",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"(800, 51)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 20\n",
"A = wave_matrix(n)\n",
"x = LinRange(-1, 1, n+1)[1:end-1]\n",
"y = x\n",
"rho0 = vec(exp.(-9*((x .+ 1e-4).^2 .+ y'.^2)))\n",
"sol0 = vcat(rho0, zero(rho0))\n",
"thist, solhist = ode_rk_explicit((t, sol) -> A * sol, sol0, h=.02)\n",
"size(solhist)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "6cb81bfe",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to \n",
"│ fn = /home/jed/cu/numpde/slides/tmp.gif\n",
"└ @ Plots /home/jed/.julia/packages/Plots/yfzIY/src/animation.jl:114\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"Plots.AnimatedGif(\"/home/jed/cu/numpde/slides/tmp.gif\")"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@gif for tstep in 1:length(thist)\n",
" rho = solhist[1:n^2, tstep]\n",
" contour(x, y, reshape(rho, n, n), title=\"\\$ t = $(thist[tstep])\\$\")\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "581ea7dc",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Accuracy, conservation of mass with RK4"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "5f19eed8",
"metadata": {
"cell_style": "split",
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.02015111748435016"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"thist, solhist = ode_rk_explicit((t, sol) -> A * sol, sol0, h=.05,\n",
" tfinal=1)\n",
"\n",
"tfinal = thist[end]\n",
"M = exp(Matrix(A*tfinal))\n",
"sol_exact = M * sol0\n",
"sol_final = solhist[:, end]\n",
"norm(sol_final - sol_exact)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "276be886",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mass = vec(sum(solhist[1:n^2, :], dims=1))\n",
"plot(thist[2:end], mass[2:end] .- mass[1])"
]
},
{
"cell_type": "markdown",
"id": "e5e98c0c",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Conservation of energy with RK4"
]
},
{
"cell_type": "markdown",
"id": "40f3b161",
"metadata": {
"cell_style": "split"
},
"source": [
"## Hamiltonian structure\n",
"\n",
"We can express the total energy for our system as a sum of kinetic and potential energy.\n",
"\n",
"$$H(\\rho, \\dot\\rho) = \\underbrace{\\frac 1 2 \\int_\\Omega (\\dot\\rho)^2}_{\\text{kinetic}} + \\underbrace{\\frac{a^2}{2} \\int_\\Omega \\lVert \\nabla \\rho \\rVert^2}_{\\text{potential}}$$\n",
"\n",
"where we identify $\\rho$ as a generalized position and $\\dot\\rho$ as generalized momentum. Hamilton's equations state that the equations of motion are\n",
"\n",
"$$ \\begin{pmatrix} \\rho \\\\ \\dot\\rho \\end{pmatrix}_t = \\begin{bmatrix}\n",
"\\frac{\\partial H}{\\partial \\dot\\rho} \\\\\n",
"-\\frac{\\partial H}{\\partial \\rho}\n",
"\\end{bmatrix} = \\begin{bmatrix}\n",
"\\dot\\rho \\\\\n",
"- a^2 L \\rho\n",
"\\end{bmatrix} $$\n",
"where we have used the weak form to associate $\\int \\nabla v \\cdot \\nabla u = v^T L u$."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "eff4eda9",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function energy(sol, n)\n",
" L = laplacian_matrix(n)\n",
" rho = sol[1:end÷2]\n",
" rhodot = sol[end÷2+1:end]\n",
" kinetic = .5 * norm(rhodot)^2\n",
" potential = .5 * rho' * L * rho\n",
" kinetic + potential\n",
"end\n",
"ehist = [energy(solhist[:,i], n) for i in 1:length(thist)]\n",
"plot(thist, ehist)"
]
},
{
"cell_type": "markdown",
"id": "10eb76d7",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Velocity Verlet integrator"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "48d5c752",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"wave_verlet (generic function with 1 method)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function wave_verlet(n, u0; tfinal=1., h=0.1)\n",
" L = laplacian_matrix(n)\n",
" u = copy(u0)\n",
" t = 0.\n",
" thist = [t]\n",
" uhist = [u0]\n",
" irho = 1:n^2\n",
" irhodot = n^2+1:2*n^2\n",
" accel = -L * u[irho]\n",
" while t < tfinal\n",
" tnext = min(t+h, tfinal)\n",
" h = tnext - t\n",
" u[irho] += h * u[irhodot] + h^2/2 * accel\n",
" accel_next = -L * u[irho]\n",
" u[irhodot] += h/2 * (accel + accel_next)\n",
" accel = accel_next\n",
" t = tnext\n",
" push!(thist, t)\n",
" push!(uhist, copy(u))\n",
" end\n",
" thist, hcat(uhist...)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "e7d63d29",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to \n",
"│ fn = /home/jed/cu/numpde/slides/tmp.gif\n",
"└ @ Plots /home/jed/.julia/packages/Plots/yfzIY/src/animation.jl:114\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"Plots.AnimatedGif(\"/home/jed/cu/numpde/slides/tmp.gif\")"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"thist, solhist = wave_verlet(n, sol0, h=.04)\n",
"@gif for tstep in 1:length(thist)\n",
" rho = solhist[1:n^2, tstep]\n",
" contour(x, y, reshape(rho, n, n), title=\"\\$ t = $(thist[tstep])\\$\")\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "7471ae72",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Accuracy and conservation for Verlet"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "8a8eb24f",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm(sol_final - sol_exact) = 6.86250099252766\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"thist, solhist = wave_verlet(n, sol0, h=.05, tfinal=50)\n",
"tfinal = thist[end]\n",
"M = exp(Matrix(A*tfinal))\n",
"sol_exact = M * sol0\n",
"sol_final = solhist[:, end]\n",
"@show norm(sol_final - sol_exact)\n",
"\n",
"mass = vec(sum(solhist[1:n^2, :], dims=1))\n",
"plot(thist[2:end], mass[2:end] .- mass[1])"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "7c8d2548",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ehist = [energy(solhist[:,i], n) for i in 1:length(thist)]\n",
"plot(thist, ehist)"
]
},
{
"cell_type": "markdown",
"id": "cb268a32",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Notes on time integrators\n",
"\n",
"* We need stability on the imaginary axis for our discretization (and the physical system)\n",
"* If the model is dissipative (e.g., we didn't make the zero-viscosity assumption), then we need stability in the left half plane.\n",
"* The split form $\\rho, \\rho\\mathbf u$ form is usually used with (nonlinear) upwinding, and thus will have dissipation.\n",
"\n",
"## Runge-Kutta methods\n",
"* Easy to use, stability region designed for spatial discretization\n",
"* Energy drift generally present\n",
"\n",
"## Verlet/leapfrog/Newmark and symplectic integrators\n",
"\n",
"* These preserve the \"geometry of the Hamiltonian\"\n",
" * energy is not exactly conserved, but it doesn't drift over time\n",
" * such methods are called \"symplectic integrators\"\n",
"* May not have stability away from the imaginary axis (for dissipation)\n",
"* Most require a generalized position/momentum split, \"canonical variables\""
]
},
{
"cell_type": "markdown",
"id": "b19fa945",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Start with advection-diffusion operator in 2D\n",
"\n",
"* Eliminate Dirichlet boundary conditions around all sides"
]
},
{
"cell_type": "code",
"execution_count": 110,
"id": "7d8f9b81",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"text/plain": [
"advdiff_matrix (generic function with 1 method)"
]
},
"execution_count": 110,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function advdiff_matrix(n; kappa=1, wind=[0, 0])\n",
" h = 2 / (n + 1)\n",
" rows = Vector{Int64}()\n",
" cols = Vector{Int64}()\n",
" vals = Vector{Float64}()\n",
" idx((i, j),) = (i-1)*n + j\n",
" in_domain((i, j),) = 1 <= i <= n && 1 <= j <= n\n",
" stencil_advect = [-wind[1], -wind[2], 0, wind[1], wind[2]] / h\n",
" stencil_diffuse = [-1, -1, 4, -1, -1] * kappa / h^2\n",
" stencil = stencil_advect + stencil_diffuse\n",
" for i in 1:n\n",
" for j in 1:n\n",
" neighbors = [(i-1, j), (i, j-1), (i, j), (i+1, j), (i, j+1)]\n",
" mask = in_domain.(neighbors)\n",
" append!(rows, idx.(repeat([(i,j)], 5))[mask])\n",
" append!(cols, idx.(neighbors)[mask])\n",
" append!(vals, stencil[mask])\n",
" end\n",
" end\n",
" sparse(rows, cols, vals)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 111,
"id": "29273984",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm(A - A') = 53.665631459994955\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 111,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function my_spy(A)\n",
" cmax = norm(vec(A), Inf)\n",
" s = max(1, ceil(120 / size(A, 1)))\n",
" spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))\n",
"end\n",
"A = advdiff_matrix(5, wind=[1, 1])\n",
"@show norm(A - A')\n",
"my_spy(A)"
]
},
{
"cell_type": "markdown",
"id": "c1d1de84",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Gaussian elimination and Cholesky"
]
},
{
"cell_type": "markdown",
"id": "550a9a20",
"metadata": {
"cell_style": "split"
},
"source": [
"## $LU = A$\n",
"Given a $2\\times 2$ block matrix, the algorithm proceeds as\n",
"\\begin{split}\n",
" \\begin{bmatrix} A & B \\\\ C & D \\end{bmatrix} &=\n",
" \\begin{bmatrix} L_A & \\\\ C U_A^{-1} & 1 \\end{bmatrix}\n",
" \\begin{bmatrix} U_A & L_A^{-1} B \\\\ & S \\end{bmatrix}\n",
"\\end{split}\n",
"where $L_A U_A = A$ and\n",
"$$S = D - C \\underbrace{U_A^{-1} L_A^{-1}}_{A^{-1}} B .$$\n",
"\n",
"## Cholesky $L L^T = A$\n",
"\\begin{split}\n",
" \\begin{bmatrix} A & B^T \\\\ B & D \\end{bmatrix} &=\n",
" \\begin{bmatrix} L_A & \\\\ B L_A^{-T} & 1 \\end{bmatrix}\n",
" \\begin{bmatrix} L_A^T & L_A^{-1} B^T \\\\ & S \\end{bmatrix}\n",
"\\end{split}\n",
"where $L_A L_A^T = A$ and\n",
"$$ S = D - B \\underbrace{L_A^{-T} L_A^{-1}}_{A^{-1}} B^T .$$"
]
},
{
"cell_type": "code",
"execution_count": 112,
"id": "4e34481d",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 112,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = advdiff_matrix(5)\n",
"N = size(A, 1)\n",
"ch = cholesky(A, perm=1:N)\n",
"my_spy(sparse(ch.L))"
]
},
{
"cell_type": "markdown",
"id": "0886cb0a",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Cost of a banded solve"
]
},
{
"cell_type": "markdown",
"id": "9253467c",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"Consider an $N\\times N$ matrix with bandwidth $b$, $1 \\le b \\le N$.\n",
"\n",
"* Work one row at a time\n",
"* Each row/column of panel has $b$ nonzeros\n",
"* Schur update affects $b\\times b$ sub-matrix\n",
"* Total compute cost $N b^2$\n",
"* Storage cost $N b$\n",
"\n",
"## Question\n",
"* What bandwidth $b$ is needed for an $N = n\\times n \\times n$ cube in 3 dimensions?\n",
"* What is the memory cost?\n",
"* What is the compute cost?"
]
},
{
"cell_type": "code",
"execution_count": 96,
"id": "2417afdd",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 96,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"my_spy(sparse(ch.L))"
]
},
{
"cell_type": "markdown",
"id": "6daa8565",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Different orderings"
]
},
{
"cell_type": "code",
"execution_count": 136,
"id": "0951307c",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 136,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 20\n",
"A = advdiff_matrix(n)\n",
"heatmap(reshape(1:n^2, n, n))"
]
},
{
"cell_type": "code",
"execution_count": 137,
"id": "2eb468a1",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 137,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import Metis\n",
"perm, iperm = Metis.permutation(A)\n",
"heatmap(reshape(iperm, n, n))"
]
},
{
"cell_type": "code",
"execution_count": 138,
"id": "59338ca0",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"SuiteSparse.CHOLMOD.Factor{Float64}\n",
"type: LLt\n",
"method: simplicial\n",
"maxnnz: 8019\n",
"nnz: 8019\n",
"success: true\n"
]
},
"execution_count": 138,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cholesky(A, perm=1:n^2)"
]
},
{
"cell_type": "code",
"execution_count": 139,
"id": "96cd8cc0",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"SuiteSparse.CHOLMOD.Factor{Float64}\n",
"type: LLt\n",
"method: simplicial\n",
"maxnnz: 4031\n",
"nnz: 4031\n",
"success: true\n"
]
},
"execution_count": 139,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cholesky(A, perm=Vector{Int64}(perm))"
]
},
{
"cell_type": "markdown",
"id": "87dfedb7",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Cholesky factors in nested dissection"
]
},
{
"cell_type": "code",
"execution_count": 107,
"id": "253f096e",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"