" ] }, { "cell_type": "markdown", "id": "05dcb88e", "metadata": {}, "source": [ "Almost all licenses require these conditions (from BSD-2-Clause)\n", "1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.\n", "2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution." ] }, { "cell_type": "markdown", "id": "66ef5772", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# License compatibility" ] }, { "cell_type": "markdown", "id": "fdecb8ef", "metadata": { "cell_style": "split" }, "source": [ "@github copilot, with "public code" blocked, emits large chunks of my copyrighted code, with no attribution, no LGPL license. For example, the simple prompt "sparse matrix transpose, cs_" produces my cs_transpose in CSparse. My code on left, github on right. Not OK. pic.twitter.com/sqpOThi8nf
— Tim Davis (@DocSparse) October 16, 2022
![](\"https://dwheeler.com/essays/floss-license-slide-image.png\")
\n",
"\n",
"[CC BY-SA-3.0](https://creativecommons.org/licenses/by-sa/3.0/) [David A. Wheeler](https://dwheeler.com/essays/floss-license-slide.html)\n",
"\n",
"## Recommendation\n",
"\n",
"Use the most standard license that meets your needs.\n",
"\n",
"License proliferation saps time and discourages users. See [OSI recommendations](https://opensource.org/licenses/category) and [Todd Gamblin's recommendations for LLNL](https://www.osti.gov/servlets/purl/1782513/)\n",
"\n",
"* [SPDX License names](https://spdx.org/licenses/)"
]
},
{
"cell_type": "markdown",
"id": "e28d6e47",
"metadata": {
"cell_style": "split"
},
"source": [
"* Permissive (MIT, BSD, Apache-2.0)\n",
" * Can put into proprietary software, as a library or with customization.\n",
" * Can't plagiarize (attribution required)\n",
"* Weak Copyleft (LGPL, MPL)\n",
" * Can use as a *dynamically linked* library in proprietary products\n",
" * Must be able to link modifications of library.\n",
" * Modifications are subject to LGPL (anyone who has the product can have the code, modify, redistribute)\n",
"* Strong Copyleft (GPL, AGPL)\n",
" * \"viral\" -- anyone who has/uses the software product can access the source code, modify, and redistribute."
]
},
{
"cell_type": "markdown",
"id": "4206c8c2",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Software patents and inbound=outbound"
]
},
{
"cell_type": "markdown",
"id": "4f11a8e5",
"metadata": {
"cell_style": "split"
},
"source": [
"* Patent litigation is hugely destructive for open source communities.\n",
"* There is no legal consensus about whether making a contribution implies a grant to any patents that may be infringed by that contribution.\n",
"* Apache-2.0 $\\approx$ MIT/BSD plus patent grant and indemnification\n",
"\n",
"## Problem\n",
"GPL/LGPL forbid adding new restrictions. The Free Software Foundation believes the Apache-2.0 indemnification clause violates GPL/LGPL.\n",
"## Solutions\n",
"* `(Apache-2.0 OR MIT)` (Rust Language)\n",
"* `Apache-2.0 WITH LLVM-exception`"
]
},
{
"cell_type": "markdown",
"id": "a0f3b8eb",
"metadata": {
"cell_style": "split"
},
"source": [
"# Contributor License Agreements\n",
"\n",
"* [Apache-2.0 CLA](https://www.apache.org/licenses/contributor-agreements.html): paper signed documents granting terms of Apache-2.0 license.\n",
" * One per person or organization (signed by employer/counsel)\n",
" * Creates friction to contribute, bookkeeping for maintainers, coarse granularity.\n",
"* [Developers Certificate of Origin](https://developercertificate.org/)\n",
" * Developed for the Linux Kernel in wake of SCO lawsuits.\n",
" * Low-friction for contributors, fine granularity with low bookkeeping.\\\n",
" * Confirms authority to contribute and **inbound=outbound**"
]
},
{
"cell_type": "markdown",
"id": "cb8e7b6b",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# PETSc experiments"
]
},
{
"cell_type": "markdown",
"id": "4d8b65e2",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"* Compare domain decomposition and multigrid preconditioning\n",
" * `-pc_type asm` (Additive Schwarz)\n",
" * `-pc_asm_type basic` (symmetric, versus `restrict`)\n",
" * `-pc_asm_overlap 2` (increase overlap)\n",
" * Effect of direct subdomain solver: `-sub_pc_type lu`"
]
},
{
"cell_type": "markdown",
"id": "ba46a7e2",
"metadata": {
"cell_style": "split"
},
"source": [
"* Symmetric example: `src/snes/examples/tutorials/ex5.c`\n",
"* Nonsymmetric example: `src/snes/examples/tutorials/ex19.c`\n",
"* Compare preconditioned versus unpreconditioned norms.\n",
"* Compare BiCG versus GMRES\n",
" * `-pc_type mg` (Geometric Multigrid)\n",
"* Use monitors:\n",
" * `-ksp_monitor_true_residual`\n",
" * `-ksp_monitor_singular_value`\n",
" * `-ksp_converged_reason`\n",
"* Explain methods: `-snes_view`\n",
"* Performance info: `-log_view`"
]
},
{
"cell_type": "markdown",
"id": "e1d6705f",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Geometric multigrid and the spectrum\n",
"\n",
"Multigrid uses a hierarchy of scales to produce an operation \"$M^{-1}$\" that can be applied in $O(n)$ (\"linear\") time and such that $\\lVert I - M^{-1} A \\rVert \\le \\rho < 1$ independent of the problem size.\n",
"We'll consider the one dimensional Laplacian using the stencil\n",
"$$ \\Big(\\frac{h^2}{2}\\Big) \\frac{1}{h^2} \\begin{bmatrix} -1 & 2 & -1 \\end{bmatrix} $$"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "f65e07be",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"laplace1d (generic function with 1 method)"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function laplace1d(n)\n",
" \"1D Laplacion with Dirichlet boundary conditions eliminated\"\n",
" h = 2 / (n + 1)\n",
" x = LinRange(-1, 1, n+2)[2:end-1]\n",
" A = diagm(0 => ones(n),\n",
" -1 => -.5*ones(n-1),\n",
" +1 => -.5*ones(n-1))\n",
" x, A # Hermitian(A)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "592a5d02",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"L[1:2] = [0.00048371770801235576, 0.0019344028664056293]\n",
"L[end - 1:end] = [1.9980655971335943, 1.999516282291988]\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x, A = laplace1d(100)\n",
"#@show A[1:4,1:4]\n",
"L, X = eigen(A)\n",
"@show L[1:2]\n",
"@show L[end-1:end]\n",
"scatter(L, title=\"cond = $(L[end]/L[1])\", legend=:none)"
]
},
{
"cell_type": "markdown",
"id": "44c367f9",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# What do the eigenvectors look like?"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "609fd56e",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x, A = laplace1d(100)\n",
"L, X = eigen(A)\n",
"plot(x, X[:,1:3], label=L[1:3]', legend=:bottomleft)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"id": "5f6ec119",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(x, X[:, end-2:end], label=L[end-2:end]', legend=:bottomleft)"
]
},
{
"cell_type": "markdown",
"id": "28ca432c",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Condition number grows with grid refinement"
]
},
{
"cell_type": "code",
"execution_count": 41,
"id": "2fb37509",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ns = 2 .^ (2:8)\n",
"eigs = vcat([eigvals(laplace1d(n)[2])[[1, end]] for n in ns]'...)\n",
"\n",
"scatter(ns, eigs, label=[\"min eig(A)\" \"max eig(A)\"])\n",
"plot!(n -> 1/n^2, label=\"\\$1/n^2\\$\", xscale=:log10, yscale=:log10)"
]
},
{
"cell_type": "markdown",
"id": "f0471c50",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Fourier analysis perspective\n",
"\n",
"Consider the basis $\\phi(x, \\theta) = e^{i \\theta x}$. If we choose the grid $x \\in h \\mathbb Z$ with grid size $h$ then we can resolve frequencies $\\lvert \\theta \\rvert \\le \\pi/h$."
]
},
{
"cell_type": "code",
"execution_count": 42,
"id": "50fe3cd4",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"plot_symbol (generic function with 2 methods)"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function symbol(S, theta)\n",
" if length(S) % 2 != 1\n",
" error(\"Length of stencil must be odd\")\n",
" end\n",
" w = length(S) ÷ 2\n",
" phi = exp.(1im * (-w:w) * theta')\n",
" vec(S * phi) # not! (S * phi)'\n",
"end\n",
"\n",
"function plot_symbol(S, deriv=2; plot_ref=true, n_theta=30)\n",
" theta = LinRange(-pi, pi, n_theta)\n",
" sym = symbol(S, theta)\n",
" rsym = real.((-1im)^deriv * sym)\n",
" fig = plot(theta, rsym, marker=:circle, label=\"stencil\")\n",
" if plot_ref\n",
" plot!(fig, th -> th^deriv, label=\"exact\")\n",
" end\n",
" fig\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 43,
"id": "765ddfb1",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_symbol([1 -2 1])\n",
"#plot!(xlims=(-1, 1))"
]
},
{
"cell_type": "markdown",
"id": "58271476",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Analytically computing smallest eigenvalue\n",
"\n",
"The longest wavelength for a domain size of 2 with Dirichlet boundary conditions is 4. The frequency is $\\theta = 2\\pi/\\lambda = \\pi/2$. The `symbol` function works on an integer grid. We can transform via $\\theta \\mapsto \\theta h$."
]
},
{
"cell_type": "code",
"execution_count": 44,
"id": "9bc3f0fe",
"metadata": {
"cell_style": "center"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scatter(ns, eigs, label=[\"min eig(A)\" \"max eig(A)\"])\n",
"plot!(n -> 1/n^2, label=\"\\$1/n^2\\$\", xscale=:log10, yscale=:log10)\n",
"theta_min = pi ./ (ns .+ 1)\n",
"symbol_min = -real(symbol([1 -2 1], theta_min))\n",
"scatter!(ns, symbol_min / 2, shape=:utriangle)"
]
},
{
"cell_type": "markdown",
"id": "9939381d",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Damping properties of Richardson/Jacobi relaxation\n",
"\n",
"Recall that we would like for $I - w A$ to have a small norm, such that powers (repeat iterations) will cause the error to decay rapidly."
]
},
{
"cell_type": "code",
"execution_count": 47,
"id": "10a1eedc",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x, A = laplace1d(80)\n",
"ws = LinRange(.3, 1.001, 50)\n",
"radius = [opnorm(I - w*A) for w in ws]\n",
"plot(ws, radius, xlabel=\"w\")"
]
},
{
"cell_type": "markdown",
"id": "087ee3fa",
"metadata": {
"cell_style": "split"
},
"source": [
"* The spectrum of $A$ runs from $\\theta_{\\min}^2$ up to 2. If $w > 1$, then $\\lVert I - w A \\rVert > 1$ because the operation amplifies the high frequencies (associated with the eigenvalue of 2).\n",
"* The value of $w$ that minimizes the norm is slightly less than 1, but the convergence rate is very slow (only barely less than 1)."
]
},
{
"cell_type": "markdown",
"id": "24691cd5",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Symbol of damping"
]
},
{
"cell_type": "code",
"execution_count": 51,
"id": "6cd7ddf5",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"w = .7\n",
"plot_symbol([0 1 0] - w * [-.5 1 -.5], 0; plot_ref=false)"
]
},
{
"cell_type": "markdown",
"id": "8b747093",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"* Evidently it is very difficult to damp low frequencies.\n",
"* This makes sense because $A$ and $I - wA$ move information only one grid point per iteration.\n",
"* It also makes sense because the polynomial needs to be 1 at the origin, and the low frequencies have eigenvalues very close to zero."
]
},
{
"cell_type": "markdown",
"id": "a3655413",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Coarse grids: Make low frequencies high (again)\n",
"\n",
"As in domain decomposition, we will express our \"coarse\" subspace, consisting of a grid $x \\in 2h\\mathbb Z$, in terms of its interpolation to the fine space. Here, we'll use linear interpolation."
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "3de3af3e",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function interpolate(m, stride=2)\n",
" s1 = (stride - 1) / stride\n",
" s2 = (stride - 2) / stride\n",
" P = diagm(0 => ones(m),\n",
" -1 => s1*ones(m-1), +1 => s1*ones(m-1),\n",
" -2 => s2*ones(m-2), +2 => s2*ones(m-2))\n",
" P[:, 1:stride:end]\n",
"end\n",
"n = 50; x, A = laplace1d(n)\n",
"P = interpolate(n, 2)\n",
"plot(x, P[:, 4:6], marker=:auto, xlims=(-1, 0))"
]
},
{
"cell_type": "code",
"execution_count": 61,
"id": "af3a414c",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L, X = eigen(A)\n",
"u_h = X[:, 5]\n",
"u_2h = .5 * P' * u_h\n",
"plot(x, [u_h, P * u_2h])"
]
},
{
"cell_type": "markdown",
"id": "c6606c70",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Galerkin approximation of $A$ in coarse space\n",
"\n",
"$$ A_{2h} u_{2h} = P^T A_h P u_{2h} $$"
]
},
{
"cell_type": "code",
"execution_count": 58,
"id": "789af8ff",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"size(A) = (50, 50)\n",
"size(A_2h) = (25, 25)\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x, A = laplace1d(n)\n",
"P = interpolate(n)\n",
"@show size(A)\n",
"A_2h = P' * A * P\n",
"@show size(A_2h)\n",
"L_2h = eigvals(A_2h)\n",
"plot(L_2h, title=\"cond $(L_2h[end]/L_2h[1])\")"
]
},
{
"cell_type": "code",
"execution_count": 59,
"id": "42bf3cf5",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"my_spy(A_2h)"
]
},
{
"cell_type": "markdown",
"id": "b2a9adf1",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Coarse grid correction\n",
"\n",
"Consider the $A$-orthogonal projection onto the range of $P$,\n",
"$$ S_c = P A_{2h}^{-1} P^T A $$"
]
},
{
"cell_type": "code",
"execution_count": 60,
"id": "629d978d",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Sc = P * (A_2h \\ P' * A)\n",
"Ls, Xs = eigen(I - Sc)\n",
"scatter(real.(Ls))"
]
},
{
"cell_type": "markdown",
"id": "dc3f5c9f",
"metadata": {
"cell_style": "split"
},
"source": [
"* This spectrum is typical for a projector. If $u$ is in the range of $P$, then $S_c u = u$. Why?\n",
"* For all vectors $v$ that are $A$-orthogonal to the range of $P$, we know that $S_c v = 0$. Why?"
]
},
{
"cell_type": "markdown",
"id": "38497d54",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# A two-grid method"
]
},
{
"cell_type": "code",
"execution_count": 66,
"id": "e4533ed8",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Warning: No strict ticks found\n",
"â”” @ PlotUtils /home/jed/.julia/packages/PlotUtils/NE7U1/src/ticks.jl:191\n",
"┌ Warning: No strict ticks found\n",
"â”” @ PlotUtils /home/jed/.julia/packages/PlotUtils/NE7U1/src/ticks.jl:191\n"
]
}
],
"source": [
"w = .67\n",
"T = (I - w*A) * (I - Sc) * (I - w*A)\n",
"Lt = eigvals(T)\n",
"scatter(Lt)"
]
},
{
"cell_type": "markdown",
"id": "089b2095",
"metadata": {
"cell_style": "split"
},
"source": [
"* Can analyze these methods in terms of frequency.\n",
"* [LFAToolkit](https://jeremylt.github.io/LFAToolkit.jl/stable/examples/diffusion/#Diffusion-operator)"
]
},
{
"cell_type": "markdown",
"id": "83dfac85",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Multigrid as factorization"
]
},
{
"cell_type": "markdown",
"id": "8de0f9b4",
"metadata": {
"cell_style": "split"
},
"source": [
"We can interpret factorization as a particular multigrid or domain decomposition method.\n",
"\n",
"We can partition an SPD matrix as\n",
"$$A = \\begin{bmatrix} F & B^T \\\\ B & D \\end{bmatrix}$$\n",
"and define the preconditioner by the factorization\n",
"$$ M = \\begin{bmatrix} I & \\\\ B F^{-1} & I \\end{bmatrix}\n",
" \\begin{bmatrix} F & \\\\ & S \\end{bmatrix}\n",
" \\begin{bmatrix} I & F^{-1}B^T \\\\ & I \\end{bmatrix}\n",
"$$\n",
"where $S = D - B F^{-1} B^T$. $M$ has an inverse\n",
"$$\n",
" \\begin{bmatrix} I & -F^{-1}B^T \\\\ & I \\end{bmatrix}\n",
" \\begin{bmatrix} F^{-1} & \\\\ & S^{-1} \\end{bmatrix}\n",
" \\begin{bmatrix} I & \\\\ -B F^{-1} & I \\end{bmatrix} .\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "e4296193",
"metadata": {
"cell_style": "split"
},
"source": [
"Define the interpolation\n",
"$$ P_f = \\begin{bmatrix} I \\\\ 0 \\end{bmatrix}, \\quad P_c = \\begin{bmatrix} -F^{-1} B^T \\\\ I \\end{bmatrix} $$ and rewrite as\n",
"$$ M^{-1} = [P_f\\ P_c] \\begin{bmatrix} F^{-1} P_f^T \\\\ S^{-1} P_c^T \\end{bmatrix} = P_f F^{-1} P_f^T + P_c S^{-1} P_c^T.$$\n",
"The iteration matrix is thus\n",
"$$ I - M^{-1} A = I - P_f F^{-1} P_f^T J - P_c S^{-1} P_c^T A .$$"
]
},
{
"cell_type": "markdown",
"id": "9ba6e1af",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Permute into C-F split"
]
},
{
"cell_type": "code",
"execution_count": 67,
"id": "75b460cf",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"