move LGMRES demo over from scipy (#50)

* move LGMRES demo over from scipy

* remove output and fix typo

Author: j-bowhay

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
+}