Standard symmetric eigenproblem for the Laplacian operator in 1-D
=================================================================

This tutorial is intended for basic use of slepc4py. For more advanced use,
the reader is referred to SLEPc tutorials as well as to slepc4py reference
documentation.

The source code for this demo can be `downloaded here
<../_static/ex1.py>`__

The first thing to do is initialize the libraries. This is normally not
required, as it is done automatically at import time. However, if you want to
gain access to the facilities for accessing command-line options, the
following lines must be executed by the main script prior to any petsc4py or
slepc4py calls:

::

  import sys, slepc4py
  slepc4py.init(sys.argv)

Next, we have to import the relevant modules. Normally, both PETSc and SLEPc
modules have to be imported in all slepc4py programs. It may be useful to
import NumPy as well:

::

  from petsc4py import PETSc
  from slepc4py import SLEPc
  import numpy

At this point, we can use any petsc4py and slepc4py operations. For instance,
the following lines allow the user to specify an integer command-line
argument n with a default value of 30 (see the next section for example usage
of command-line options):

::

  opts = PETSc.Options()
  n = opts.getInt('n', 30)

It is necessary to build a matrix to define an eigenproblem (or two in the
case of generalized eigenproblems). The following fragment of code creates
the matrix object and then fills the non-zero elements one by one. The matrix
of this particular example is tridiagonal, with value 2 in the diagonal, and
-1 in off-diagonal positions. See petsc4py documentation for details about
matrix objects:

::

  A = PETSc.Mat(); A.create()
  A.setSizes([n, n])
  A.setFromOptions()

  rstart, rend = A.getOwnershipRange()

  # first row
  if rstart == 0:
    A[0, :2] = [2, -1]
    rstart += 1
  # last row
  if rend == n:
    A[n-1, -2:] = [-1, 2]
    rend -= 1
  # other rows
  for i in range(rstart, rend):
    A[i, i-1:i+2] = [-1, 2, -1]

  A.assemble()

The solver object is created in a similar way as other objects in petsc4py:

::

  E = SLEPc.EPS(); E.create()

Once the object is created, the eigenvalue problem must be specified. At
least one matrix must be provided. The problem type must be indicated as
well, in this case it is HEP (Hermitian eigenvalue problem). Apart from
these, other settings could be provided here (for instance, the tolerance for
the computation). After all options have been set, the user should call the
`setFromOptions() <EPS.setFromOptions()>` operation, so that any options
specified at run time in the command line are passed to the solver object:

::

  E.setOperators(A)
  E.setProblemType(SLEPc.EPS.ProblemType.HEP)

  history = []
  def monitor(eps, its, nconv, eig, err):
      if nconv<len(err): history.append(err[nconv])
  E.setMonitor(monitor)

  E.setFromOptions()

After that, the `solve() <EPS.solve()>` method will run the selected
eigensolver, keeping the solution stored internally:

::

  E.solve()

Once the computation has finished, we are ready to print the results. First,
some informative data can be retrieved from the solver object:

::

  Print = PETSc.Sys.Print

  Print()
  Print("******************************")
  Print("*** SLEPc Solution Results ***")
  Print("******************************")
  Print()

  its = E.getIterationNumber()
  Print( "Number of iterations of the method: %d" % its )

  eps_type = E.getType()
  Print( "Solution method: %s" % eps_type )

  nev, ncv, mpd = E.getDimensions()
  Print( "Number of requested eigenvalues: %d" % nev )

  tol, maxit = E.getTolerances()
  Print( "Stopping condition: tol=%.4g, maxit=%d" % (tol, maxit) )

For retrieving the solution, it is necessary to find out how many eigenpairs
have converged to the requested precision:

::

  nconv = E.getConverged()
  Print( "Number of converged eigenpairs %d" % nconv )

For each of the ``nconv`` eigenpairs, we can retrieve the eigenvalue ``k``,
and the eigenvector, which is represented by means of two petsc4py vectors
``vr`` and ``vi`` (the real and imaginary part of the eigenvector, since for
real matrices the eigenvalue and eigenvector may be complex). We also compute
the corresponding relative errors in order to make sure that the computed
solution is indeed correct:

::

  if nconv > 0:
    # Create the results vectors
    v, _ = A.createVecs()
    #
    Print()
    Print("        k          ||Ax-kx||/||kx|| ")
    Print("----------------- ------------------")
    for i in range(nconv):
      k = E.getEigenpair(i, v)
      error = E.computeError(i)
      Print( " %12f       %12g" % (k, error) )
    Print()


Example of command-line usage
-----------------------------

Now we illustrate how to specify command-line options in order to extract the
full potential of slepc4py.

A simple execution of the ``demo/ex1.py`` script will result in the following
output:

.. code-block:: console

  $ python demo/ex1.py

  ******************************
  *** SLEPc Solution Results ***
  ******************************

  Number of iterations of the method: 4
  Solution method: krylovschur
  Number of requested eigenvalues: 1
  Stopping condition: tol=1e-07, maxit=100
  Number of converged eigenpairs 4

      k          ||Ax-kx||/||kx||
  ----------------- ------------------
       3.989739        5.76012e-09
       3.959060        1.41957e-08
       3.908279        6.74118e-08
       3.837916        8.34269e-08

For specifying different setting for the solver parameters, we can use SLEPc
command-line options with the -eps prefix. For instance, to change the number
of requested eigenvalues and the tolerance:

.. code-block:: console

  $ python demo/ex1.py -eps_nev 10 -eps_tol 1e-11

The method used by the solver object can also be set at run time:

.. code-block:: console

  $ python demo/ex1.py -eps_type subspace

All the above settings can also be changed within the source code by making
use of the appropriate slepc4py method. Since options can be set from within
the code and the command-line, it is often useful to view the particular
settings that are currently being used:

.. code-block:: console

  $ python demo/ex1.py -eps_view

  EPS Object: 1 MPI process
    type: krylovschur
      50% of basis vectors kept after restart
      using the locking variant
    problem type: symmetric eigenvalue problem
    selected portion of the spectrum: largest eigenvalues in magnitude
    number of eigenvalues (nev): 1
    number of column vectors (ncv): 16
    maximum dimension of projected problem (mpd): 16
    maximum number of iterations: 100
    tolerance: 1e-08
    convergence test: relative to the eigenvalue
  BV Object: 1 MPI process
    type: mat
    17 columns of global length 30
    orthogonalization method: classical Gram-Schmidt
    orthogonalization refinement: if needed (eta: 0.7071)
    block orthogonalization method: GS
    doing matmult as a single matrix-matrix product
  DS Object: 1 MPI process
    type: hep
    solving the problem with: Implicit QR method (_steqr)
  ST Object: 1 MPI process
    type: shift
    shift: 0
    number of matrices: 1

Note that for computing eigenvalues of smallest magnitude we can use the
option ``-eps_smallest_magnitude``, but for interior eigenvalues things are
not so straightforward. One possibility is to try with harmonic extraction,
for instance to get the eigenvalues closest to 0.6:

.. code-block:: console

  $ python demo/ex1.py -eps_harmonic -eps_target 0.6

Depending on the problem, harmonic extraction may fail to converge. In those
cases, it is necessary to specify a spectral transformation other than the
default. In the command-line, this is indicated with the ``-st_`` prefix. For
example, shift-and-invert with a value of the shift equal to 0.6 would be:

.. code-block:: console

  $ python demo/ex1.py -st_type sinvert -eps_target 0.6