diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..763513e
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+.ipynb_checkpoints
diff --git a/tutorials/linear-operators/LinearOperators.ipynb b/tutorials/linear-operators/LinearOperators.ipynb
new file mode 100644
index 0000000..713d642
--- /dev/null
+++ b/tutorials/linear-operators/LinearOperators.ipynb
@@ -0,0 +1,6380 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m environment at `~/Documents/streaming/notebooks/tutorials/linear-operators/Project.toml`\n",
+ "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
+ "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/Documents/streaming/notebooks/tutorials/linear-operators/Project.toml`\n",
+ "\u001b[90m [no changes]\u001b[39m\n",
+ "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/Documents/streaming/notebooks/tutorials/linear-operators/Manifest.toml`\n",
+ "\u001b[90m [no changes]\u001b[39m\n",
+ "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
+ "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Plots ─ v1.0.5\n",
+ "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/Documents/streaming/notebooks/tutorials/linear-operators/Project.toml`\n",
+ " \u001b[90m [91a5bcdd]\u001b[39m\u001b[92m + Plots v1.0.5\u001b[39m\n",
+ "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/Documents/streaming/notebooks/tutorials/linear-operators/Manifest.toml`\n",
+ " \u001b[90m [6e34b625]\u001b[39m\u001b[92m + Bzip2_jll v1.0.6+2\u001b[39m\n",
+ " \u001b[90m [3da002f7]\u001b[39m\u001b[92m + ColorTypes v0.10.0\u001b[39m\n",
+ " \u001b[90m [5ae59095]\u001b[39m\u001b[92m + Colors v0.12.0\u001b[39m\n",
+ " \u001b[90m [d38c429a]\u001b[39m\u001b[92m + Contour v0.5.2\u001b[39m\n",
+ " \u001b[90m [9a962f9c]\u001b[39m\u001b[92m + DataAPI v1.1.0\u001b[39m\n",
+ " \u001b[90m [864edb3b]\u001b[39m\u001b[92m + DataStructures v0.17.11\u001b[39m\n",
+ " \u001b[90m [c87230d0]\u001b[39m\u001b[92m + FFMPEG v0.3.0\u001b[39m\n",
+ " \u001b[90m [b22a6f82]\u001b[39m\u001b[92m + FFMPEG_jll v4.1.0+2\u001b[39m\n",
+ " \u001b[90m [53c48c17]\u001b[39m\u001b[92m + FixedPointNumbers v0.8.0\u001b[39m\n",
+ " \u001b[90m [d7e528f0]\u001b[39m\u001b[92m + FreeType2_jll v2.10.1+2\u001b[39m\n",
+ " \u001b[90m [559328eb]\u001b[39m\u001b[92m + FriBidi_jll v1.0.5+2\u001b[39m\n",
+ " \u001b[90m [28b8d3ca]\u001b[39m\u001b[92m + GR v0.48.0\u001b[39m\n",
+ " \u001b[90m [4d00f742]\u001b[39m\u001b[92m + GeometryTypes v0.8.2\u001b[39m\n",
+ " \u001b[90m [682c06a0]\u001b[39m\u001b[92m + JSON v0.21.0\u001b[39m\n",
+ " \u001b[90m [c1c5ebd0]\u001b[39m\u001b[92m + LAME_jll v3.100.0+0\u001b[39m\n",
+ " \u001b[90m [dd192d2f]\u001b[39m\u001b[92m + LibVPX_jll v1.8.1+1\u001b[39m\n",
+ " \u001b[90m [442fdcdd]\u001b[39m\u001b[92m + Measures v0.3.1\u001b[39m\n",
+ " \u001b[90m [e1d29d7a]\u001b[39m\u001b[92m + Missings v0.4.3\u001b[39m\n",
+ " \u001b[90m [77ba4419]\u001b[39m\u001b[92m + NaNMath v0.3.3\u001b[39m\n",
+ " \u001b[90m [e7412a2a]\u001b[39m\u001b[92m + Ogg_jll v1.3.3+0\u001b[39m\n",
+ " \u001b[90m [458c3c95]\u001b[39m\u001b[92m + OpenSSL_jll v1.1.1+2\u001b[39m\n",
+ " \u001b[90m [91d4177d]\u001b[39m\u001b[92m + Opus_jll v1.3.1+0\u001b[39m\n",
+ " \u001b[90m [bac558e1]\u001b[39m\u001b[92m + OrderedCollections v1.1.0\u001b[39m\n",
+ " \u001b[90m [69de0a69]\u001b[39m\u001b[92m + Parsers v1.0.1\u001b[39m\n",
+ " \u001b[90m [ccf2f8ad]\u001b[39m\u001b[92m + PlotThemes v1.0.3\u001b[39m\n",
+ " \u001b[90m [995b91a9]\u001b[39m\u001b[92m + PlotUtils v0.6.5\u001b[39m\n",
+ " \u001b[90m [91a5bcdd]\u001b[39m\u001b[92m + Plots v1.0.5\u001b[39m\n",
+ " \u001b[90m [3cdcf5f2]\u001b[39m\u001b[92m + RecipesBase v1.0.0\u001b[39m\n",
+ " \u001b[90m [189a3867]\u001b[39m\u001b[92m + Reexport v0.2.0\u001b[39m\n",
+ " \u001b[90m [ae029012]\u001b[39m\u001b[92m + Requires v1.0.1\u001b[39m\n",
+ " \u001b[90m [992d4aef]\u001b[39m\u001b[92m + Showoff v0.3.1\u001b[39m\n",
+ " \u001b[90m [a2af1166]\u001b[39m\u001b[92m + SortingAlgorithms v0.3.1\u001b[39m\n",
+ " \u001b[90m [90137ffa]\u001b[39m\u001b[92m + StaticArrays v0.12.1\u001b[39m\n",
+ " \u001b[90m [2913bbd2]\u001b[39m\u001b[92m + StatsBase v0.32.2\u001b[39m\n",
+ " \u001b[90m [83775a58]\u001b[39m\u001b[92m + Zlib_jll v1.2.11+9\u001b[39m\n",
+ " \u001b[90m [0ac62f75]\u001b[39m\u001b[92m + libass_jll v0.14.0+1\u001b[39m\n",
+ " \u001b[90m [f638f0a6]\u001b[39m\u001b[92m + libfdk_aac_jll v0.1.6+1\u001b[39m\n",
+ " \u001b[90m [f27f6e37]\u001b[39m\u001b[92m + libvorbis_jll v1.3.6+2\u001b[39m\n",
+ " \u001b[90m [1270edf5]\u001b[39m\u001b[92m + x264_jll v2019.5.25+1\u001b[39m\n",
+ " \u001b[90m [dfaa095f]\u001b[39m\u001b[92m + x265_jll v3.0.0+0\u001b[39m\n",
+ " \u001b[90m [2a0f44e3]\u001b[39m\u001b[92m + Base64 \u001b[39m\n",
+ " \u001b[90m [ade2ca70]\u001b[39m\u001b[92m + Dates \u001b[39m\n",
+ " \u001b[90m [8bb1440f]\u001b[39m\u001b[92m + DelimitedFiles \u001b[39m\n",
+ " \u001b[90m [8ba89e20]\u001b[39m\u001b[92m + Distributed \u001b[39m\n",
+ " \u001b[90m [b77e0a4c]\u001b[39m\u001b[92m + InteractiveUtils \u001b[39m\n",
+ " \u001b[90m [76f85450]\u001b[39m\u001b[92m + LibGit2 \u001b[39m\n",
+ " \u001b[90m [56ddb016]\u001b[39m\u001b[92m + Logging \u001b[39m\n",
+ " \u001b[90m [d6f4376e]\u001b[39m\u001b[92m + Markdown \u001b[39m\n",
+ " \u001b[90m [a63ad114]\u001b[39m\u001b[92m + Mmap \u001b[39m\n",
+ " \u001b[90m [44cfe95a]\u001b[39m\u001b[92m + Pkg \u001b[39m\n",
+ " \u001b[90m [3fa0cd96]\u001b[39m\u001b[92m + REPL \u001b[39m\n",
+ " \u001b[90m [ea8e919c]\u001b[39m\u001b[92m + SHA \u001b[39m\n",
+ " \u001b[90m [6462fe0b]\u001b[39m\u001b[92m + Sockets \u001b[39m\n",
+ " \u001b[90m [10745b16]\u001b[39m\u001b[92m + Statistics \u001b[39m\n",
+ " \u001b[90m [8dfed614]\u001b[39m\u001b[92m + Test \u001b[39m\n",
+ " \u001b[90m [cf7118a7]\u001b[39m\u001b[92m + UUIDs \u001b[39m\n",
+ "\u001b[32m\u001b[1m Building\u001b[22m\u001b[39m Plots → `~/.julia/packages/Plots/GzWxB/deps/build.log`\n",
+ "┌ Error: Error building `Plots`: \n",
+ "│ ERROR: LoadError: ArgumentError: Package OhMyREPL not found in current path:\n",
+ "│ - Run `import Pkg; Pkg.add(\"OhMyREPL\")` to install the OhMyREPL package.\n",
+ "│ \n",
+ "│ Stacktrace:\n",
+ "│ [1] require(::Module, ::Symbol) at ./loading.jl:892\n",
+ "│ [2] include(::Module, ::String) at ./Base.jl:377\n",
+ "│ [3] include_ifexists at ./client.jl:212 [inlined]\n",
+ "│ [4] load_julia_startup() at ./client.jl:320\n",
+ "│ [5] exec_options(::Base.JLOptions) at ./client.jl:259\n",
+ "│ [6] _start() at ./client.jl:484\n",
+ "│ in expression starting at /home/abel/.julia/config/startup.jl:1\n",
+ "└ @ Pkg.Operations /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/Pkg/src/Operations.jl:892\n",
+ "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
+ "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/Documents/streaming/notebooks/tutorials/linear-operators/Project.toml`\n",
+ " \u001b[90m [28b8d3ca]\u001b[39m\u001b[92m + GR v0.48.0\u001b[39m\n",
+ "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/Documents/streaming/notebooks/tutorials/linear-operators/Manifest.toml`\n",
+ "\u001b[90m [no changes]\u001b[39m\n"
+ ]
+ }
+ ],
+ "source": [
+ "using Pkg\n",
+ "pkg\"activate .\"\n",
+ "pkg\"add LinearOperators\"\n",
+ "pkg\"add Plots\"\n",
+ "pkg\"add GR\""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\u001b[32m\u001b[1mStatus\u001b[22m\u001b[39m `~/Documents/streaming/notebooks/tutorials/linear-operators/Project.toml`\n",
+ " \u001b[90m [28b8d3ca]\u001b[39m\u001b[37m GR v0.48.0\u001b[39m\n",
+ " \u001b[90m [5c8ed15e]\u001b[39m\u001b[37m LinearOperators v1.0.1\u001b[39m\n",
+ " \u001b[90m [91a5bcdd]\u001b[39m\u001b[37m Plots v1.0.5\u001b[39m\n"
+ ]
+ }
+ ],
+ "source": [
+ "pkg\"status\""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "┌ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]\n",
+ "└ @ Base loading.jl:1260\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n"
+ ]
+ },
+ "execution_count": 7,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using Plots\n",
+ "gr(size=(400,300))\n",
+ "plot(rand(3), rand(3))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# LinearOperators"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "- v -> Av\n",
+ "- v -> Aᵀv\n",
+ "- v -> A*v"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "2-element Array{Float64,1}:\n",
+ " 0.7\n",
+ " 1.9"
+ ]
+ },
+ "execution_count": 10,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using LinearOperators\n",
+ "\n",
+ "prod(v) = [v[1] + v[2]; 2v[1] + 3v[2]] # [1 1; 2 3] * [v[1]; v[2]]\n",
+ "tprod(v) = [v[1] + 2v[2]; v[1] + 3v[2]]\n",
+ "\n",
+ "A = LinearOperator(Float64, 2, 2, false, false, prod, tprod, tprod)\n",
+ "\n",
+ "A * [0.2; 0.5]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "5×2 Array{Float64,2}:\n",
+ " 0.374079 0.492553\n",
+ " 0.187629 0.600294\n",
+ " 0.183658 0.398463\n",
+ " 0.736983 0.434411\n",
+ " 0.779344 0.406829"
+ ]
+ },
+ "execution_count": 11,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A = rand(5, 2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Linear operator\n",
+ " nrow: 5\n",
+ " ncol: 2\n",
+ " eltype: Float64\n",
+ " symmetric: false\n",
+ " hermitian: false\n",
+ " nprod: 0\n",
+ " ntprod: 0\n",
+ " nctprod: 0\n",
+ "\n"
+ ]
+ },
+ "execution_count": 12,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "opA = LinearOperator(A)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " 4.008001 seconds (2 allocations: 190.735 MiB)\n"
+ ]
+ }
+ ],
+ "source": [
+ "A = rand(5000, 5000);\n",
+ "B = rand(5000, 5000);\n",
+ "@time C = A * B;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " 0.000010 seconds (4 allocations: 176 bytes)\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "Linear operator\n",
+ " nrow: 5000\n",
+ " ncol: 5000\n",
+ " eltype: Float64\n",
+ " symmetric: false\n",
+ " hermitian: false\n",
+ " nprod: 0\n",
+ " ntprod: 0\n",
+ " nctprod: 0\n",
+ "\n"
+ ]
+ },
+ "execution_count": 25,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "opA = LinearOperator(A)\n",
+ "opB = LinearOperator(B)\n",
+ "@time opA * opB"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " 4.006204 seconds (4 allocations: 190.773 MiB)\n",
+ " 0.041024 seconds (4 allocations: 78.281 KiB)\n",
+ " 0.043637 seconds (8 allocations: 78.453 KiB)\n",
+ " 0.051469 seconds (4 allocations: 78.281 KiB)\n"
+ ]
+ }
+ ],
+ "source": [
+ "v = rand(5000)\n",
+ "\n",
+ "@time (A * B) * v\n",
+ "@time A * (B * v)\n",
+ "@time (opA * opB) * v\n",
+ "@time opA * (opB * v);"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**Should I access an operator's element?**\n",
+ "\n",
+ "No"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 31,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "2×2 Array{Float64,2}:\n",
+ " 1.0 1.0\n",
+ " 2.0 3.0"
+ ]
+ },
+ "execution_count": 31,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "prod(v) = [v[1] + v[2]; 2v[1] + 3v[2]] # [1 1; 2 3] * [v[1]; v[2]]\n",
+ "tprod(v) = [v[1] + 2v[2]; v[1] + 3v[2]]\n",
+ "\n",
+ "A = LinearOperator(Float64, 2, 2, false, false, prod, tprod, tprod)\n",
+ "Matrix(A)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Linear operator\n",
+ " nrow: 1\n",
+ " ncol: 1\n",
+ " eltype: Float64\n",
+ " symmetric: false\n",
+ " hermitian: false\n",
+ " nprod: 0\n",
+ " ntprod: 0\n",
+ " nctprod: 0\n",
+ "\n"
+ ]
+ },
+ "execution_count": 34,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A[2,2] # Compute Z' * A * P"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 35,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "1×1 Array{Float64,2}:\n",
+ " 3.0"
+ ]
+ },
+ "execution_count": 35,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "Matrix(A[2,2])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 36,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Using opExtension, opRestriction"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 37,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Linear operator\n",
+ " nrow: 5\n",
+ " ncol: 2\n",
+ " eltype: Float64\n",
+ " symmetric: false\n",
+ " hermitian: false\n",
+ " nprod: 0\n",
+ " ntprod: 0\n",
+ " nctprod: 0\n",
+ "\n"
+ ]
+ },
+ "execution_count": 37,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A = rand(5, 2)\n",
+ "opA = LinearOperator(A)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 39,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "true"
+ ]
+ },
+ "execution_count": 39,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A * ones(2) == opA * ones(2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 41,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "2-element Array{Float64,1}:\n",
+ " 2.9377342310451233\n",
+ " 1.1810979848574872"
+ ]
+ },
+ "execution_count": 41,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "opA' * ones(5)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 42,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
+ "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/Documents/streaming/notebooks/tutorials/linear-operators/Project.toml`\n",
+ " \u001b[90m [7a1cc6ca]\u001b[39m\u001b[92m + FFTW v1.2.0\u001b[39m\n",
+ "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `~/Documents/streaming/notebooks/tutorials/linear-operators/Manifest.toml`\n",
+ " \u001b[90m [621f4979]\u001b[39m\u001b[92m + AbstractFFTs v0.5.0\u001b[39m\n",
+ " \u001b[90m [7a1cc6ca]\u001b[39m\u001b[92m + FFTW v1.2.0\u001b[39m\n",
+ " \u001b[90m [f5851436]\u001b[39m\u001b[92m + FFTW_jll v3.3.9+5\u001b[39m\n",
+ " \u001b[90m [1d5cc7b8]\u001b[39m\u001b[92m + IntelOpenMP_jll v2018.0.3+0\u001b[39m\n",
+ " \u001b[90m [856f044c]\u001b[39m\u001b[92m + MKL_jll v2019.0.117+2\u001b[39m\n"
+ ]
+ }
+ ],
+ "source": [
+ "pkg\"add FFTW\""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 43,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "16-element Array{Complex{Float64},1}:\n",
+ " 8.946161193394072 + 8.045984280505243im\n",
+ " 0.1742357191505106 + 0.3086013649576881im\n",
+ " -0.7998530430635715 + 0.1115711495074454im\n",
+ " -0.571398395752621 + 0.4544987711531042im\n",
+ " -0.2704035645753835 - 1.9213665753970317im\n",
+ " 2.324686823971292 + 0.5140642148647618im\n",
+ " -0.007250822847700414 - 0.846160591182173im\n",
+ " 1.176282988707905 - 0.9399125991116004im\n",
+ " -0.9915573471002972 - 1.2934188837560376im\n",
+ " -1.287207915709439 + 1.4402476202412147im\n",
+ " -0.08436271058601133 + 2.053446184649336im\n",
+ " 1.0016829875186235 + 1.0366293915288958im\n",
+ " -1.043324233945478 - 0.4355521279299075im\n",
+ " 1.54306565608972 + 2.213251096982505im\n",
+ " 0.23391290302568207 + 0.6780221260155858im\n",
+ " 0.20842004995391392 - 0.5432979505790254im"
+ ]
+ },
+ "execution_count": 43,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using FFTW, LinearAlgebra\n",
+ "\n",
+ "A = LinearOperator(16, 16, false, false, fft, nothing, ifft)\n",
+ "\n",
+ "v = rand(16) + im * rand(16)\n",
+ "A * v"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 44,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "16-element Array{Complex{Float64},1}:\n",
+ " 0.5591350745871295 + 0.5028740175315777im\n",
+ " 0.01302625312211962 - 0.03395612191118909im\n",
+ " 0.01461955643910513 + 0.04237638287597411im\n",
+ " 0.0964416035056075 + 0.13832819356140658im\n",
+ " -0.06520776462159238 - 0.02722200799561922im\n",
+ " 0.06260518671991397 + 0.06478933697055599im\n",
+ " -0.005272669411625708 + 0.1283403865405835im\n",
+ " -0.08045049473183993 + 0.09001547626507592im\n",
+ " -0.061972334193768575 - 0.08083868023475235im\n",
+ " 0.07351768679424406 - 0.058744537444475026im\n",
+ " -0.00045317642798127587 - 0.052885036948885814im\n",
+ " 0.14529292649820574 + 0.032129013429047615im\n",
+ " -0.01690022278596147 - 0.12008541096231448im\n",
+ " -0.03571239973453881 + 0.028406173197069014im\n",
+ " -0.04999081519147322 + 0.0069731968442153375im\n",
+ " 0.010889732446906912 + 0.019287585309855505im"
+ ]
+ },
+ "execution_count": 44,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A' * v"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 49,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "true"
+ ]
+ },
+ "execution_count": 49,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A' * (A * v) ≈ v"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 47,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "4.0697662735997785e-16"
+ ]
+ },
+ "execution_count": 47,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "norm(A' * (A * v) - v)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 50,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "true"
+ ]
+ },
+ "execution_count": 50,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "Matrix(A' * A) ≈ Matrix(I, 16, 16)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 51,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Linear operator\n",
+ " nrow: 3\n",
+ " ncol: 3\n",
+ " eltype: Float64\n",
+ " symmetric: false\n",
+ " hermitian: false\n",
+ " nprod: 0\n",
+ " ntprod: 0\n",
+ " nctprod: 0\n",
+ "\n"
+ ]
+ },
+ "execution_count": 51,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "opA = LinearOperator(rand(3, 3))\n",
+ "opB = LinearOperator(rand(3, 3))\n",
+ "\n",
+ "opA + opB"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 56,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Linear operator\n",
+ " nrow: 3\n",
+ " ncol: 3\n",
+ " eltype: Float64\n",
+ " symmetric: false\n",
+ " hermitian: false\n",
+ " nprod: 0\n",
+ " ntprod: 0\n",
+ " nctprod: 0\n",
+ "\n"
+ ]
+ },
+ "execution_count": 56,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "transpose(conj(opA * 2 - 3 * opB')')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**Misuse of nonlinear in LinearOperators**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 84,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "2×2 Array{Float64,2}:\n",
+ " 1.0 1.0\n",
+ " 0.0 0.0"
+ ]
+ },
+ "execution_count": 84,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "prod(v) = [v[1]^2 + v[2]^2; v[1] * v[2]]\n",
+ "tprod(v) = [v[1] + 2v[2]; v[1] + 3v[2]]\n",
+ "\n",
+ "A = LinearOperator(Float64, 2, 2, false, false, prod, tprod, tprod)\n",
+ "Matrix(A)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Large problems in optimization"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 57,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Hv (generic function with 1 method)"
+ ]
+ },
+ "execution_count": 57,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using ForwardDiff\n",
+ "\n",
+ "f(x) = (x[1] - 1)^2 + 100 * (x[2] - x[1]^2)^2\n",
+ "∇f(x) = ForwardDiff.gradient(f, x)\n",
+ "Hv(x, v) = ForwardDiff.derivative(t -> ∇f(x + t * v), 0.0)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 58,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "(24.199999999999996, [-215.59999999999997, -87.99999999999999])"
+ ]
+ },
+ "execution_count": 58,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "x = [-1.2; 1.0]\n",
+ "\n",
+ "f(x), ∇f(x)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 59,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "2-element Array{Float64,1}:\n",
+ " 1809.9999999999998\n",
+ " 680.0"
+ ]
+ },
+ "execution_count": 59,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "v = ones(2)\n",
+ "Hv(x, v)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 61,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "2-element Array{Float64,1}:\n",
+ " 1810.0\n",
+ " 680.0"
+ ]
+ },
+ "execution_count": 61,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "ForwardDiff.hessian(f, x) * v"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 62,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "opHatx (generic function with 1 method)"
+ ]
+ },
+ "execution_count": 62,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "opHatx(x) = LinearOperator(Float64, 2, 2, true, true, v -> Hv(x, v))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 63,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Linear operator\n",
+ " nrow: 2\n",
+ " ncol: 2\n",
+ " eltype: Float64\n",
+ " symmetric: true\n",
+ " hermitian: true\n",
+ " nprod: 0\n",
+ " ntprod: 0\n",
+ " nctprod: 0\n",
+ "\n"
+ ]
+ },
+ "execution_count": 63,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "opH = opHatx(x)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 64,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "2-element Array{Float64,1}:\n",
+ " 1809.9999999999998\n",
+ " 680.0"
+ ]
+ },
+ "execution_count": 64,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "opH * v"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 66,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "([0.9999999999999999, 0.9999999999814724], [7.410960511933238e-9, -3.7054803669889225e-9])"
+ ]
+ },
+ "execution_count": 66,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using Krylov # Implement Conjugate Gradient - more on future tutorials\n",
+ "\n",
+ "while norm(∇f(x)) > 1e-6\n",
+ " opH = opHatx(x)\n",
+ " d, ks = Krylov.cg(opH, -∇f(x)) # Solution of the system H(x)d = -∇f(x)\n",
+ " x += d\n",
+ "end\n",
+ "\n",
+ "x, ∇f(x)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Preallocated operators"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 71,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m, n = 50, 30\n",
+ "A = rand(50, 30)\n",
+ "op1 = LinearOperator(A)\n",
+ "op2 = PreallocatedLinearOperator(A) # create a vector for storing Av\n",
+ "v = rand(n);"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 73,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "49600"
+ ]
+ },
+ "execution_count": 73,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "al = @allocated for i = 1:100\n",
+ " op1 * v\n",
+ "end\n",
+ "al"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 74,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "0"
+ ]
+ },
+ "execution_count": 74,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "al = @allocated for i = 1:100\n",
+ " op2 * v\n",
+ "end\n",
+ "al"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 83,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Aones = op2 * ones(n)\n",
+ "Atwos = op2 * (2 * ones(n)); # Changes the internal memory"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 85,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "true"
+ ]
+ },
+ "execution_count": 85,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "Aones === Atwos"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Inverse operators"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 93,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Linear operator\n",
+ " nrow: 5\n",
+ " ncol: 5\n",
+ " eltype: Float64\n",
+ " symmetric: true\n",
+ " hermitian: true\n",
+ " nprod: 0\n",
+ " ntprod: 0\n",
+ " nctprod: 0\n",
+ "\n"
+ ]
+ },
+ "execution_count": 93,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A = rand(5, 5)\n",
+ "A = A' * A + I\n",
+ "b = A * ones(5)\n",
+ "\n",
+ "opA = LinearOperator(A)\n",
+ "opAinv = opInverse(A)\n",
+ "opAchol = opCholesky(A)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 91,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "5-element Array{Float64,1}:\n",
+ " 1.0000000000000004\n",
+ " 1.0000000000000007\n",
+ " 0.9999999999999988\n",
+ " 1.0000000000000004\n",
+ " 1.0000000000000004"
+ ]
+ },
+ "execution_count": 91,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "opAinv * b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 92,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "5-element Array{Float64,1}:\n",
+ " 1.0000000000000004\n",
+ " 1.0\n",
+ " 0.999999999999999\n",
+ " 1.0000000000000002\n",
+ " 1.0000000000000007"
+ ]
+ },
+ "execution_count": 92,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "opAchol * b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 94,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "([0.9999999999999422, 0.9999999999999255, 0.9999999999999524, 0.9999999999999439, 0.9999999999999344], \n",
+ "Simple stats\n",
+ " solved: true\n",
+ " inconsistent: false\n",
+ " residuals: [ 1.8e+01 3.8e-01 9.8e-02 6.8e-03 1.3e-03 1.1e-12 ]\n",
+ " Aresiduals: []\n",
+ " status: solution good enough given atol and rtol\n",
+ ")"
+ ]
+ },
+ "execution_count": 94,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "Krylov.cg(opA, b)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 95,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "([1.0000000000000007, 1.0000000000000002, 0.9999999999999987, 1.0000000000000004, 0.9999999999999997], \n",
+ "Simple stats\n",
+ " solved: true\n",
+ " inconsistent: false\n",
+ " residuals: [ 6.4e+00 6.7e-16 ]\n",
+ " Aresiduals: []\n",
+ " status: solution good enough given atol and rtol\n",
+ ")"
+ ]
+ },
+ "execution_count": 95,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "Krylov.cg(opA, b, M=opAchol)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## BFGS Operators"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 98,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Linear operator\n",
+ " nrow: 5\n",
+ " ncol: 5\n",
+ " eltype: Float64\n",
+ " symmetric: true\n",
+ " hermitian: true\n",
+ " nprod: 0\n",
+ " ntprod: 0\n",
+ " nctprod: 0\n",
+ "\n"
+ ]
+ },
+ "execution_count": 98,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "B = LBFGSOperator(5, scaling=false)\n",
+ "H = InverseLBFGSOperator(5, scaling=false)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 99,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "5-element Array{Float64,1}:\n",
+ " 1.0\n",
+ " 1.0\n",
+ " 1.0\n",
+ " 1.0\n",
+ " 1.0"
+ ]
+ },
+ "execution_count": 99,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "B * ones(5)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 100,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Linear operator\n",
+ " nrow: 5\n",
+ " ncol: 5\n",
+ " eltype: Float64\n",
+ " symmetric: true\n",
+ " hermitian: true\n",
+ " nprod: 1\n",
+ " ntprod: 0\n",
+ " nctprod: 0\n",
+ "\n"
+ ]
+ },
+ "execution_count": 100,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "s = rand(5)\n",
+ "y = rand(5)\n",
+ "push!(B, s, y) # s = xₖ₊₁ - xₖ, y = ∇f(xₖ₊₁) - ∇f(xₖ)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 101,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "5-element Array{Float64,1}:\n",
+ " 0.30087051250127916\n",
+ " 0.2741514306687156\n",
+ " 0.15426237166249435\n",
+ " 1.4713507119542513\n",
+ " 2.1279737389023086"
+ ]
+ },
+ "execution_count": 101,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "B * ones(5)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 102,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Linear operator\n",
+ " nrow: 5\n",
+ " ncol: 5\n",
+ " eltype: Float64\n",
+ " symmetric: true\n",
+ " hermitian: true\n",
+ " nprod: 0\n",
+ " ntprod: 0\n",
+ " nctprod: 0\n",
+ "\n"
+ ]
+ },
+ "execution_count": 102,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "push!(H, s, y)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 103,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "5-element Array{Float64,1}:\n",
+ " 1.0000000000000002\n",
+ " 1.0000000000000002\n",
+ " 1.0000000000000002\n",
+ " 1.0000000000000004\n",
+ " 1.0"
+ ]
+ },
+ "execution_count": 103,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "H * (B * ones(5))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 115,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "([1.0000000051838849, 1.0000000106534566], [-1.0390695385617519e-7, 5.713736150880777e-8], 677)"
+ ]
+ },
+ "execution_count": 115,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "H = InverseLBFGSOperator(2, 2) # nvar, mem\n",
+ "\n",
+ "iter = 0\n",
+ "x = [-1.2; 1.0]\n",
+ "gx = ∇f(x)\n",
+ "while !(norm(gx) < 1e-6 || iter > 10_000)\n",
+ " #d = -gx\n",
+ " d = -H * gx\n",
+ " slope = dot(d, gx)\n",
+ " t = 1.0\n",
+ " while !(f(x + t * d) < f(x) + 1e-4 * t * slope)\n",
+ " t /= 2\n",
+ " end\n",
+ " x += t * d\n",
+ " s = t * d\n",
+ " gxp = ∇f(x)\n",
+ " y = gxp - gx\n",
+ " gx = gxp\n",
+ " push!(H, s, y)\n",
+ " iter += 1\n",
+ "end\n",
+ "\n",
+ "x, ∇f(x), iter"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 119,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "([1.0000000337283783, 1.0000000676093637], [6.4143834083636564e-9, 3.0521185578891163e-8], 100)"
+ ]
+ },
+ "execution_count": 119,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "B = LBFGSOperator(2, 2) # nvar, mem\n",
+ "\n",
+ "iter = 0\n",
+ "x = [-1.2; 1.0]\n",
+ "gx = ∇f(x)\n",
+ "Δ = 1.0\n",
+ "while !(norm(gx) < 1e-6 || iter > 10_000)\n",
+ " d, ks = Krylov.cg(B, -gx, radius=Δ)\n",
+ " # Update Δ\n",
+ " x += d\n",
+ " s = d\n",
+ " gxp = ∇f(x)\n",
+ " y = gxp - gx\n",
+ " gx = gxp\n",
+ " push!(B, s, y)\n",
+ " iter += 1\n",
+ "end\n",
+ "\n",
+ "x, ∇f(x), iter"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Application: Heat Equation\n",
+ "\n",
+ "https://juliasmoothoptimizers.github.io/JSOTutorials.jl/linear-operators/introduction-to-linear-operators/introduction-to-linear-operators.html"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "a plate with length L, with no heat exchange on the boundaries."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\\begin{align}\n",
+ "\\begin{bmatrix}\n",
+ "T & 2D \\\\\n",
+ "D & T & D \\\\\n",
+ " & D & T & D \\\\\n",
+ " & & 2D & T\n",
+ "\\end{bmatrix}\n",
+ "\\end{align}\n",
+ "where\n",
+ "\\begin{align}\n",
+ "T = \\begin{bmatrix}\n",
+ "\\kappa & 2\\gamma \\\\\n",
+ "\\gamma & \\kappa & \\gamma \\\\\n",
+ " & \\gamma & \\kappa & \\gamma \\\\\n",
+ " & & 2\\gamma & \\kappa\n",
+ "\\end{bmatrix}\n",
+ "\\end{align}\n",
+ "and $D = \\text{diag}(\\gamma,\\dots,\\gamma)$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 121,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "HeatEquationOp (generic function with 1 method)"
+ ]
+ },
+ "execution_count": 121,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "function HeatEquationOp(L, # length\n",
+ " m, # number of points in the grid\n",
+ " δ, # The time step\n",
+ " α, # heat coefficient\n",
+ " )\n",
+ " h = L / (m - 1) # spatial step\n",
+ " γ = α * δ / h^2\n",
+ " κ = 1 - 4γ\n",
+ " \n",
+ " Tprod(v) = [κ * v[1] + 2γ * v[2];\n",
+ " [γ * v[i-1] + κ * v[i] + γ * v[i+1] for i = 2:m-1]\n",
+ " κ * v[m] + 2γ * v[m-1]]\n",
+ " \n",
+ " T = LinearOperator(Float64, m, m, false, false, Tprod)\n",
+ " D = opEye(m) * γ\n",
+ " \n",
+ " function prod(v)\n",
+ " Hv = zeros(m^2)\n",
+ " Hv[1:m] .= [T 2D] * v[1:2m]\n",
+ " for i = 2:m-1\n",
+ " Hv[(i-1)*m+1:i*m] .= [D T D] * v[(i-2)*m+1:(i+1)*m]\n",
+ " end\n",
+ " Hv[end-m+1:end] .= [2D T] * v[end-2m+1:end]\n",
+ " return Hv\n",
+ " end\n",
+ " \n",
+ " return LinearOperator(Float64, m^2, m^2, false, false, prod)\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 127,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "9×9 Array{Float64,2}:\n",
+ " 0.84 0.08 0.0 0.08 0.0 0.0 0.0 0.0 0.0\n",
+ " 0.04 0.84 0.04 0.0 0.08 0.0 0.0 0.0 0.0\n",
+ " 0.0 0.08 0.84 0.0 0.0 0.08 0.0 0.0 0.0\n",
+ " 0.04 0.0 0.0 0.84 0.08 0.0 0.04 0.0 0.0\n",
+ " 0.0 0.04 0.0 0.04 0.84 0.04 0.0 0.04 0.0\n",
+ " 0.0 0.0 0.04 0.0 0.08 0.84 0.0 0.0 0.04\n",
+ " 0.0 0.0 0.0 0.08 0.0 0.0 0.84 0.08 0.0\n",
+ " 0.0 0.0 0.0 0.0 0.08 0.0 0.04 0.84 0.04\n",
+ " 0.0 0.0 0.0 0.0 0.0 0.08 0.0 0.08 0.84"
+ ]
+ },
+ "execution_count": 127,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A = HeatEquationOp(1.0, 3, 0.1, 0.1)\n",
+ "Matrix(A)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 131,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "