diff --git a/ipython/lgmres.ipynb b/ipython/lgmres.ipynb new file mode 100644 index 0000000..ac9e577 --- /dev/null +++ b/ipython/lgmres.ipynb @@ -0,0 +1,111 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "f43700ac", + "metadata": {}, + "source": [ + "# LGMRES" + ] + }, + { + "cell_type": "markdown", + "id": "6d17c3dc", + "metadata": {}, + "source": [ + "Example showing how LGMRES avoids some problems in the convergence of restarted GMRES." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9f4264f8", + "metadata": {}, + "outputs": [], + "source": [ + "import scipy.sparse.linalg as la\n", + "import scipy.io as io\n", + "import numpy as np\n", + "import sys\n", + "\n", + "#problem = \"SPARSKIT/drivcav/e05r0100\"\n", + "problem = \"SPARSKIT/drivcav/e05r0200\"\n", + "#problem = \"Harwell-Boeing/sherman/sherman1\"\n", + "#problem = \"misc/hamm/add32\"\n", + "\n", + "mm = np.lib._datasource.Repository('https://math.nist.gov/pub/MatrixMarket2/')\n", + "f = mm.open(f'{problem}.mtx.gz')\n", + "Am = io.mmread(f).tocsr()\n", + "f.close()\n", + "\n", + "f = mm.open(f'{problem}_rhs1.mtx.gz')\n", + "b = np.array(io.mmread(f)).ravel()\n", + "f.close()\n", + "\n", + "bnorm = np.linalg.norm(b)\n", + "count = [0]\n", + "\n", + "\n", + "def matvec(v):\n", + " count[0] += 1\n", + " sys.stderr.write(f\"{count[0]}\\r\")\n", + " return Am@v\n", + "\n", + "\n", + "A = la.LinearOperator(matvec=matvec, shape=Am.shape, dtype=Am.dtype)\n", + "\n", + "M = 100\n", + "\n", + "print(\"MatrixMarket problem %s\" % problem)\n", + "print(f\"Invert {Am.shape[0]} x {Am.shape[1]} matrix; nnz = {Am.nnz}\")\n", + "\n", + "count[0] = 0\n", + "x0, info = la.gmres(A, b, restrt=M, tol=1e-14)\n", + "count_0 = count[0]\n", + "err0 = np.linalg.norm(Am@x0 - b) / bnorm\n", + "print(f\"GMRES({M}): {count_0} matvecs, relative residual: {err0}\")\n", + "if info != 0:\n", + " print(\"Didn't converge\")\n", + "\n", + "count[0] = 0\n", + "x1, info = la.lgmres(A, b, inner_m=M-6*2, outer_k=6, tol=1e-14)\n", + "count_1 = count[0]\n", + "err1 = np.linalg.norm(Am@x1 - b) / bnorm\n", + "print(f\"LGMRES({M - 2*6}, 6) [same memory req.]: {count_1} \"\n", + " f\"matvecs, relative residual: {err1}\")\n", + "if info != 0:\n", + " print(\"Didn't converge\")\n", + "\n", + "count[0] = 0\n", + "x2, info = la.lgmres(A, b, inner_m=M-6, outer_k=6, tol=1e-14)\n", + "count_2 = count[0]\n", + "err2 = np.linalg.norm(Am@x2 - b) / bnorm\n", + "print(f\"LGMRES({M - 6}, 6) [same subspace size]: {count_2} \"\n", + " f\"matvecs, relative residual: {err2}\")\n", + "if info != 0:\n", + " print(\"Didn't converge\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}