{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.2 Programming an eigenvalue solver\n",
    "\n",
    "We solve the generalized eigenvalue problem\n",
    "\n",
    "$$\n",
    "A u = \\lambda M u,\n",
    "$$\n",
    "\n",
    "where $A$ comes from $\\int \\nabla u \\nabla v$, and $M$ from $\\int u v$, on the space $H_0^1$.\n",
    "\n",
    "This tutorial shows how to implement linear algebra algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "\n",
    "import scipy.linalg\n",
    "\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We setup a stiffness matrix $A$ and a mass matrix $M$, and declare a preconditioner for $A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=4, dirichlet=\".*\")\n",
    "u,v = fes.TnT()\n",
    "\n",
    "a = BilinearForm(grad(u)*grad(v)*dx)\n",
    "pre = Preconditioner(a, \"multigrid\")\n",
    "a.Assemble()\n",
    "\n",
    "m = BilinearForm(u*v*dx).Assemble()\n",
    "\n",
    "u = GridFunction(fes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inverse iteration is\n",
    "\n",
    "$$\n",
    "u_{n+1} = A^{-1} M u_n,\n",
    "$$\n",
    "\n",
    "where the Rayleigh quotient\n",
    "$$\n",
    "\\rho_n = \\frac{\\left \\langle A u_n, u_n\\right \\rangle}{\\left \n",
    "\\langle M u_n, u_n\\right \\rangle}\n",
    "$$\n",
    "\n",
    "converges to the smallest eigenvalue $\\lambda_1$, with rate of convergence $\\lambda_1 / \\lambda_2,$ where $\\lambda_2$ is the next smallest eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The preconditioned inverse iteration (PINVIT), see [Knyazef+Neymeyr], replaces $A^{-1}$ by an approximate inverse $C^{-1}$:\n",
    "\n",
    "$$\n",
    "\\rho_n = \\frac{\\left \\langle A u_n, u_n\\right \\rangle}{\\left \\langle M u_n, u_n\\right \\rangle} \\\\\n",
    "w_n = C^{-1} (A u_n - \\rho_n M u_n) \\\\\n",
    "u_{n+1} = u_n + \\alpha w_n\n",
    "$$\n",
    "\n",
    "The optimal step-size $\\alpha$ is found by minimizing the Rayleigh-quotient on a two-dimensional space:\n",
    "\n",
    "$$\n",
    "u_{n+1} = \\operatorname{arg} \\min_{v \\in \\{ u_n, w_n\\}} \\frac{\\left \\langle A v, v\\right \\rangle}{\\left \\langle M v, v\\right \\rangle} \n",
    "$$\n",
    "\n",
    "This minimization problem can be solved by a small eigenvalue problem\n",
    "\n",
    "$$\n",
    "a y = \\lambda m y\n",
    "$$\n",
    "\n",
    "with matrices\n",
    "\n",
    "$$\n",
    "a = \n",
    "\\left( \\begin{array}{cc}\n",
    "\\left \\langle  A u_n, u_n \\right \\rangle &\n",
    "\\left \\langle  A u_n, w_n \\right \\rangle \\\\\n",
    "\\left \\langle  A w_n, u_n \\right \\rangle &\n",
    "\\left \\langle  A w_n, w_n \\right \\rangle \n",
    "\\end{array} \\right), \\quad\n",
    "m = \n",
    "\\left( \\begin{array}{cc}\n",
    "\\left \\langle  M u_n, u_n \\right \\rangle &\n",
    "\\left \\langle  M u_n, w_n \\right \\rangle \\\\\n",
    "\\left \\langle  M w_n, u_n \\right \\rangle &\n",
    "\\left \\langle  M w_n, w_n \\right \\rangle \n",
    "\\end{array} \\right).\n",
    "$$\n",
    "\n",
    "Then, the new iterate is\n",
    "\n",
    "$$\n",
    "u_{n+1} = y_1 u_n + y_2 w_n\n",
    "$$\n",
    "\n",
    "where $y$ is the eigenvector corresponding to the smaller eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation in NGSolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we create some help vectors. `CreateVector` creates new vectors of the same format as the existing vector, i.e., same dimension, same real/ complex type, same entry-size, and same MPI-parallel distribution if any."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r, w, Mu, Au, Mw, Aw = u.vec.CreateVectors(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we pick a random initial vector, which is  zeroed on the Dirichlet boundary.\n",
    "\n",
    "Below, the `FV` method (short for `FlatVector`)  lets us access the abstract vector's  linear memory chunk, which in turn provides a \"numpy\" view of the memory. The projector clears the entries at the Dirichlet boundary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# r.FV().NumPy()[:] = random.rand(fes.ndof)\n",
    "r.SetRandom()\n",
    "u.vec.data = Projector(fes.FreeDofs(), True) * r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we run the PINVIT algorithm. Note that the small matrices $a$ and $m$ defined above are called `asmall` and `msmall` below. They are of type `Matrix`, a class provided by  NGSolve for  dense matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(20):\n",
    "    Au.data = a.mat * u.vec\n",
    "    Mu.data = m.mat * u.vec\n",
    "    auu = InnerProduct(Au, u.vec)\n",
    "    muu = InnerProduct(Mu, u.vec)\n",
    "    # Rayleigh quotient\n",
    "    lam = auu/muu\n",
    "    print (lam / (pi**2))\n",
    "    # residual\n",
    "    r.data = Au - lam * Mu\n",
    "    w.data = pre.mat * r.data\n",
    "    # w.data = 1/Norm(w) * w\n",
    "    w *= 1/Norm(w)\n",
    "    Aw.data = a.mat * w\n",
    "    Mw.data = m.mat * w\n",
    "\n",
    "    # setup and solve 2x2 small eigenvalue problem\n",
    "    asmall = Matrix(2,2)\n",
    "    asmall[0,0] = auu\n",
    "    asmall[0,1] = asmall[1,0] = InnerProduct(Au, w)\n",
    "    asmall[1,1] = InnerProduct(Aw, w)\n",
    "    msmall = Matrix(2,2)\n",
    "    msmall[0,0] = muu\n",
    "    msmall[0,1] = msmall[1,0] = InnerProduct(Mu, w)\n",
    "    msmall[1,1] = InnerProduct(Mw, w)\n",
    "    # print (\"asmall =\", asmall, \", msmall = \", msmall)\n",
    "    \n",
    "    \n",
    "    eval,evec = scipy.linalg.eigh(a=asmall, b=msmall)\n",
    "    # print (eval, evec)\n",
    "    u.vec.data = float(evec[0,0]) * u.vec + float(evec[1,0]) * w\n",
    "    \n",
    "Draw (u);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simultaneous iteration for several eigenvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the steps for extending the above to `num` vectors.\n",
    "\n",
    "Declare a `GridFunction` with multiple components to store several eigenvectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num = 5\n",
    "u = GridFunction(fes, multidim=num)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create list of help vectors, and a set of random initial vectors in `u`, with zero boundary conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = u.vec.CreateVector()\n",
    "Av = u.vec.CreateVector()\n",
    "Mv = u.vec.CreateVector()\n",
    "\n",
    "vecs = []\n",
    "for v in u.vec.CreateVectors(2*num):\n",
    "    vecs.append(v)\n",
    "\n",
    "for v in u.vecs:\n",
    "    r.SetRandom()\n",
    "    v.data = Projector(fes.FreeDofs(), True) * r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute `num` residuals, and solve a small eigenvalue problem on a 2 $\\times$ `num` dimensional space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "asmall = Matrix(2*num, 2*num)\n",
    "msmall = Matrix(2*num, 2*num)\n",
    "lams = num * [1]\n",
    "\n",
    "for i in range(20):\n",
    "    \n",
    "    for j in range(num):\n",
    "        vecs[j].data = u.vecs[j]\n",
    "        r.data = a.mat * vecs[j] - lams[j] * m.mat * vecs[j]\n",
    "        vecs[num+j].data = pre.mat * r\n",
    "\n",
    "    for j in range(2*num):\n",
    "        Av.data = a.mat * vecs[j]\n",
    "        Mv.data = m.mat * vecs[j]\n",
    "        for k in range(2*num):\n",
    "            asmall[j,k] = InnerProduct(Av, vecs[k])\n",
    "            msmall[j,k] = InnerProduct(Mv, vecs[k])\n",
    "\n",
    "    ev,evec = scipy.linalg.eigh(a=asmall, b=msmall)\n",
    "    lams[:] = ev[0:num]\n",
    "    print (i, \":\", [lam/pi**2 for lam in lams])\n",
    "    \n",
    "    for j in range(num):\n",
    "        u.vecs[j][:] = 0.0\n",
    "        for k in range(2*num):\n",
    "            u.vecs[j].data += float(evec[k,j]) * vecs[k]\n",
    "            \n",
    "Draw (u);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *multidim-component* select in the *Visualization* dialog box allows to switch between the components of the multidim-GridFunction (Netgen-gui only)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation using MultiVector\n",
    "\n",
    "The simultaneous iteration can be optimized by using `MultiVector`s. These are arrays of vectors of the same format. You can think of a MultiVector with m components of vector size n as an $n \\times m$ matrix. \n",
    "\n",
    "  - a `MultiVector` consisting of num vectors of the same format as an existing vector vec is created via\n",
    "MultiVector(vec, num). \n",
    "  - we can iterate over the components of a MultiVector, and the bracket operator allows to access a subset of vectors\n",
    "  - linear operator application is optimized for MultiVector\n",
    "  - vector operations are optimized and called as mv * densematrix: \n",
    "    $x = y * mat$ results in x[i] = sum_j y[j] * mat[j,i]  (where x and y are 'MultiVector's, and mat is a dense matrix) \n",
    "  - pair-wise inner products of two MultiVectors is available, the result is a dense matrix: `InnerProduct(x,y)[i,j] = InnerProduct(x_j, y_i)`. Before ngsolve-version  6.2.2301-242, the result was the transposed matrix.\n",
    "  - mv.Orthogonalize() uses modified Gram-Schmidt to orthogonalize the vectors. Optionally, a matrix defining the inner product can be provided.\n",
    "  - with mv.Append(vec) we can add another vector to the array of vectos. A new vector is created, and the values are copied\n",
    "  - mv.AppendOrthogonalize(vec) appends a new vector, and orthogonalizes it against the existing vectors, which are assumed to be orthogonal. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uvecs = MultiVector(u.vec, num)\n",
    "vecs = MultiVector(u.vec, 2*num)\n",
    "\n",
    "for v in vecs[0:num]:\n",
    "    v.SetRandom()\n",
    "uvecs[:] = pre * vecs[0:num]\n",
    "lams = Vector(num * [1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(20):\n",
    "    vecs[0:num] = a.mat * uvecs - (m.mat * uvecs).Scale (lams)\n",
    "    vecs[num:2*num] = pre * vecs[0:num]\n",
    "    vecs[0:num] = uvecs\n",
    "        \n",
    "    vecs.Orthogonalize() # m.mat)\n",
    "\n",
    "    asmall = InnerProduct (vecs, a.mat*vecs)\n",
    "    msmall = InnerProduct (vecs, m.mat*vecs)\n",
    "    \n",
    "    ev,evec = scipy.linalg.eigh(a=asmall, b=msmall)\n",
    "    lams = Vector(ev[0:num])\n",
    "    print (i, \":\", [l/pi**2 for l in lams])\n",
    "\n",
    "    uvecs[:] = vecs * Matrix(evec[:,0:num])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The operations are implemented using late evaluation. The operations return expressions, and computation is done within the assignment operator. The advantage is to avoid dynamic allocation. An exception is InnerProduct, which allows an expression in the second argument (and then needs vector allocation in every call)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}