.. Description of the jdsym module
.. module:: jdsym
.. _jdsym-page:

=================
Eigenvalue Solver
=================

The :mod:`jdsym` Module
=======================

The ``jdsym`` module provides an implementation of the JDSYM algorithm, that is
conveniently callable from Python. JDSYM is an eigenvalue solver to compute
eigenpairs of a generalised matrix eigenvalue problem of the form

.. math:: \mathbf{A} \mathbf{x} = \lambda \mathbf{M} \mathbf{x}
   :label: eq:python:2

or a standard eigenvalue problem of the form

.. math:: \mathbf{A} \mathbf{x} = \lambda \mathbf{x}
   :label: eq:python:3
  
where :math:`\mathbf{A}` is symmetric and :math:`\mathbf{M}` is symmetric
positive definite.

The module exports a single function:

.. function:: jdsym(A, M, K, kmax, tau, jdtol, itmax, linsolver, **kwargs)

   Implements Jacobi-Davidson iterative method to identify a given number of
   eigenvalues near a target value.

   :parameters:

        :A: the matrix :math:`\mathbf{A}` in :eq:`eq:python:2` or
             :eq:`eq:python:3`. ``A`` must provide the ``shape`` attribute and
             the ``matvec`` and ``matvec_transp`` methods.
        :M: the matrix :math:`\mathbf{M}`
             in :eq:`eq:python:2`. :math:`\mathbf{M}` must provide the ``shape``
             attribute and the ``matvec`` and ``matvec_transp`` methods. If the
             standard eigenvalue problem :eq:`eq:python:3` is to be solved,
             ``M`` should be set to ``None``.
        :K: a preconditioner object that supplies the ``shape`` attribute and
             the ``precon`` method. If no preconditioner is to used, then the
             ``None`` value can be passed for this parameter.
        :kmax: the number of eigenpairs to be computed.
        :tau: the target value :math:`\tau`. Eigenvalues in the vicinity
             of :math:`\tau` will be computed.
        :jdtol: the convergence tolerance for
             eigenpairs :math:`(\lambda,\mathbf{x})`. The converged eigenpairs
             have a residual :math:`\|\mathbf{A} \mathbf{x} - \lambda \mathbf{M}
             \mathbf{x}\|_2` less than ``jdtol``.
        :itmax: an integer that specifies the maximum number of Jacobi-Davidson
             iterations to perform.
        :linsolver: a function that implements an iterative method for solving
             linear systems of equations. The function ``linsolver`` is required
             to conform to the standards mentioned in :ref:`itsolvers-page`.

   :keywords:

        :jmax: the maximum dimension of the search subspace (default: 25).
        :jmin: the dimension of the search subspace after a restart (default:
             10).
        :blksize: the block size used in the JDSYM algorithm (default: 1).
        :blkwise: is an integer that affects the convergence criterion if
             ``blksize`` is larger than 1 (default: 0).
        :V0: a NumPy array of rank one or two. It specifies the initial
             search subspace (default: a randomly generated initial search
             subspace).
        :optype: is an integer specifying the operator type used in the
             correction equation. If ``optype=1``, the non-symmetric version is
             used. If ``optype=2``, the  symmetric version is used (default: 2).
        :linitmax: the maximum number of steps taken in the inner iteration
             (iterative linear solver) (default: 200).
        :eps_tr: the tracking parameter (default: 1.0e-3).
        :toldecay: is a float value that influences the dynamic adaptation of
             the stopping criterion of the inner iteration (default: 1.5).
        :clvl: verbosity level. The higher the ``clvl`` parameter, the more
             output is sent to the standard output. ``clvl=0`` produces no
             output (default: 0).
        :strategy: is an integer specifying shifting and sorting strategy of
             JDSYM. ``strategy=0`` enables the default JDSYM
             algorithm. ``strategy=1`` enables JDSYM to avoid convergence to
             eigenvalues smaller than :math:`\tau` (default: 0).
        :projector: is used to keep the search subspace and the eigenvectors
             in a certain subspace. The parameter ``projector`` can be any
             Python object that has a ``shape`` attribute and a ``project``
             method. The ``project`` method takes a vector (a rank-1 NumPy
             array) as its sole argument and projects that vector in-place. This
             parameter can be used to implement the DIRPROJ and SAUG methods
             (default: no projection).
   :returns:
        :kconv: the number of converged eigenpairs.
        :lambda: a rank-1 NumPy array containing the converged eigenvalues.
        :Q: a rank-2 NumPy array containing the converged eigenvectors. The
           i-th eigenvector is accessed by ``Q[:,i]``.
        :it: an integer indicating the number of Jacobi-Davidson steps
           (outer iteration steps) performed.



Example: Maxwell Problem
------------------------

.. todo:: Update the timings below.

.. warning:: The timings below are Roman's old benchmarks. We should run them
   again.

The following code illustrates the use of the ``jdsym`` module.  Two
matrices :math:`\mathbf{A}` and :math:`\mathbf{M}` are read from files. A Jacobi
preconditioner from :math:`\mathbf{A} - \tau\mathbf{M}` is built. Then the JDSYM
eigensolver is called, calculating 5 eigenvalues near 25.0 and the associated
eigenvalues to an accuracy of :math:`10^{-10}`.  We set ``strategy=1``
to avoid convergence to the high-dimensional null space of
(:math:`\mathbf{A}`, :math:`\mathbf{M}`)::

    from pysparse.sparse import spmatrix
    from pysparse.itsolvers import krylov
    from pysparse.eig import jdsym
    from pysparse.precon import precon

    A = spmatrix.ll_mat_from_mtx('edge6x3x5_A.mtx')
    M = spmatrix.ll_mat_from_mtx('edge6x3x5_B.mtx')
    tau = 25.0

    Atau = A.copy()
    Atau.shift(-tau, M)
    K = precon.jacobi(Atau)

    A = A.to_sss(); M = M.to_sss()
    k_conv, lmbd, Q, it  = jdsym.jdsym(A, M, K, 5, tau,
                                       1e-10, 150, krylov.qmrs,
                                       jmin=5, jmax=10, clvl=1, strategy=1)

This code takes 33.71 seconds to compute the five eigenpairs.  A native
C version, using the same computational kernels, takes 35.64 for the same
task. We expected the Python version to be slower due to the overhead generated
when calling the matrix-vector multiplication and the preconditioner, but
surprisingly the Python code was even a bit faster.
