{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Julia Version 1.6.2\n",
"Commit 1b93d53fc4 (2021-07-14 15:36 UTC)\n",
"Platform Info:\n",
" OS: macOS (x86_64-apple-darwin18.7.0)\n",
" CPU: Intel(R) Core(TM) i5-8279U CPU @ 2.40GHz\n",
" WORD_SIZE: 64\n",
" LIBM: libopenlibm\n",
" LLVM: libLLVM-11.0.1 (ORCJIT, skylake)\n"
]
}
],
"source": [
"versioninfo()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m environment at `~/Dropbox/class/M1399.000200/2021/M1399_000200-2021fall/Project.toml`\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[32m\u001b[1m Status\u001b[22m\u001b[39m `~/Dropbox/class/M1399.000200/2021/M1399_000200-2021fall/Project.toml`\n",
" \u001b[90m [7d9fca2a] \u001b[39mArpack v0.4.0\n",
" \u001b[90m [6e4b80f9] \u001b[39mBenchmarkTools v1.2.0\n",
" \u001b[90m [1e616198] \u001b[39mCOSMO v0.8.1\n",
" \u001b[90m [f65535da] \u001b[39mConvex v0.14.16\n",
" \u001b[90m [a93c6f00] \u001b[39mDataFrames v1.2.2\n",
" \u001b[90m [31a5f54b] \u001b[39mDebugger v0.6.8\n",
" \u001b[90m [31c24e10] \u001b[39mDistributions v0.24.18\n",
" \u001b[90m [e2685f51] \u001b[39mECOS v0.12.3\n",
" \u001b[90m [f6369f11] \u001b[39mForwardDiff v0.10.21\n",
" \u001b[90m [28b8d3ca] \u001b[39mGR v0.61.0\n",
" \u001b[90m [c91e804a] \u001b[39mGadfly v1.3.3\n",
" \u001b[90m [bd48cda9] \u001b[39mGraphRecipes v0.5.7\n",
" \u001b[90m [f67ccb44] \u001b[39mHDF5 v0.14.3\n",
" \u001b[90m [82e4d734] \u001b[39mImageIO v0.5.8\n",
" \u001b[90m [6218d12a] \u001b[39mImageMagick v1.2.1\n",
" \u001b[90m [916415d5] \u001b[39mImages v0.24.1\n",
" \u001b[90m [b6b21f68] \u001b[39mIpopt v0.7.0\n",
" \u001b[90m [42fd0dbc] \u001b[39mIterativeSolvers v0.9.1\n",
" \u001b[90m [4076af6c] \u001b[39mJuMP v0.21.10\n",
" \u001b[90m [b51810bb] \u001b[39mMatrixDepot v1.0.4\n",
" \u001b[90m [1ec41992] \u001b[39mMosekTools v0.9.4\n",
" \u001b[90m [76087f3c] \u001b[39mNLopt v0.6.3\n",
" \u001b[90m [47be7bcc] \u001b[39mORCA v0.5.0\n",
" \u001b[90m [a03496cd] \u001b[39mPlotlyBase v0.4.3\n",
" \u001b[90m [f0f68f2c] \u001b[39mPlotlyJS v0.15.0\n",
" \u001b[90m [91a5bcdd] \u001b[39mPlots v1.22.6\n",
" \u001b[90m [438e738f] \u001b[39mPyCall v1.92.3\n",
" \u001b[90m [d330b81b] \u001b[39mPyPlot v2.10.0\n",
" \u001b[90m [dca85d43] \u001b[39mQuartzImageIO v0.7.3\n",
" \u001b[90m [6f49c342] \u001b[39mRCall v0.13.12\n",
" \u001b[90m [ce6b1742] \u001b[39mRDatasets v0.7.5\n",
" \u001b[90m [c946c3f1] \u001b[39mSCS v0.7.1\n",
" \u001b[90m [276daf66] \u001b[39mSpecialFunctions v1.7.0\n",
" \u001b[90m [2913bbd2] \u001b[39mStatsBase v0.33.12\n",
" \u001b[90m [f3b207a7] \u001b[39mStatsPlots v0.14.28\n",
" \u001b[90m [b8865327] \u001b[39mUnicodePlots v2.4.6\n",
" \u001b[90m [0f1e0344] \u001b[39mWebIO v0.8.16\n",
" \u001b[90m [8f399da3] \u001b[39mLibdl\n",
" \u001b[90m [2f01184e] \u001b[39mSparseArrays\n",
" \u001b[90m [10745b16] \u001b[39mStatistics\n"
]
}
],
"source": [
"using Pkg\n",
"Pkg.activate(\"../..\")\n",
"Pkg.status()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Iterative Methods for Solving Linear Equations\n",
"\n",
"## Introduction\n",
"\n",
"So far we have considered direct methods for solving linear equations. \n",
"\n",
"* **Direct methods** (flops fixed _a priori_) vs **iterative methods**:\n",
" - Direct method (GE/LU, Cholesky, QR, SVD): accurate. good for dense, small to moderate sized $\\mathbf{A}$ \n",
" - Iterative methods (Jacobi, Gauss-Seidal, SOR, conjugate-gradient, GMRES): accuracy depends on number of iterations. good for large, sparse, structured linear system, parallel computing, warm start (reasonable accuracy after, say, 100 iterations)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Motivation: PageRank \n",
"\n",
"\n",
"\n",
"Source: [Wikepedia](https://en.wikipedia.org/wiki/PageRank)\n",
"\n",
"* $\\mathbf{A} \\in \\{0,1\\}^{n \\times n}$ the connectivity matrix of webpages with entries\n",
"$$\n",
"\\begin{eqnarray*}\n",
"\ta_{ij} = \\begin{cases}\n",
"\t1 & \\text{if page $i$ links to page $j$} \\\\\n",
"\t0 & \\text{otherwise}\n",
"\t\\end{cases}. \n",
"\\end{eqnarray*}\n",
"$$\n",
"$n \\approx 10^9$ in May 2017.\n",
"\n",
"* $r_i = \\sum_j a_{ij}$ is the *out-degree* of page $i$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* [Larry Page](https://en.wikipedia.org/wiki/PageRank) imagined a random surfer wandering on internet according to following rules:\n",
" - From a page $i$ with $r_i>0$\n",
" * with probability $p$, (s)he randomly chooses a link $j$ from page $i$ (uniformly) and follows that link to the next page \n",
" * with probability $1-p$, (s)he randomly chooses one page from the set of all $n$ pages (uniformly) and proceeds to that page \n",
" - From a page $i$ with $r_i=0$ (a dangling page), (s)he randomly chooses one page from the set of all $n$ pages (uniformly) and proceeds to that page\n",
" \n",
"* The process defines an $n$-state Markov chain, where each state corresponds to each page. \n",
"$$\n",
" p_{ij} = (1-p)\\frac{1}{n} + p\\frac{a_{ij}}{r_i}\n",
"$$\n",
"with interpretation $a_{ij}/r_i = 1/n$ if $r_i = 0$.\n",
"\n",
"* Stationary distribution of this Markov chain gives the score (long term probability of visit) of each page.\n",
"\n",
"* Stationary distribution can be obtained as the top *left* eigenvector of the transition matrix $\\mathbf{P}=(p_{ij})$ corresponding to eigenvalue 1. \n",
"$$\n",
" \\mathbf{x}^T\\mathbf{P} = \\mathbf{x}^T.\n",
"$$\n",
"Equivalently it can be cast as a linear equation.\n",
"$$\n",
" (\\mathbf{I} - \\mathbf{P}^T) \\mathbf{x} = \\mathbf{0}.\n",
"$$\n",
"\n",
"* You've got to solve a linear equation with $10^9$ variables!\n",
"\n",
"* GE/LU will take $2 \\times (10^9)^3/3/10^{12} \\approx 6.66 \\times 10^{14}$ seconds $\\approx 20$ million years on a tera-flop supercomputer!\n",
"\n",
"* Iterative methods come to the rescue."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Splitting and fixed point iteration\n",
"\n",
"* The key idea of iterative method for solving $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$ is to **split** the matrix $\\mathbf{A}$ so that\n",
"$$\n",
" \\mathbf{A} = \\mathbf{M} - \\mathbf{K}\n",
"$$\n",
"where $\\mathbf{M}$ is invertible and easy to invert.\n",
"\n",
"* Then $\\mathbf{A}\\mathbf{x} = \\mathbf{M}\\mathbf{x} - \\mathbf{K}\\mathbf{x} = \\mathbf{b}$ or\n",
"$$\n",
" \\mathbf{x} = \\mathbf{M}^{-1}\\mathbf{K}\\mathbf{x} - \\mathbf{M}^{-1}\\mathbf{b}\n",
" .\n",
"$$\n",
"Thus a solution to $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$ is a fixed point of iteration\n",
"$$\n",
" \\mathbf{x}^{(t+1)} = \\mathbf{M}^{-1}\\mathbf{K}\\mathbf{x}^{(t)} - \\mathbf{M}^{-1}\\mathbf{b}\n",
" = \\mathbf{R}\\mathbf{x}^{(k)} - \\mathbf{c}\n",
" .\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Under a suitable choice of $\\mathbf{R}$, i.e., splitting of $\\mathbf{A}$, the sequence $\\mathbf{x}^{(k)}$ generated by the above iteration converges to a solution $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$:\n",
"\n",
"> Let $\\rho(\\mathbf{R})=\\max_{i=1,\\dotsc,n}|\\lambda_i(R)|$, where $\\lambda_i(R)$ is the $i$th (complex) eigenvalue of $\\mathbf{R}$. The iteration $\\mathbf{x}^{(t+1)}=\\mathbf{R}\\mathbf{x}^{(t)} - \\mathbf{c}$ converges to a solution to $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$ if and only if $\\rho(\\mathbf{R}) < 1$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Jacobi's method\n",
"\n",
"$$\n",
"x_i^{(t+1)} = \\frac{b_i - \\sum_{j=1}^{i-1} a_{ij} x_j^{(t)} - \\sum_{j=i+1}^n a_{ij} x_j^{(t)}}{a_{ii}}.\n",
"$$\n",
"\n",
"\n",
"* Split $\\mathbf{A} = \\mathbf{L} + \\mathbf{D} + \\mathbf{U}$ = (strictly lower triangular) + (diagonal) + (strictly upper triangular).\n",
"\n",
"\n",
"* Take $\\mathbf{M}=\\mathbf{D}$ (easy to invert!) and $\\mathbf{K}=-(\\mathbf{L} + \\mathbf{U})$:\n",
"$$\n",
" \\mathbf{D} \\mathbf{x}^{(t+1)} = - (\\mathbf{L} + \\mathbf{U}) \\mathbf{x}^{(t)} + \\mathbf{b},\n",
"$$\n",
"i.e., \n",
"$$\n",
"\t\\mathbf{x}^{(t+1)} = - \\mathbf{D}^{-1} (\\mathbf{L} + \\mathbf{U}) \\mathbf{x}^{(t)} + \\mathbf{D}^{-1} \\mathbf{b} = - \\mathbf{D}^{-1} \\mathbf{A} \\mathbf{x}^{(t)} + \\mathbf{x}^{(t)} + \\mathbf{D}^{-1} \\mathbf{b}.\n",
"$$\n",
"\n",
"* Convergence is guaranteed if $\\mathbf{A}$ is striclty row diagonally dominant: $|a_{ii}| > \\sum_{j\\neq i}|a_{ij}|$.\n",
"\n",
"* One round costs $2n^2$ flops with an **unstructured** $\\mathbf{A}$. Gain over GE/LU if converges in $o(n)$ iterations. \n",
"\n",
"* Saving is huge for **sparse** or **structured** $\\mathbf{A}$. By structured, we mean the matrix-vector multiplication $\\mathbf{A} \\mathbf{v}$ is fast ($O(n)$ or $O(n \\log n)$).\n",
" - Often the multiplication is implicit and $\\mathbf{A}$ is not even stored, e.g., finite difference: $(\\mathbf{A}\\mathbf{v})_i = v_i - v_{i+1}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gauss-Seidel method\n",
"\n",
"$$\n",
"x_i^{(t+1)} = \\frac{b_i - \\sum_{j=1}^{i-1} a_{ij} x_j^{(t+1)} - \\sum_{j=i+1}^n a_{ij} x_j^{(t)}}{a_{ii}}.\n",
"$$\n",
"\n",
"* Split $\\mathbf{A} = \\mathbf{L} + \\mathbf{D} + \\mathbf{U}$ = (strictly lower triangular) + (diagonal) + (strictly upper triangular) as Jacobi.\n",
"\n",
"* Take $\\mathbf{M}=\\mathbf{D}+\\mathbf{L}$ (easy to invert, why?) and $\\mathbf{K}=-\\mathbf{U}$:\n",
"$$\n",
"(\\mathbf{D} + \\mathbf{L}) \\mathbf{x}^{(t+1)} = - \\mathbf{U} \\mathbf{x}^{(t)} + \\mathbf{b}\n",
"$$\n",
"i.e., \n",
"$$\n",
"\\mathbf{x}^{(t+1)} = - (\\mathbf{D} + \\mathbf{L})^{-1} \\mathbf{U} \\mathbf{x}^{(t)} + (\\mathbf{D} + \\mathbf{L})^{-1} \\mathbf{b}.\n",
"$$\n",
"\n",
"* Equivalent to\n",
"$$\n",
"\\mathbf{D}\\mathbf{x}^{(t+1)} = - \\mathbf{L} \\mathbf{x}^{(t+1)} - \\mathbf{U} \\mathbf{x}^{(t)} + \\mathbf{b}\n",
"$$\n",
"or\n",
"$$\n",
"\\mathbf{x}^{(t+1)} = \\mathbf{D}^{-1}(- \\mathbf{L} \\mathbf{x}^{(t+1)} - \\mathbf{U} \\mathbf{x}^{(t)} + \\mathbf{b})\n",
"$$\n",
"leading to the iteration.\n",
"\n",
"* \"Coordinate descent\" version of Jacobi.\n",
"\n",
"* Convergence is guaranteed if $\\mathbf{A}$ is striclty row diagonally dominant."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Successive over-relaxation (SOR)\n",
"\n",
"$$\n",
"x_i^{(t+1)} = \\frac{\\omega}{a_{ii}} \\left( b_i - \\sum_{j=1}^{i-1} a_{ij} x_j^{(t+1)} - \\sum_{j=i+1}^n a_{ij} x_j^{(t)} \\right)+ (1-\\omega) x_i^{(t)},\n",
"$$\n",
"\n",
"* $\\omega=1$: Gauss-Seidel; $\\omega \\in (0, 1)$: underrelaxation; $\\omega > 1$: overrelaxation\n",
" - Relaxation in hope of faster convergence\n",
" \n",
"* Split $\\mathbf{A} = \\mathbf{L} + \\mathbf{D} + \\mathbf{U}$ = (strictly lower triangular) + (diagonal) + (strictly upper triangular) as before.\n",
"\n",
"* Take $\\mathbf{M}=\\frac{1}{\\omega}\\mathbf{D}+\\mathbf{L}$ (easy to invert, why?) and $\\mathbf{K}=\\frac{1-\\omega}{\\omega}\\mathbf{D}-\\mathbf{U}$:\n",
"$$\n",
"\\begin{split}\n",
"(\\mathbf{D} + \\omega \\mathbf{L})\\mathbf{x}^{(t+1)} &= [(1-\\omega) \\mathbf{D} - \\omega \\mathbf{U}] \\mathbf{x}^{(t)} + \\omega \\mathbf{b} \n",
"\\\\\n",
"\\mathbf{D}\\mathbf{x}^{(t+1)} &= (1-\\omega) \\mathbf{D}\\mathbf{x}^{(t)} + \\omega ( -\\mathbf{U}\\mathbf{x}^{(t)} - \\mathbf{L}\\mathbf{x}^{(t+1)} + \\mathbf{b} )\n",
"\\\\\n",
"\\mathbf{x}^{(t+1)} &= (1-\\omega)\\mathbf{x}^{(t)} + \\omega \\mathbf{D}^{-1} ( -\\mathbf{U}\\mathbf{x}^{(t)} - \\mathbf{L}\\mathbf{x}^{(t+1)} + \\mathbf{b} )\n",
"\\end{split}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conjugate gradient method\n",
"\n",
"* **Conjugate gradient and its variants are the top-notch iterative methods for solving huge, structured linear systems.**\n",
"\n",
"* Key idea: solving $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$ is equivalent to minimizing the quadratic function $\\frac{1}{2} \\mathbf{x}^T \\mathbf{A} \\mathbf{x} - \\mathbf{b}^T \\mathbf{x}$ if $\\mathbf{A}$ *is positive definite*.\n",
"\n",
"[Application to a fusion problem in physics](http://www.sciencedirect.com/science/article/pii/0021999178900980?via%3Dihub):\n",
"\n",
"| Method | Number of Iterations |\n",
"|----------------------------------------|----------------------|\n",
"| Gauss Seidel | 208,000 |\n",
"| Block SOR methods | 765 |\n",
"| Incomplete Cholesky conjugate gradient | 25 |\n",
"\n",
"\n",
"* More on CG later."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## MatrixDepot.jl\n",
"\n",
"MatrixDepot.jl is an extensive collection of test matrices in Julia. After installation, we can check available test matrices by"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"include group.jl for user defined matrix generators\n",
"verify download of index files...\n",
"reading database\n",
"adding metadata...\n",
"adding svd data...\n",
"writing database\n",
"used remote sites are sparse.tamu.edu with MAT index and math.nist.gov with HTML index\n"
]
},
{
"data": {
"text/latex": [
"\\subsubsection{Currently loaded Matrices}\n",
"\\begin{tabular}\n",
"{l | l | l | l | l}\n",
"builtin(\\#) & & & & \\\\\n",
"\\hline\n",
"1 baart & 13 fiedler & 25 invhilb & 37 parter & 49 shaw \\\\\n",
"2 binomial & 14 forsythe & 26 invol & 38 pascal & 50 smallworld \\\\\n",
"3 blur & 15 foxgood & 27 kahan & 39 pei & 51 spikes \\\\\n",
"4 cauchy & 16 frank & 28 kms & 40 phillips & 52 toeplitz \\\\\n",
"5 chebspec & 17 gilbert & 29 lehmer & 41 poisson & 53 tridiag \\\\\n",
"6 chow & 18 golub & 30 lotkin & 42 prolate & 54 triw \\\\\n",
"7 circul & 19 gravity & 31 magic & 43 randcorr & 55 ursell \\\\\n",
"8 clement & 20 grcar & 32 minij & 44 rando & 56 vand \\\\\n",
"9 companion & 21 hadamard & 33 moler & 45 randsvd & 57 wathen \\\\\n",
"10 deriv2 & 22 hankel & 34 neumann & 46 rohess & 58 wilkinson \\\\\n",
"11 dingdong & 23 heat & 35 oscillate & 47 rosser & 59 wing \\\\\n",
"12 erdrey & 24 hilb & 36 parallax & 48 sampling & \\\\\n",
"\\end{tabular}\n",
"\\begin{tabular}\n",
"{l}\n",
"user(\\#) \\\\\n",
"\\hline\n",
"\\end{tabular}\n",
"\\begin{tabular}\n",
"{l | l | l | l | l | l | l | l | l | l | l | l}\n",
"Groups & & & & & & & & & & & \\\\\n",
"\\hline\n",
"all & local & eigen & illcond & posdef & regprob & symmetric & & & & & \\\\\n",
"builtin & user & graph & inverse & random & sparse & & & & & & \\\\\n",
"\\end{tabular}\n",
"\\begin{tabular}\n",
"{l | l}\n",
"Suite Sparse & of \\\\\n",
"\\hline\n",
"1 & 2893 \\\\\n",
"\\end{tabular}\n",
"\\begin{tabular}\n",
"{l | l}\n",
"MatrixMarket & of \\\\\n",
"\\hline\n",
"0 & 498 \\\\\n",
"\\end{tabular}\n"
],
"text/markdown": [
"### Currently loaded Matrices\n",
"\n",
"| builtin(#) | | | | |\n",
"|:----------- |:----------- |:------------ |:----------- |:------------- |\n",
"| 1 baart | 13 fiedler | 25 invhilb | 37 parter | 49 shaw |\n",
"| 2 binomial | 14 forsythe | 26 invol | 38 pascal | 50 smallworld |\n",
"| 3 blur | 15 foxgood | 27 kahan | 39 pei | 51 spikes |\n",
"| 4 cauchy | 16 frank | 28 kms | 40 phillips | 52 toeplitz |\n",
"| 5 chebspec | 17 gilbert | 29 lehmer | 41 poisson | 53 tridiag |\n",
"| 6 chow | 18 golub | 30 lotkin | 42 prolate | 54 triw |\n",
"| 7 circul | 19 gravity | 31 magic | 43 randcorr | 55 ursell |\n",
"| 8 clement | 20 grcar | 32 minij | 44 rando | 56 vand |\n",
"| 9 companion | 21 hadamard | 33 moler | 45 randsvd | 57 wathen |\n",
"| 10 deriv2 | 22 hankel | 34 neumann | 46 rohess | 58 wilkinson |\n",
"| 11 dingdong | 23 heat | 35 oscillate | 47 rosser | 59 wing |\n",
"| 12 erdrey | 24 hilb | 36 parallax | 48 sampling | |\n",
"\n",
"| user(#) |\n",
"|:------- |\n",
"\n",
"| Groups | | | | | | | | | | | |\n",
"|:------- |:----- |:----- |:------- |:------ |:------- |:--------- |:--- |:--- |:--- |:--- |:--- |\n",
"| all | local | eigen | illcond | posdef | regprob | symmetric | | | | | |\n",
"| builtin | user | graph | inverse | random | sparse | | | | | | |\n",
"\n",
"| Suite Sparse | of |\n",
"|:------------ |:---- |\n",
"| 1 | 2893 |\n",
"\n",
"| MatrixMarket | of |\n",
"|:------------ |:--- |\n",
"| 0 | 498 |\n"
],
"text/plain": [
"\u001b[1m Currently loaded Matrices\u001b[22m\n",
"\u001b[1m –––––––––––––––––––––––––––\u001b[22m\n",
"\n",
" builtin(#) \n",
" ––––––––––– ––––––––––– –––––––––––– ––––––––––– –––––––––––––\n",
" 1 baart 13 fiedler 25 invhilb 37 parter 49 shaw \n",
" 2 binomial 14 forsythe 26 invol 38 pascal 50 smallworld\n",
" 3 blur 15 foxgood 27 kahan 39 pei 51 spikes \n",
" 4 cauchy 16 frank 28 kms 40 phillips 52 toeplitz \n",
" 5 chebspec 17 gilbert 29 lehmer 41 poisson 53 tridiag \n",
" 6 chow 18 golub 30 lotkin 42 prolate 54 triw \n",
" 7 circul 19 gravity 31 magic 43 randcorr 55 ursell \n",
" 8 clement 20 grcar 32 minij 44 rando 56 vand \n",
" 9 companion 21 hadamard 33 moler 45 randsvd 57 wathen \n",
" 10 deriv2 22 hankel 34 neumann 46 rohess 58 wilkinson \n",
" 11 dingdong 23 heat 35 oscillate 47 rosser 59 wing \n",
" 12 erdrey 24 hilb 36 parallax 48 sampling \n",
"\n",
" user(#)\n",
" –––––––\n",
"\n",
" Groups \n",
" ––––––– ––––– ––––– ––––––– –––––– ––––––– ––––––––– \n",
" all local eigen illcond posdef regprob symmetric \n",
" builtin user graph inverse random sparse \n",
"\n",
" Suite Sparse of \n",
" –––––––––––– ––––\n",
" 1 2893\n",
"\n",
" MatrixMarket of \n",
" –––––––––––– –––\n",
" 0 498"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using MatrixDepot\n",
"\n",
"mdinfo()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Vector{String}:\n",
" \"poisson\"\n",
" \"wathen\""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# List matrices that are positive definite and sparse:\n",
"mdlist(:symmetric & :posdef & :sparse)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\section{Poisson Matrix (poisson)}\n",
"A block tridiagonal matrix from Poisson’s equation. This matrix is sparse, symmetric positive definite and has known eigenvalues. The result is of type \\texttt{SparseMatrixCSC}.\n",
"\n",
"\\emph{Input options:}\n",
"\n",
"\\begin{itemize}\n",
"\\item [type,] dim: the dimension of the matirx is \\texttt{dim\\^{}2}.\n",
"\n",
"\\end{itemize}\n",
"\\emph{Groups:} [\"inverse\", \"symmetric\", \"posdef\", \"eigen\", \"sparse\"]\n",
"\n",
"\\emph{References:}\n",
"\n",
"\\textbf{G. H. Golub and C. F. Van Loan}, Matrix Computations, second edition, Johns Hopkins University Press, Baltimore, Maryland, 1989 (Section 4.5.4).\n",
"\n"
],
"text/markdown": [
"# Poisson Matrix (poisson)\n",
"\n",
"A block tridiagonal matrix from Poisson’s equation. This matrix is sparse, symmetric positive definite and has known eigenvalues. The result is of type `SparseMatrixCSC`.\n",
"\n",
"*Input options:*\n",
"\n",
" * [type,] dim: the dimension of the matirx is `dim^2`.\n",
"\n",
"*Groups:* [\"inverse\", \"symmetric\", \"posdef\", \"eigen\", \"sparse\"]\n",
"\n",
"*References:*\n",
"\n",
"**G. H. Golub and C. F. Van Loan**, Matrix Computations, second edition, Johns Hopkins University Press, Baltimore, Maryland, 1989 (Section 4.5.4).\n"
],
"text/plain": [
"\u001b[1m Poisson Matrix (poisson)\u001b[22m\n",
"\u001b[1m ≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡\u001b[22m\n",
"\n",
" A block tridiagonal matrix from Poisson’s equation. This matrix is sparse,\n",
" symmetric positive definite and has known eigenvalues. The result is of type\n",
" \u001b[36mSparseMatrixCSC\u001b[39m.\n",
"\n",
" \u001b[4mInput options:\u001b[24m\n",
"\n",
" • [type,] dim: the dimension of the matirx is \u001b[36mdim^2\u001b[39m.\n",
"\n",
" \u001b[4mGroups:\u001b[24m [\"inverse\", \"symmetric\", \"posdef\", \"eigen\", \"sparse\"]\n",
"\n",
" \u001b[4mReferences:\u001b[24m\n",
"\n",
" \u001b[1mG. H. Golub and C. F. Van Loan\u001b[22m, Matrix Computations, second edition, Johns\n",
" Hopkins University Press, Baltimore, Maryland, 1989 (Section 4.5.4)."
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Get help on Poisson matrix\n",
"mdinfo(\"poisson\")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"100×100 SparseArrays.SparseMatrixCSC{Float64, Int64} with 460 stored entries:\n",
"⠻⣦⡀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠈⠻⣦⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠑⢄⠀⠀⠻⣦⡀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠑⢄⠀⠈⠻⣦⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠑⢄⠀⠀⠻⣦⡀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠑⢄⠀⠈⠻⣦⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠻⣦⡀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠈⠻⣦⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠻⣦⡀⠀⠱⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠈⠻⣦⠀⠀⠱⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢆⠀⠀⠻⣦⡀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢆⠀⠈⠻⣦⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠻⣦⡀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠈⠻⣦⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠻⣦⡀⠀⠑⢄⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠈⠻⣦⠀⠀⠑⢄⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠻⣦⡀⠀⠑⢄⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠈⠻⣦⠀⠀⠑⢄\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠻⣦⡀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠈⠻⣦"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Generate a Poisson matrix of dimension n=10\n",
"A = matrixdepot(\"poisson\", 10)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\u001b[1mSparsity Pattern\u001b[22m⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ \n",
" \u001b[90m┌──────────────────────────────────────────┐\u001b[39m \n",
" \u001b[90m1\u001b[39m \u001b[90m│\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠘\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \u001b[31m> 0\u001b[39m\n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠉\u001b[39m\u001b[38;5;4m⢆\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \u001b[34m< 0\u001b[39m\n",
" \u001b[90m│\u001b[39m\u001b[38;5;4m⠒\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠱\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⢢\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠣\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⠤\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠱\u001b[39m\u001b[38;5;4m⣀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠱\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⡄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠘\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠱\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠉\u001b[39m\u001b[38;5;4m⢆\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;4m⠒\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⢢\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠣\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⠤\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠱\u001b[39m\u001b[38;5;4m⣀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⡄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠘\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠉\u001b[39m\u001b[38;5;4m⢆\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;4m⠒\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⢢\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠣\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⠤\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠱\u001b[39m\u001b[38;5;4m⣀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⢆\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⡄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠘\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⢆\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠉\u001b[39m\u001b[38;5;4m⢆\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;4m⠒\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⢢\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠣\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⢆\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⠤\u001b[39m\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠱\u001b[39m\u001b[38;5;4m⣀\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m100\u001b[39m \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;4m⡄\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[90m│\u001b[39m \n",
" \u001b[90m└──────────────────────────────────────────┘\u001b[39m \n",
" ⠀\u001b[90m1\u001b[39m⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\u001b[90m100\u001b[39m⠀ \n",
" ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\u001b[0mnz = 460⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ "
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using UnicodePlots\n",
"spy(A)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\section{Wathen Matrix (wathen)}\n",
"Wathen Matrix is a sparse, symmetric positive, random matrix arose from the finite element method. The generated matrix is the consistent mass matrix for a regular nx-by-ny grid of 8-nodes.\n",
"\n",
"\\emph{Input options:}\n",
"\n",
"\\begin{itemize}\n",
"\\item [type,] nx, ny: the dimension of the matrix is equal to \\texttt{3 * nx * ny + 2 * nx * ny + 1}.\n",
"\n",
"\n",
"\\item [type,] n: \\texttt{nx = ny = n}.\n",
"\n",
"\\end{itemize}\n",
"\\emph{Groups:} [\"symmetric\", \"posdef\", \"eigen\", \"random\", \"sparse\"]\n",
"\n",
"\\emph{References:}\n",
"\n",
"\\textbf{A. J. Wathen}, Realistic eigenvalue bounds for the Galerkin mass matrix, IMA J. Numer. Anal., 7 (1987), pp. 449-457.\n",
"\n"
],
"text/markdown": [
"# Wathen Matrix (wathen)\n",
"\n",
"Wathen Matrix is a sparse, symmetric positive, random matrix arose from the finite element method. The generated matrix is the consistent mass matrix for a regular nx-by-ny grid of 8-nodes.\n",
"\n",
"*Input options:*\n",
"\n",
" * [type,] nx, ny: the dimension of the matrix is equal to `3 * nx * ny + 2 * nx * ny + 1`.\n",
" * [type,] n: `nx = ny = n`.\n",
"\n",
"*Groups:* [\"symmetric\", \"posdef\", \"eigen\", \"random\", \"sparse\"]\n",
"\n",
"*References:*\n",
"\n",
"**A. J. Wathen**, Realistic eigenvalue bounds for the Galerkin mass matrix, IMA J. Numer. Anal., 7 (1987), pp. 449-457.\n"
],
"text/plain": [
"\u001b[1m Wathen Matrix (wathen)\u001b[22m\n",
"\u001b[1m ≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡\u001b[22m\n",
"\n",
" Wathen Matrix is a sparse, symmetric positive, random matrix arose from the\n",
" finite element method. The generated matrix is the consistent mass matrix\n",
" for a regular nx-by-ny grid of 8-nodes.\n",
"\n",
" \u001b[4mInput options:\u001b[24m\n",
"\n",
" • [type,] nx, ny: the dimension of the matrix is equal to \u001b[36m3 * nx *\n",
" ny + 2 * nx * ny + 1\u001b[39m.\n",
"\n",
" • [type,] n: \u001b[36mnx = ny = n\u001b[39m.\n",
"\n",
" \u001b[4mGroups:\u001b[24m [\"symmetric\", \"posdef\", \"eigen\", \"random\", \"sparse\"]\n",
"\n",
" \u001b[4mReferences:\u001b[24m\n",
"\n",
" \u001b[1mA. J. Wathen\u001b[22m, Realistic eigenvalue bounds for the Galerkin mass matrix, IMA\n",
" J. Numer. Anal., 7 (1987), pp. 449-457."
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Get help on Wathen matrix\n",
"mdinfo(\"wathen\")"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"96×96 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1256 stored entries:\n",
"⠿⣧⣀⠀⠸⣧⡀⠿⣧⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠘⢻⣶⡀⠘⣇⠀⠘⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠶⣦⣀⠈⠱⣦⡉⠶⣦⣀⠈⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⣤⡌⠉⠙⢣⡌⠛⣤⡌⠉⠙⢣⡄⠀⣤⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠉⢻⣶⣀⠈⢻⡆⠉⢻⣶⣀⠈⢻⣆⠉⠛⣷⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠘⠻⠆⠀⠷⣀⡀⠘⠻⢆⡀⠻⣀⡀⠘⠻⠶⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠉⠻⢶⣤⡈⠻⣦⠉⠛⢷⣤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠿⣧⠀⠀⠸⣧⠀⠿⣧⡄⠀⠸⣧⠀⠿⣧⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠙⢻⣶⡀⠙⣷⠀⠉⢻⣶⣀⠙⣷⡀⠉⢻⣶⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠃⠀⠘⠶⣦⣄⠘⠻⣦⡘⠷⣦⣄⠘⠛⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣤⡄⠙⠻⢶⡌⠻⣦⡄⠙⠻⠶⡄⠀⢠⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠿⣧⣀⠈⢿⣄⠉⠿⣧⣀⠀⢿⣄⠈⠿⣧⣄⠀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢻⣶⠀⢻⡆⠀⠘⢻⣶⠀⢻⡆⠀⠀⢻⣶⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠛⢷⣤⣀⠻⣦⡈⠛⠷⣦⣀⠀⠀⠀⠀⠀⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠶⣦⡄⠈⠉⣦⠈⠱⣦⡄⠈⠉⢶⠀⠰⣦⡄⠀⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢿⣤⣀⠹⣧⡀⠉⠿⣧⣀⠸⣧⡀⠉⠿⣧⣀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠛⠀⠘⢣⣄⣀⡘⠛⣤⡘⢣⣄⣀⡘⠛\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡀⠉⠻⠶⣈⠻⢆⡀⠉⠻⠶\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⡄⠀⢹⡄⠈⠿⣧⡄⠀\n",
"⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢻⣶⠈⢻⡆⠀⠉⢻⣶"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Generate a Wathen matrix of dimension n=5\n",
"A = matrixdepot(\"wathen\", 5)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\u001b[1mSparsity Pattern\u001b[22m⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ \n",
" \u001b[90m┌──────────────────────────────────────────┐\u001b[39m \n",
" \u001b[90m1\u001b[39m \u001b[90m│\u001b[39m\u001b[38;5;5m⠿\u001b[39m\u001b[38;5;5m⣧\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⢿\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[38;5;5m⠸\u001b[39m\u001b[38;5;5m⢿\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \u001b[31m> 0\u001b[39m\n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⠿\u001b[39m\u001b[38;5;5m⣧\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⣧\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[38;5;1m⠈\u001b[39m\u001b[38;5;5m⠹\u001b[39m\u001b[38;5;5m⣧\u001b[39m\u001b[38;5;5m⣄\u001b[39m\u001b[38;5;1m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \u001b[34m< 0\u001b[39m\n",
" \u001b[90m│\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;5m⣄\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;1m⣤\u001b[39m\u001b[38;5;5m⡘\u001b[39m\u001b[38;5;5m⢣\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⠃\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[38;5;5m⡉\u001b[39m\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⠶\u001b[39m\u001b[38;5;5m⣈\u001b[39m\u001b[38;5;1m⠻\u001b[39m\u001b[38;5;1m⢆\u001b[39m\u001b[38;5;5m⣈\u001b[39m\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⠷\u001b[39m\u001b[38;5;4m⢆\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⣀\u001b[39m\u001b[38;5;1m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⣷\u001b[39m\u001b[38;5;5m⣆\u001b[39m\u001b[38;5;5m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;5m⢻\u001b[39m\u001b[38;5;5m⡆\u001b[39m\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⢻\u001b[39m\u001b[38;5;5m⣶\u001b[39m\u001b[38;5;5m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;5m⢸\u001b[39m\u001b[38;5;5m⣆\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⣷\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⢿\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[0m⠀\u001b[38;5;5m⢿\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠈\u001b[39m\u001b[38;5;5m⠿\u001b[39m\u001b[38;5;5m⣧\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[38;5;5m⠹\u001b[39m\u001b[38;5;5m⣧\u001b[39m\u001b[0m⠀\u001b[38;5;1m⠈\u001b[39m\u001b[38;5;5m⠹\u001b[39m\u001b[38;5;5m⢿\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;1m⠈\u001b[39m\u001b[38;5;5m⠉\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;4m⠑\u001b[39m\u001b[38;5;5m⠲\u001b[39m\u001b[38;5;5m⢶\u001b[39m\u001b[38;5;5m⣄\u001b[39m\u001b[38;5;5m⡉\u001b[39m\u001b[38;5;1m⠱\u001b[39m\u001b[38;5;1m⣦\u001b[39m\u001b[38;5;5m⡉\u001b[39m\u001b[38;5;5m⠒\u001b[39m\u001b[38;5;5m⢶\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;5m⣈\u001b[39m\u001b[38;5;5m⠁\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⢠\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⢣\u001b[39m\u001b[38;5;1m⠈\u001b[39m\u001b[38;5;1m⠛\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠙\u001b[39m\u001b[38;5;5m⢣\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[0m⠀\u001b[38;5;5m⢠\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;1m⠈\u001b[39m\u001b[38;5;5m⠙\u001b[39m\u001b[38;5;5m⢻\u001b[39m\u001b[38;5;5m⣆\u001b[39m\u001b[38;5;1m⡀\u001b[39m\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⣷\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⢻\u001b[39m\u001b[38;5;5m⣶\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⢻\u001b[39m\u001b[38;5;5m⣆\u001b[39m\u001b[38;5;5m⠈\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⣷\u001b[39m\u001b[38;5;5m⣆\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⠷\u001b[39m\u001b[38;5;5m⠆\u001b[39m\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⠷\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⢆\u001b[39m\u001b[38;5;1m⡀\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⢆\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⠶\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⢶\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;5m⡈\u001b[39m\u001b[38;5;1m⠻\u001b[39m\u001b[38;5;1m⣦\u001b[39m\u001b[38;5;5m⡈\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⠷\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⠶\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠱\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;1m⠈\u001b[39m\u001b[38;5;5m⠱\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⢶\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[38;5;5m⠰\u001b[39m\u001b[38;5;5m⢶\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⠹\u001b[39m\u001b[38;5;5m⢿\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;5m⡀\u001b[39m\u001b[38;5;5m⠹\u001b[39m\u001b[38;5;5m⣧\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⠿\u001b[39m\u001b[38;5;5m⣧\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⢿\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[38;5;1m⠈\u001b[39m\u001b[38;5;5m⠹\u001b[39m\u001b[38;5;5m⣧\u001b[39m\u001b[38;5;5m⣄\u001b[39m\u001b[38;5;1m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⠃\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⢣\u001b[39m\u001b[38;5;5m⣄\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;1m⣤\u001b[39m\u001b[38;5;1m⡀\u001b[39m\u001b[38;5;5m⢣\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⠃\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⢀\u001b[39m\u001b[38;5;5m⡉\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⠷\u001b[39m\u001b[38;5;5m⠤\u001b[39m\u001b[38;5;5m⣈\u001b[39m\u001b[38;5;1m⠻\u001b[39m\u001b[38;5;1m⢆\u001b[39m\u001b[38;5;5m⣈\u001b[39m\u001b[38;5;5m⠙\u001b[39m\u001b[38;5;5m⠷\u001b[39m\u001b[38;5;5m⠦\u001b[39m\u001b[38;5;4m⢄\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;5m⣀\u001b[39m\u001b[38;5;1m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⣷\u001b[39m\u001b[38;5;5m⣆\u001b[39m\u001b[38;5;1m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;5m⢻\u001b[39m\u001b[38;5;5m⣆\u001b[39m\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⢻\u001b[39m\u001b[38;5;5m⣶\u001b[39m\u001b[38;5;5m⡀\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⣷\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⣷\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⢿\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠹\u001b[39m\u001b[38;5;5m⡇\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠈\u001b[39m\u001b[38;5;5m⠿\u001b[39m\u001b[38;5;5m⣧\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[38;5;5m⠸\u001b[39m\u001b[38;5;5m⣧\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠈\u001b[39m\u001b[38;5;5m⠹\u001b[39m\u001b[38;5;5m⢿\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;1m⠈\u001b[39m\u001b[38;5;5m⠉\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠱\u001b[39m\u001b[38;5;5m⢶\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[38;5;5m⡉\u001b[39m\u001b[38;5;1m⠱\u001b[39m\u001b[38;5;1m⣦\u001b[39m\u001b[38;5;5m⡉\u001b[39m\u001b[38;5;5m⠶\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[38;5;5m⣈\u001b[39m\u001b[38;5;5m⠉\u001b[39m\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⢠\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[0m⠀\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⢣\u001b[39m\u001b[38;5;5m⡌\u001b[39m\u001b[38;5;1m⠛\u001b[39m\u001b[38;5;5m⣤\u001b[39m\u001b[38;5;5m⡄\u001b[39m\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠙\u001b[39m\u001b[38;5;5m⠛\u001b[39m\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;1m⠈\u001b[39m\u001b[38;5;5m⠙\u001b[39m\u001b[38;5;5m⢻\u001b[39m\u001b[38;5;5m⣆\u001b[39m\u001b[38;5;1m⡀\u001b[39m\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⢻\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[38;5;5m⠉\u001b[39m\u001b[38;5;5m⢻\u001b[39m\u001b[38;5;5m⣶\u001b[39m\u001b[38;5;5m⣀\u001b[39m\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m96\u001b[39m \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⠛\u001b[39m\u001b[38;5;5m⣷\u001b[39m\u001b[38;5;5m⡆\u001b[39m\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⣷\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[38;5;5m⠘\u001b[39m\u001b[38;5;5m⢻\u001b[39m\u001b[38;5;5m⣶\u001b[39m\u001b[90m│\u001b[39m \n",
" \u001b[90m└──────────────────────────────────────────┘\u001b[39m \n",
" ⠀\u001b[90m1\u001b[39m⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\u001b[90m96\u001b[39m⠀ \n",
" ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\u001b[0mnz = 1256⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ "
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"spy(A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Numerical examples\n",
"\n",
"The [`IterativeSolvers.jl`](https://github.com/JuliaMath/IterativeSolvers.jl) package implements most commonly used iterative solvers."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Generate test matrix"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(A) = SparseArrays.SparseMatrixCSC{Float64, Int64}\n",
"typeof(Afull) = Matrix{Float64}\n"
]
},
{
"data": {
"text/plain": [
"0.000496"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using BenchmarkTools, IterativeSolvers, LinearAlgebra, MatrixDepot, Random\n",
"\n",
"Random.seed!(280)\n",
"\n",
"n = 100\n",
"# Poisson matrix of dimension n^2=10000, pd and sparse\n",
"A = matrixdepot(\"poisson\", n)\n",
"@show typeof(A)\n",
"# dense matrix representation of A\n",
"Afull = convert(Matrix, A)\n",
"@show typeof(Afull)\n",
"# sparsity level\n",
"count(!iszero, A) / length(A)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\u001b[1mSparsity Pattern\u001b[22m⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ \n",
" \u001b[90m┌──────────────────────────────────────────┐\u001b[39m \n",
" \u001b[90m1\u001b[39m \u001b[90m│\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \u001b[31m> 0\u001b[39m\n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \u001b[34m< 0\u001b[39m\n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[38;5;4m⡀\u001b[39m\u001b[0m⠀\u001b[90m│\u001b[39m \n",
" \u001b[90m10000\u001b[39m \u001b[90m│\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[38;5;4m⠈\u001b[39m\u001b[38;5;5m⠻\u001b[39m\u001b[38;5;5m⣦\u001b[39m\u001b[90m│\u001b[39m \n",
" \u001b[90m└──────────────────────────────────────────┘\u001b[39m \n",
" ⠀\u001b[90m1\u001b[39m⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\u001b[90m10000\u001b[39m⠀ \n",
" ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀\u001b[0mnz = 49600⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ "
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"spy(A)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(873768, 800000040)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# storage difference\n",
"Base.summarysize(A), Base.summarysize(Afull)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Matrix-vector muliplication"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"BenchmarkTools.Trial: 116 samples with 1 evaluation.\n",
" Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m): \u001b[39m\u001b[36m\u001b[1m41.223 ms\u001b[22m\u001b[39m … \u001b[35m 45.530 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.00% … 0.00%\n",
" Time \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m): \u001b[39m\u001b[34m\u001b[1m43.316 ms \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m): \u001b[39m0.00%\n",
" Time \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m): \u001b[39m\u001b[32m\u001b[1m43.407 ms\u001b[22m\u001b[39m ± \u001b[32m846.034 μs\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m): \u001b[39m0.00% ± 0.00%\n",
"\n",
" \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▂\u001b[39m▄\u001b[39m \u001b[39m▂\u001b[39m▄\u001b[39m█\u001b[39m▄\u001b[39m▂\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[34m▄\u001b[39m\u001b[39m \u001b[32m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▆\u001b[39m█\u001b[39m▄\u001b[39m▂\u001b[39m▆\u001b[39m \u001b[39m \u001b[39m▂\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▂\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \n",
" \u001b[39m▄\u001b[39m▁\u001b[39m▄\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▄\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▄\u001b[39m▆\u001b[39m▆\u001b[39m▄\u001b[34m█\u001b[39m\u001b[39m▄\u001b[32m▄\u001b[39m\u001b[39m▄\u001b[39m▄\u001b[39m▁\u001b[39m▆\u001b[39m▆\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▄\u001b[39m█\u001b[39m▁\u001b[39m▄\u001b[39m▆\u001b[39m▄\u001b[39m█\u001b[39m█\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▄\u001b[39m \u001b[39m▄\n",
" 41.2 ms\u001b[90m Histogram: frequency by time\u001b[39m 45.3 ms \u001b[0m\u001b[1m<\u001b[22m\n",
"\n",
" Memory estimate\u001b[90m: \u001b[39m\u001b[33m78.20 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m2\u001b[39m."
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# randomly generated vector of length n^2\n",
"b = randn(n^2)\n",
"# dense matrix-vector multiplication\n",
"@benchmark $Afull * $b"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"BenchmarkTools.Trial: 10000 samples with 1 evaluation.\n",
" Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m): \u001b[39m\u001b[36m\u001b[1m63.610 μs\u001b[22m\u001b[39m … \u001b[35m 44.016 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.00% … 99.70%\n",
" Time \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m): \u001b[39m\u001b[34m\u001b[1m85.806 μs \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m): \u001b[39m0.00%\n",
" Time \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m): \u001b[39m\u001b[32m\u001b[1m94.585 μs\u001b[22m\u001b[39m ± \u001b[32m439.588 μs\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m): \u001b[39m4.64% ± 1.00%\n",
"\n",
" \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▄\u001b[39m▇\u001b[39m█\u001b[39m█\u001b[34m▇\u001b[39m\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[32m▄\u001b[39m\u001b[39m▃\u001b[39m▂\u001b[39m▂\u001b[39m \u001b[39m▁\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▂\n",
" \u001b[39m▆\u001b[39m▅\u001b[39m▆\u001b[39m▄\u001b[39m▃\u001b[39m▁\u001b[39m▄\u001b[39m▅\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[34m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[32m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m▆\u001b[39m▆\u001b[39m▇\u001b[39m▇\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m█\u001b[39m▇\u001b[39m▅\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▆\u001b[39m▇\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▆\u001b[39m▇\u001b[39m▇\u001b[39m▆\u001b[39m▇\u001b[39m \u001b[39m█\n",
" 63.6 μs\u001b[90m \u001b[39m\u001b[90mHistogram: \u001b[39m\u001b[90m\u001b[1mlog(\u001b[22m\u001b[39m\u001b[90mfrequency\u001b[39m\u001b[90m\u001b[1m)\u001b[22m\u001b[39m\u001b[90m by time\u001b[39m 179 μs \u001b[0m\u001b[1m<\u001b[22m\n",
"\n",
" Memory estimate\u001b[90m: \u001b[39m\u001b[33m78.20 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m2\u001b[39m."
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# sparse matrix-vector multiplication\n",
"@benchmark $A * $b"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Dense solve via Cholesky"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"# record the Cholesky solution\n",
"xchol = cholesky(Afull) \\ b;"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"BenchmarkTools.Trial: 2 samples with 1 evaluation.\n",
" Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m): \u001b[39m\u001b[36m\u001b[1m2.933 s\u001b[22m\u001b[39m … \u001b[35m 3.085 s\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.02% … 1.46%\n",
" Time \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m): \u001b[39m\u001b[34m\u001b[1m3.009 s \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m): \u001b[39m0.76%\n",
" Time \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m): \u001b[39m\u001b[32m\u001b[1m3.009 s\u001b[22m\u001b[39m ± \u001b[32m107.518 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m): \u001b[39m0.76% ± 1.02%\n",
"\n",
" \u001b[34m█\u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[32m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m \u001b[39m \n",
" \u001b[34m█\u001b[39m\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[32m▁\u001b[39m\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m \u001b[39m▁\n",
" 2.93 s\u001b[90m Histogram: frequency by time\u001b[39m 3.09 s \u001b[0m\u001b[1m<\u001b[22m\n",
"\n",
" Memory estimate\u001b[90m: \u001b[39m\u001b[33m763.02 MiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m6\u001b[39m."
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# dense solve via Cholesky\n",
"@benchmark cholesky($Afull) \\ $b"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Jacobi solver\n",
"\n",
"It seems that Jacobi solver doesn't give the correct answer."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm(xjacobi - xchol) = 553.3696146722076\n"
]
},
{
"data": {
"text/plain": [
"553.3696146722076"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xjacobi = jacobi(A, b)\n",
"@show norm(xjacobi - xchol)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Documentation](https://juliamath.github.io/IterativeSolvers.jl/dev/linear_systems/stationary/#Jacobi-1) reveals that the default value of `maxiter` is 10. A couple trial runs shows that 30000 Jacobi iterations are enough to get an accurate solution."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm(xjacobi - xchol) = 0.0001726309300532125\n"
]
},
{
"data": {
"text/plain": [
"0.0001726309300532125"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xjacobi = jacobi(A, b, maxiter=30000)\n",
"@show norm(xjacobi - xchol)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"BenchmarkTools.Trial: 3 samples with 1 evaluation.\n",
" Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m): \u001b[39m\u001b[36m\u001b[1m2.087 s\u001b[22m\u001b[39m … \u001b[35m 2.209 s\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.00% … 0.00%\n",
" Time \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m): \u001b[39m\u001b[34m\u001b[1m2.124 s \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m): \u001b[39m0.00%\n",
" Time \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m): \u001b[39m\u001b[32m\u001b[1m2.140 s\u001b[22m\u001b[39m ± \u001b[32m62.377 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m): \u001b[39m0.00% ± 0.00%\n",
"\n",
" \u001b[34m█\u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[32m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m \u001b[39m \n",
" \u001b[34m█\u001b[39m\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[32m▁\u001b[39m\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m \u001b[39m▁\n",
" 2.09 s\u001b[90m Histogram: frequency by time\u001b[39m 2.21 s \u001b[0m\u001b[1m<\u001b[22m\n",
"\n",
" Memory estimate\u001b[90m: \u001b[39m\u001b[33m2.06 MiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m30008\u001b[39m."
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@benchmark jacobi($A, $b, maxiter=30000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Gauss-Seidel"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm(xgs - xchol) = 0.00016756826216903108\n"
]
},
{
"data": {
"text/plain": [
"0.00016756826216903108"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Gauss-Seidel solution is fairly close to Cholesky solution after 15000 iters\n",
"xgs = gauss_seidel(A, b, maxiter=15000)\n",
"@show norm(xgs - xchol)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"BenchmarkTools.Trial: 4 samples with 1 evaluation.\n",
" Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m): \u001b[39m\u001b[36m\u001b[1m1.596 s\u001b[22m\u001b[39m … \u001b[35m 1.610 s\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.00% … 0.00%\n",
" Time \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m): \u001b[39m\u001b[34m\u001b[1m1.598 s \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m): \u001b[39m0.00%\n",
" Time \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m): \u001b[39m\u001b[32m\u001b[1m1.601 s\u001b[22m\u001b[39m ± \u001b[32m6.595 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m): \u001b[39m0.00% ± 0.00%\n",
"\n",
" \u001b[39m█\u001b[39m \u001b[39m \u001b[39m█\u001b[34m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[32m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m \u001b[39m \n",
" \u001b[39m█\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[34m▁\u001b[39m\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[32m▁\u001b[39m\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m \u001b[39m▁\n",
" 1.6 s\u001b[90m Histogram: frequency by time\u001b[39m 1.61 s \u001b[0m\u001b[1m<\u001b[22m\n",
"\n",
" Memory estimate\u001b[90m: \u001b[39m\u001b[33m156.64 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m7\u001b[39m."
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@benchmark gauss_seidel($A, $b, maxiter=15000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### SOR"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm(xsor - xchol) = 0.003162012427719623\n"
]
},
{
"data": {
"text/plain": [
"0.003162012427719623"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# symmetric SOR with ω=0.75\n",
"xsor = ssor(A, b, 0.75, maxiter=10000)\n",
"@show norm(xsor - xchol)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"BenchmarkTools.Trial: 4 samples with 1 evaluation.\n",
" Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m): \u001b[39m\u001b[36m\u001b[1m1.262 s\u001b[22m\u001b[39m … \u001b[35m 1.283 s\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.00% … 0.00%\n",
" Time \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m): \u001b[39m\u001b[34m\u001b[1m1.271 s \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m): \u001b[39m0.00%\n",
" Time \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m): \u001b[39m\u001b[32m\u001b[1m1.272 s\u001b[22m\u001b[39m ± \u001b[32m8.957 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m): \u001b[39m0.00% ± 0.00%\n",
"\n",
" \u001b[39m█\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[34m█\u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[32m \u001b[39m\u001b[39m█\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m \u001b[39m \n",
" \u001b[39m█\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[34m█\u001b[39m\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[32m▁\u001b[39m\u001b[39m█\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m \u001b[39m▁\n",
" 1.26 s\u001b[90m Histogram: frequency by time\u001b[39m 1.28 s \u001b[0m\u001b[1m<\u001b[22m\n",
"\n",
" Memory estimate\u001b[90m: \u001b[39m\u001b[33m547.36 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m20009\u001b[39m."
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@benchmark sor($A, $b, 0.75, maxiter=10000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Conjugate Gradient"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm(xcg - xchol) = 1.4922213222422449e-5\n"
]
},
{
"data": {
"text/plain": [
"1.4922213222422449e-5"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# conjugate gradient\n",
"xcg = cg(A, b)\n",
"@show norm(xcg - xchol)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"BenchmarkTools.Trial: 237 samples with 1 evaluation.\n",
" Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m): \u001b[39m\u001b[36m\u001b[1m20.092 ms\u001b[22m\u001b[39m … \u001b[35m 25.142 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.00% … 0.00%\n",
" Time \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m): \u001b[39m\u001b[34m\u001b[1m21.049 ms \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m): \u001b[39m0.00%\n",
" Time \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m): \u001b[39m\u001b[32m\u001b[1m21.167 ms\u001b[22m\u001b[39m ± \u001b[32m678.097 μs\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m): \u001b[39m0.00% ± 0.00%\n",
"\n",
" \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m▁\u001b[39m \u001b[39m▁\u001b[39m▂\u001b[39m▄\u001b[39m▄\u001b[34m▄\u001b[39m\u001b[39m▄\u001b[32m▁\u001b[39m\u001b[39m \u001b[39m \u001b[39m▅\u001b[39m▁\u001b[39m▂\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \n",
" \u001b[39m▅\u001b[39m▆\u001b[39m▄\u001b[39m▅\u001b[39m▅\u001b[39m▆\u001b[39m▆\u001b[39m▇\u001b[39m█\u001b[39m▇\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[34m█\u001b[39m\u001b[39m█\u001b[32m█\u001b[39m\u001b[39m▆\u001b[39m▆\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▆\u001b[39m▄\u001b[39m▅\u001b[39m▄\u001b[39m▇\u001b[39m▅\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▅\u001b[39m▄\u001b[39m▄\u001b[39m▁\u001b[39m▅\u001b[39m▃\u001b[39m▃\u001b[39m▃\u001b[39m▁\u001b[39m▄\u001b[39m▁\u001b[39m▁\u001b[39m▃\u001b[39m▃\u001b[39m▃\u001b[39m▃\u001b[39m▁\u001b[39m▁\u001b[39m▃\u001b[39m▄\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▃\u001b[39m \u001b[39m▄\n",
" 20.1 ms\u001b[90m Histogram: frequency by time\u001b[39m 23.4 ms \u001b[0m\u001b[1m<\u001b[22m\n",
"\n",
" Memory estimate\u001b[90m: \u001b[39m\u001b[33m313.61 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m16\u001b[39m."
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@benchmark cg($A, $b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* A basic tenet in numerical analysis: \n",
"\n",
"> **The structure should be exploited whenever possible in solving a problem.** "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Acknowledgment\n",
"\n",
"Many parts of this lecture note is based on [Dr. Hua Zhou](http://hua-zhou.github.io)'s 2019 Spring Statistical Computing course notes available at ."
]
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"kernelspec": {
"display_name": "Julia 1.6.2",
"language": "julia",
"name": "julia-1.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.6.2"
},
"toc": {
"colors": {
"hover_highlight": "#DAA520",
"running_highlight": "#FF0000",
"selected_highlight": "#FFD700"
},
"moveMenuLeft": true,
"nav_menu": {
"height": "143px",
"width": "252px"
},
"navigate_menu": true,
"number_sections": false,
"sideBar": true,
"skip_h1_title": true,
"threshold": 4,
"toc_cell": true,
"toc_position": {
"height": "443.5px",
"left": "0px",
"right": "796px",
"top": "67px",
"width": "164px"
},
"toc_section_display": "block",
"toc_window_display": true,
"widenNotebook": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}