{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "

1  Iterative Methods for Solving Linear Equations
1.1  Introduction
1.2  PageRank problem
1.3  Jacobi method
1.4  Gauss-Seidel method
1.5  Successive over-relaxation (SOR)
1.6  Conjugate gradient method
1.7  MatrixDepot.jl
1.7.1  Currently loaded Matrices
2  Poisson Matrix (poisson)
3  Wathen Matrix (wathen)
3.1  Numerical examples
3.1.1  Generate test matrix
3.1.2  Matrix-vector muliplication
3.1.3  Dense solve via Cholesky
3.1.4  Jacobi solver
3.2  Gauss-Seidal
3.2.1  SOR
3.2.2  Conjugate Gradient
" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Julia Version 1.1.0\n", "Commit 80516ca202 (2019-01-21 21:24 UTC)\n", "Platform Info:\n", " OS: macOS (x86_64-apple-darwin14.5.0)\n", " CPU: Intel(R) Core(TM) i7-6920HQ CPU @ 2.90GHz\n", " WORD_SIZE: 64\n", " LIBM: libopenlibm\n", " LLVM: libLLVM-6.0.1 (ORCJIT, skylake)\n", "Environment:\n", " JULIA_EDITOR = code\n" ] } ], "source": [ "versioninfo()" ] }, { "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): good for dense, small to moderate sized $\\mathbf{A}$ \n", " - Iterative methods (Jacobi, Gauss-Seidal, SOR, conjugate-gradient, GMRES): good for large, sparse, structured linear system, parallel computing, warm start\n", "\n", "\n", "## PageRank problem\n", "\n", "\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$. \n", "\n", "* [Larry Page](https://en.wikipedia.org/wiki/PageRank) imagines 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 on 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 a Markov chain on the space of $n$ pages. Stationary distribution of this Markov chain gives the ranks (probabilities) of each page.\n", "\n", "* Stationary distribution is the top left eigenvector of the transition matrix $\\mathbf{P}$ corresponding to eigenvalue 1. Equivalently it can be cast as a linear equation.\n", "$$\n", " (\\mathbf{I} - \\mathbf{P}^T) \\mathbf{x} = \\mathbf{0}.\n", "$$\n", "\n", "* Gene Golub: Largest matrix computation in world. \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": [ "## Jacobi method\n", "\n", "\n", "\n", "Solve $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$.\n", "\n", "* Jacobi iteration: \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", "* With splitting: $\\mathbf{A} = \\mathbf{L} + \\mathbf{D} + \\mathbf{U}$, Jacobi iteration can be written as\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", "* One round costs $2n^2$ flops with an **unstructured** $\\mathbf{A}$. Gain over GE/LU if converges in $o(n)$ iterations. 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)$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gauss-Seidel method\n", "\n", "\n", "\n", "\n", "\n", "* Gauss-Seidel (GS) iteration:\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", "* With splitting, $(\\mathbf{D} + \\mathbf{L}) \\mathbf{x}^{(t+1)} = - \\mathbf{U} \\mathbf{x}^{(t)} + \\mathbf{b}$, 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", "* GS converges for any $\\mathbf{x}^{(0)}$ for symmetric and pd $\\mathbf{A}$.\n", "\n", "* Convergence rate of Gauss-Seidel is the spectral radius of the $(\\mathbf{D} + \\mathbf{L})^{-1} \\mathbf{U}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Successive over-relaxation (SOR)\n", "\n", "* SOR: \n", "$$\n", "x_i^{(t+1)} = \\omega \\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) / a_{ii} + (1-\\omega) x_i^{(t)},\n", "$$\n", "i.e., \n", "$$\n", "\\mathbf{x}^{(t+1)} = (\\mathbf{D} + \\omega \\mathbf{L})^{-1} [(1-\\omega) \\mathbf{D} - \\omega \\mathbf{U}] \\mathbf{x}^{(t)} + (\\mathbf{D} + \\omega \\mathbf{L})^{-1} (\\mathbf{D} + \\mathbf{L})^{-1} \\omega \\mathbf{b}.\n", "$$\n", "\n", "* Need to pick $\\omega \\in [0,1]$ beforehand, with the goal of improving convergence rate." ] }, { "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", "* A UCLA invention! Hestenes and Steifel in 60s.\n", "\n", "* 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}$. \n", "\n", "[Kershaw's results](http://www.sciencedirect.com/science/article/pii/0021999178900980?via%3Dihub) for a fusion problem.\n", "\n", "| Method | Number of Iterations |\n", "|----------------------------------------|----------------------|\n", "| Gauss Seidel | 208,000 |\n", "| Block SOR methods | 765 |\n", "| Incomplete Cholesky conjugate gradient | 25 |" ] }, { "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": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "include group.jl for user defined matrix generators\n", "verify download of index files...\n", "used remote site is https://sparse.tamu.edu/?per_page=All\n", "populating internal database...\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 & 2833 \\\\\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 | 2833 |\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 2833\n", "\n", "MatrixMarket of \n", "–––––––––––– –––\n", "0 498" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using MatrixDepot\n", "\n", "mdinfo()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Array{String,1}:\n", " \"poisson\"\n", " \"wathen\" " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# List matrices that are positive definite and sparse:\n", "mdlist(:symmetric & :posdef & :sparse)" ] }, { "cell_type": "code", "execution_count": 4, "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": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Get help on Poisson matrix\n", "mdinfo(\"poisson\")" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "100×100 SparseArrays.SparseMatrixCSC{Float64,Int64} with 460 stored entries:\n", " [1 , 1] = 4.0\n", " [2 , 1] = -1.0\n", " [11 , 1] = -1.0\n", " [1 , 2] = -1.0\n", " [2 , 2] = 4.0\n", " [3 , 2] = -1.0\n", " [12 , 2] = -1.0\n", " [2 , 3] = -1.0\n", " [3 , 3] = 4.0\n", " [4 , 3] = -1.0\n", " [13 , 3] = -1.0\n", " [3 , 4] = -1.0\n", " ⋮\n", " [98 , 97] = -1.0\n", " [88 , 98] = -1.0\n", " [97 , 98] = -1.0\n", " [98 , 98] = 4.0\n", " [99 , 98] = -1.0\n", " [89 , 99] = -1.0\n", " [98 , 99] = -1.0\n", " [99 , 99] = 4.0\n", " [100, 99] = -1.0\n", " [90 , 100] = -1.0\n", " [99 , 100] = -1.0\n", " [100, 100] = 4.0" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Generate a Poisson matrix of dimension n=10\n", "A = matrixdepot(\"poisson\", 10)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[1m Sparsity Pattern\u001b[22m\n", "\u001b[90m ┌──────────────────────────────────────────┐\u001b[39m \n", " \u001b[90m1\u001b[39m\u001b[90m │\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠘\u001b[39m\u001b[34m⢄\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠉\u001b[39m\u001b[34m⢆\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[34m⠒\u001b[39m\u001b[34m⢄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠱\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠑\u001b[39m\u001b[34m⢢\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[34m⠣\u001b[39m\u001b[34m⢄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠑\u001b[39m\u001b[34m⠤\u001b[39m\u001b[34m⡀\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[34m⠱\u001b[39m\u001b[34m⣀\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠱\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠑\u001b[39m\u001b[34m⢄\u001b[39m\u001b[34m⡀\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[34m⠑\u001b[39m\u001b[34m⡄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠘\u001b[39m\u001b[34m⢄\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[34m⠈\u001b[39m\u001b[34m⠑\u001b[39m\u001b[34m⢄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠱\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠉\u001b[39m\u001b[34m⢆\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[34m⠈\u001b[39m\u001b[34m⠒\u001b[39m\u001b[34m⢄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠑\u001b[39m\u001b[34m⢢\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[34m⠣\u001b[39m\u001b[34m⢄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⣤\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠑\u001b[39m\u001b[34m⠤\u001b[39m\u001b[34m⡀\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[34m⠱\u001b[39m\u001b[34m⣀\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠑\u001b[39m\u001b[34m⢄\u001b[39m\u001b[34m⡀\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[34m⠑\u001b[39m\u001b[34m⡄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⣤\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠘\u001b[39m\u001b[34m⢄\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[34m⠈\u001b[39m\u001b[34m⠑\u001b[39m\u001b[34m⢄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠉\u001b[39m\u001b[34m⢆\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[34m⠈\u001b[39m\u001b[34m⠒\u001b[39m\u001b[34m⢄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⣤\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠑\u001b[39m\u001b[34m⢢\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[34m⠣\u001b[39m\u001b[34m⢄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠑\u001b[39m\u001b[34m⠤\u001b[39m\u001b[34m⡀\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[34m⠱\u001b[39m\u001b[34m⣀\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⢆\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠑\u001b[39m\u001b[34m⢄\u001b[39m\u001b[34m⡀\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[34m⠑\u001b[39m\u001b[34m⡄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠘\u001b[39m\u001b[34m⢄\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[34m⠈\u001b[39m\u001b[34m⠑\u001b[39m\u001b[34m⢄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⢆\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠉\u001b[39m\u001b[34m⢆\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[34m⠈\u001b[39m\u001b[34m⠒\u001b[39m\u001b[34m⢄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠑\u001b[39m\u001b[34m⢢\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[34m⠣\u001b[39m\u001b[34m⢄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⢆\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[34m⠑\u001b[39m\u001b[34m⠤\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[34m⠱\u001b[39m\u001b[34m⣀\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠑\u001b[39m\u001b[34m⡄\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[90m│\u001b[39m \n", "\u001b[90m └──────────────────────────────────────────┘\u001b[39m \n", "\u001b[90m 1\u001b[39m\u001b[90m \u001b[39m\u001b[90m 100\u001b[39m\n", "\u001b[0m nz = 460" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using UnicodePlots\n", "spy(A)" ] }, { "cell_type": "code", "execution_count": 7, "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": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Get help on Wathen matrix\n", "mdinfo(\"wathen\")" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "96×96 SparseArrays.SparseMatrixCSC{Float64,Int64} with 1256 stored entries:\n", " [1 , 1] = 1.81938\n", " [2 , 1] = -1.81938\n", " [3 , 1] = 0.606459\n", " [12, 1] = -1.81938\n", " [13, 1] = -2.42583\n", " [18, 1] = 0.606459\n", " [19, 1] = -2.42583\n", " [20, 1] = 0.909688\n", " [1 , 2] = -1.81938\n", " [2 , 2] = 9.70334\n", " [3 , 2] = -1.81938\n", " [12, 2] = 6.06459\n", " ⋮\n", " [85, 95] = 20.6997\n", " [94, 95] = -6.2099\n", " [95, 95] = 33.1194\n", " [96, 95] = -6.2099\n", " [77, 96] = 3.10495\n", " [78, 96] = -8.27986\n", " [79, 96] = 2.06997\n", " [84, 96] = -8.27986\n", " [85, 96] = -6.2099\n", " [94, 96] = 2.06997\n", " [95, 96] = -6.2099\n", " [96, 96] = 6.2099" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Generate a Wathen matrix of dimension n=5\n", "A = matrixdepot(\"wathen\", 5)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[1m Sparsity Pattern\u001b[22m\n", "\u001b[90m ┌──────────────────────────────────────────┐\u001b[39m \n", " \u001b[90m1\u001b[39m\u001b[90m │\u001b[39m\u001b[35m⠿\u001b[39m\u001b[35m⣧\u001b[39m\u001b[35m⡄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[35m⢿\u001b[39m\u001b[35m⡄\u001b[39m\u001b[35m⠸\u001b[39m\u001b[35m⢿\u001b[39m\u001b[35m⣤\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[35m⠉\u001b[39m\u001b[35m⠿\u001b[39m\u001b[35m⣧\u001b[39m\u001b[35m⣀\u001b[39m\u001b[34m⠈\u001b[39m\u001b[35m⣧\u001b[39m\u001b[34m⡀\u001b[39m\u001b[31m⠈\u001b[39m\u001b[35m⠹\u001b[39m\u001b[35m⣧\u001b[39m\u001b[35m⣄\u001b[39m\u001b[31m⡀\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[35m⣤\u001b[39m\u001b[35m⣄\u001b[39m\u001b[34m⡀\u001b[39m\u001b[35m⠘\u001b[39m\u001b[35m⠛\u001b[39m\u001b[31m⣤\u001b[39m\u001b[35m⡘\u001b[39m\u001b[35m⢣\u001b[39m\u001b[35m⣤\u001b[39m\u001b[35m⣀\u001b[39m\u001b[0m⠀\u001b[35m⠛\u001b[39m\u001b[35m⠃\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[35m⣀\u001b[39m\u001b[35m⡉\u001b[39m\u001b[35m⠉\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⠶\u001b[39m\u001b[35m⣈\u001b[39m\u001b[31m⠻\u001b[39m\u001b[31m⢆\u001b[39m\u001b[35m⣈\u001b[39m\u001b[35m⠉\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⠷\u001b[39m\u001b[34m⢆\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[35m⣀\u001b[39m\u001b[31m⡀\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[35m⠛\u001b[39m\u001b[35m⣷\u001b[39m\u001b[35m⣆\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[35m⢻\u001b[39m\u001b[35m⡆\u001b[39m\u001b[35m⠘\u001b[39m\u001b[35m⢻\u001b[39m\u001b[35m⣶\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[35m⢸\u001b[39m\u001b[35m⣆\u001b[39m\u001b[0m⠀\u001b[35m⠛\u001b[39m\u001b[35m⣷\u001b[39m\u001b[35m⣀\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[35m⠉\u001b[39m\u001b[35m⢿\u001b[39m\u001b[35m⣤\u001b[39m\u001b[0m⠀\u001b[35m⢿\u001b[39m\u001b[35m⡄\u001b[39m\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠿\u001b[39m\u001b[35m⣧\u001b[39m\u001b[35m⡄\u001b[39m\u001b[35m⠹\u001b[39m\u001b[35m⣧\u001b[39m\u001b[0m⠀\u001b[31m⠈\u001b[39m\u001b[35m⠹\u001b[39m\u001b[35m⢿\u001b[39m\u001b[35m⡄\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[31m⠈\u001b[39m\u001b[35m⠉\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[34m⠑\u001b[39m\u001b[35m⠲\u001b[39m\u001b[35m⢶\u001b[39m\u001b[35m⣄\u001b[39m\u001b[35m⡉\u001b[39m\u001b[31m⠱\u001b[39m\u001b[31m⣦\u001b[39m\u001b[35m⡉\u001b[39m\u001b[35m⠒\u001b[39m\u001b[35m⢶\u001b[39m\u001b[35m⣤\u001b[39m\u001b[35m⣈\u001b[39m\u001b[35m⠁\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[35m⢠\u001b[39m\u001b[35m⣤\u001b[39m\u001b[0m⠀\u001b[35m⠉\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⢣\u001b[39m\u001b[31m⠈\u001b[39m\u001b[31m⠛\u001b[39m\u001b[35m⣤\u001b[39m\u001b[35m⡄\u001b[39m\u001b[34m⠈\u001b[39m\u001b[35m⠙\u001b[39m\u001b[35m⢣\u001b[39m\u001b[35m⡄\u001b[39m\u001b[0m⠀\u001b[35m⢠\u001b[39m\u001b[35m⡄\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[31m⠈\u001b[39m\u001b[35m⠙\u001b[39m\u001b[35m⢻\u001b[39m\u001b[35m⣆\u001b[39m\u001b[31m⡀\u001b[39m\u001b[35m⠘\u001b[39m\u001b[35m⣷\u001b[39m\u001b[34m⡀\u001b[39m\u001b[35m⠉\u001b[39m\u001b[35m⢻\u001b[39m\u001b[35m⣶\u001b[39m\u001b[35m⣀\u001b[39m\u001b[34m⠈\u001b[39m\u001b[35m⢻\u001b[39m\u001b[35m⣆\u001b[39m\u001b[35m⠈\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⣷\u001b[39m\u001b[35m⣆\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[35m⠛\u001b[39m\u001b[35m⠷\u001b[39m\u001b[35m⠆\u001b[39m\u001b[35m⠘\u001b[39m\u001b[35m⠷\u001b[39m\u001b[35m⣀\u001b[39m\u001b[34m⡀\u001b[39m\u001b[35m⠘\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⢆\u001b[39m\u001b[31m⡀\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⢆\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[35m⠻\u001b[39m\u001b[35m⠶\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[35m⠉\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⢶\u001b[39m\u001b[35m⣤\u001b[39m\u001b[35m⡈\u001b[39m\u001b[31m⠻\u001b[39m\u001b[31m⣦\u001b[39m\u001b[35m⡈\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⠷\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⣀\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[35m⠶\u001b[39m\u001b[35m⣦\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠱\u001b[39m\u001b[35m⣦\u001b[39m\u001b[31m⠈\u001b[39m\u001b[35m⠱\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡄\u001b[39m\u001b[34m⠈\u001b[39m\u001b[35m⠉\u001b[39m\u001b[35m⢶\u001b[39m\u001b[35m⡄\u001b[39m\u001b[35m⠰\u001b[39m\u001b[35m⢶\u001b[39m\u001b[35m⣤\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[35m⠹\u001b[39m\u001b[35m⢿\u001b[39m\u001b[35m⣤\u001b[39m\u001b[35m⡀\u001b[39m\u001b[35m⠹\u001b[39m\u001b[35m⣧\u001b[39m\u001b[34m⡀\u001b[39m\u001b[35m⠉\u001b[39m\u001b[35m⠿\u001b[39m\u001b[35m⣧\u001b[39m\u001b[35m⣀\u001b[39m\u001b[34m⠈\u001b[39m\u001b[35m⢿\u001b[39m\u001b[35m⡄\u001b[39m\u001b[31m⠈\u001b[39m\u001b[35m⠹\u001b[39m\u001b[35m⣧\u001b[39m\u001b[35m⣄\u001b[39m\u001b[31m⡀\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[35m⠘\u001b[39m\u001b[35m⠃\u001b[39m\u001b[0m⠀\u001b[35m⠘\u001b[39m\u001b[35m⢣\u001b[39m\u001b[35m⣄\u001b[39m\u001b[34m⡀\u001b[39m\u001b[35m⠘\u001b[39m\u001b[35m⠛\u001b[39m\u001b[31m⣤\u001b[39m\u001b[31m⡀\u001b[39m\u001b[35m⢣\u001b[39m\u001b[35m⣤\u001b[39m\u001b[35m⣀\u001b[39m\u001b[0m⠀\u001b[35m⠛\u001b[39m\u001b[35m⠃\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[35m⢀\u001b[39m\u001b[35m⡉\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⠷\u001b[39m\u001b[35m⠤\u001b[39m\u001b[35m⣈\u001b[39m\u001b[31m⠻\u001b[39m\u001b[31m⢆\u001b[39m\u001b[35m⣈\u001b[39m\u001b[35m⠙\u001b[39m\u001b[35m⠷\u001b[39m\u001b[35m⠦\u001b[39m\u001b[34m⢄\u001b[39m\u001b[34m⡀\u001b[39m\u001b[0m⠀\u001b[35m⣀\u001b[39m\u001b[31m⡀\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[35m⠘\u001b[39m\u001b[35m⣷\u001b[39m\u001b[35m⣆\u001b[39m\u001b[31m⡀\u001b[39m\u001b[0m⠀\u001b[35m⢻\u001b[39m\u001b[35m⣆\u001b[39m\u001b[35m⠘\u001b[39m\u001b[35m⢻\u001b[39m\u001b[35m⣶\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[35m⠘\u001b[39m\u001b[35m⣷\u001b[39m\u001b[0m⠀\u001b[35m⠛\u001b[39m\u001b[35m⣷\u001b[39m\u001b[35m⣀\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[35m⠉\u001b[39m\u001b[35m⢿\u001b[39m\u001b[35m⣤\u001b[39m\u001b[0m⠀\u001b[35m⠹\u001b[39m\u001b[35m⡇\u001b[39m\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠿\u001b[39m\u001b[35m⣧\u001b[39m\u001b[35m⡄\u001b[39m\u001b[35m⠸\u001b[39m\u001b[35m⣧\u001b[39m\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠹\u001b[39m\u001b[35m⢿\u001b[39m\u001b[35m⣤\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[31m⠈\u001b[39m\u001b[35m⠉\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[34m⠱\u001b[39m\u001b[35m⢶\u001b[39m\u001b[35m⣤\u001b[39m\u001b[35m⣀\u001b[39m\u001b[35m⡉\u001b[39m\u001b[31m⠱\u001b[39m\u001b[31m⣦\u001b[39m\u001b[35m⡉\u001b[39m\u001b[35m⠶\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⣀\u001b[39m\u001b[35m⣈\u001b[39m\u001b[35m⠉\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[35m⢠\u001b[39m\u001b[35m⣤\u001b[39m\u001b[0m⠀\u001b[35m⠉\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⢣\u001b[39m\u001b[35m⡌\u001b[39m\u001b[31m⠛\u001b[39m\u001b[35m⣤\u001b[39m\u001b[35m⡄\u001b[39m\u001b[34m⠈\u001b[39m\u001b[35m⠙\u001b[39m\u001b[35m⠛\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[31m⠈\u001b[39m\u001b[35m⠙\u001b[39m\u001b[35m⢻\u001b[39m\u001b[35m⣆\u001b[39m\u001b[31m⡀\u001b[39m\u001b[34m⠈\u001b[39m\u001b[35m⢻\u001b[39m\u001b[34m⡀\u001b[39m\u001b[35m⠉\u001b[39m\u001b[35m⢻\u001b[39m\u001b[35m⣶\u001b[39m\u001b[35m⣀\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[35m⠛\u001b[39m\u001b[35m⣷\u001b[39m\u001b[35m⡆\u001b[39m\u001b[35m⠘\u001b[39m\u001b[35m⣷\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[35m⠘\u001b[39m\u001b[35m⢻\u001b[39m\u001b[35m⣶\u001b[39m\u001b[90m│\u001b[39m \n", "\u001b[90m └──────────────────────────────────────────┘\u001b[39m \n", "\u001b[90m 1\u001b[39m\u001b[90m \u001b[39m\u001b[90m 96\u001b[39m\n", "\u001b[0m nz = 1256" ] }, "execution_count": 9, "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": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "typeof(A) = SparseArrays.SparseMatrixCSC{Float64,Int64}\n", "typeof(Afull) = Array{Float64,2}\n" ] }, { "data": { "text/plain": [ "0.000496" ] }, "execution_count": 10, "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": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[1m Sparsity Pattern\u001b[22m\n", "\u001b[90m ┌──────────────────────────────────────────┐\u001b[39m \n", " \u001b[90m1\u001b[39m\u001b[90m │\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[34m⡀\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[34m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[90m│\u001b[39m \n", "\u001b[90m └──────────────────────────────────────────┘\u001b[39m \n", "\u001b[90m 1\u001b[39m\u001b[90m \u001b[39m\u001b[90m 10000\u001b[39m\n", "\u001b[0m nz = 49600" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "spy(A)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(873768, 800000040)" ] }, "execution_count": 12, "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": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 78.20 KiB\n", " allocs estimate: 2\n", " --------------\n", " minimum time: 35.722 ms (0.00% GC)\n", " median time: 43.154 ms (0.00% GC)\n", " mean time: 44.512 ms (0.00% GC)\n", " maximum time: 51.161 ms (0.00% GC)\n", " --------------\n", " samples: 113\n", " evals/sample: 1" ] }, "execution_count": 13, "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": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 78.20 KiB\n", " allocs estimate: 2\n", " --------------\n", " minimum time: 51.993 μs (0.00% GC)\n", " median time: 93.657 μs (0.00% GC)\n", " mean time: 105.492 μs (10.11% GC)\n", " maximum time: 4.137 ms (97.35% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 1" ] }, "execution_count": 14, "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": 15, "metadata": {}, "outputs": [], "source": [ "# record the Cholesky solution\n", "xchol = cholesky(Afull) \\ b;" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 763.02 MiB\n", " allocs estimate: 8\n", " --------------\n", " minimum time: 3.505 s (3.20% GC)\n", " median time: 3.508 s (1.77% GC)\n", " mean time: 3.508 s (1.77% GC)\n", " maximum time: 3.511 s (0.34% GC)\n", " --------------\n", " samples: 2\n", " evals/sample: 1" ] }, "execution_count": 16, "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": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "norm(xjacobi - xchol) = 533.315106586665\n" ] }, { "data": { "text/plain": [ "533.315106586665" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xjacobi = jacobi(A, b)\n", "@show norm(xjacobi - xchol)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reading [documentation](https://juliamath.github.io/IterativeSolvers.jl/dev/linear_systems/stationary/#Jacobi-1) we found 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": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "norm(xjacobi - xchol) = 0.00012526024863889176\n" ] }, { "data": { "text/plain": [ "0.00012526024863889176" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xjacobi = jacobi(A, b, maxiter=30000)\n", "@show norm(xjacobi - xchol)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 234.72 KiB\n", " allocs estimate: 9\n", " --------------\n", " minimum time: 2.261 s (0.00% GC)\n", " median time: 2.348 s (0.00% GC)\n", " mean time: 2.329 s (0.00% GC)\n", " maximum time: 2.377 s (0.00% GC)\n", " --------------\n", " samples: 3\n", " evals/sample: 1" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@benchmark jacobi($A, $b, maxiter=30000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gauss-Seidal" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "norm(xgs - xchol) = 0.0001245846635425542\n" ] }, { "data": { "text/plain": [ "0.0001245846635425542" ] }, "execution_count": 20, "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": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 156.55 KiB\n", " allocs estimate: 8\n", " --------------\n", " minimum time: 1.801 s (0.00% GC)\n", " median time: 1.853 s (0.00% GC)\n", " mean time: 1.849 s (0.00% GC)\n", " maximum time: 1.892 s (0.00% GC)\n", " --------------\n", " samples: 3\n", " evals/sample: 1" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@benchmark gauss_seidel($A, $b, maxiter=15000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### SOR" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "norm(xsor - xchol) = 0.0022945523203956853\n" ] }, { "data": { "text/plain": [ "0.0022945523203956853" ] }, "execution_count": 22, "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": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 547.27 KiB\n", " allocs estimate: 20010\n", " --------------\n", " minimum time: 1.433 s (0.00% GC)\n", " median time: 1.441 s (0.00% GC)\n", " mean time: 1.443 s (0.00% GC)\n", " maximum time: 1.455 s (0.00% GC)\n", " --------------\n", " samples: 4\n", " evals/sample: 1" ] }, "execution_count": 23, "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": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "norm(xcg - xchol) = 1.1005386920768871e-5\n" ] }, { "data": { "text/plain": [ "1.1005386920768871e-5" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# conjugate gradient\n", "xcg = cg(A, b)\n", "@show norm(xcg - xchol)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 391.89 KiB\n", " allocs estimate: 23\n", " --------------\n", " minimum time: 19.614 ms (0.00% GC)\n", " median time: 23.401 ms (0.00% GC)\n", " mean time: 23.440 ms (0.18% GC)\n", " maximum time: 28.527 ms (0.00% GC)\n", " --------------\n", " samples: 214\n", " evals/sample: 1" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@benchmark cg($A, $b)" ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.1.0", "language": "julia", "name": "julia-1.1" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.1.0" }, "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": 2 }