{ "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 }