# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 29.3 `eigen` Chapter.
``eigen`` - Eigenvalues and Eigenvectors
This module provides functions for various types of matrix eigenvalue problem:
standard eigenvalue problems (finding eigenvalues and eigenvectors of a square matrix :math:`A`);
singular value problems (finding singular values and singular vectors of a rectangular matrix :math:`A`);
generalized eigenvalue problems (finding eigenvalues and eigenvectors of a matrix pencil :math:`A-\lambda B`).
quadratic eigenvalue problems (finding eigenvalues and eigenvectors of the quadratic :math:`\lambda^2A+\lambda B+C`).
Functions are provided for both `real` and `complex` data.
The majority of functions for these problems can be found in submodule :mod:`~naginterfaces.library.lapackeig` which contains software derived from LAPACK (see Anderson `et al.` (1999)).
However, you should read `the F02 Introduction <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02intro.html>`__ before turning to submodule :mod:`~naginterfaces.library.lapackeig`, especially if you are a new user. Submodule :mod:`~naginterfaces.library.sparseig` contains functions for large sparse eigenvalue problems, although one such function is also available in this module.
Submodule ``eigen`` and submodule :mod:`~naginterfaces.library.lapackeig` contain **Black Box** (or **Driver**) functions that enable many problems to be solved by a call to a single function, and the decision trees in `Decision Trees <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02intro.html#dtree>`__ direct you to the most appropriate functions in submodule ``eigen`` and submodule :mod:`~naginterfaces.library.lapackeig`.
The submodule ``eigen`` functions call functions in submodule :mod:`~naginterfaces.library.lapacklin` and submodule :mod:`~naginterfaces.library.lapackeig` wherever possible to perform the computations, and there are pointers in `Decision Trees <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02intro.html#dtree>`__ to the relevant decision trees in submodule :mod:`~naginterfaces.library.lapackeig`.
Functionality Index
-------------------
**Black Box functions**
complex eigenproblem
selected eigenvalues and eigenvectors: :meth:`complex_gen_eigsys`
complex quadratic eigenproblem
all eigenvalues and optionally eigenvectors, backward
errors and eigenvalue condition numbers: :meth:`complex_gen_quad`
complex upper triangular matrix
singular values and, optionally, left and/or right singular vectors: :meth:`complex_triang_svd`
generalized real sparse symmetric-definite eigenproblem
selected eigenvalues and eigenvectors: :meth:`real_symm_sparse_eigsys`
real eigenproblem
selected eigenvalues and eigenvectors: :meth:`real_gen_eigsys`
real quadratic eigenproblem
all eigenvalues and optionally eigenvectors, backward
errors and eigenvalue condition numbers: :meth:`real_gen_quad`
real sparse eigenproblem
selected eigenvalues and eigenvectors: :meth:`real_gen_sparse_arnoldi`
real sparse symmetric matrix
driver
selected eigenvalues and eigenvectors: :meth:`real_symm_sparse_arnoldi`
selected eigenvalues and eigenvectors: :meth:`real_symm_sparse_eigsys`
real upper triangular matrix
singular values and, optionally, left and/or right singular vectors: :meth:`real_triang_svd`
**General Purpose functions (see also** submodule :mod:`~naginterfaces.library.sparseig` **)**
real :math:`m\times n` matrix, leading terms SVD: :meth:`real_gen_partialsvd`
For full information please refer to the NAG Library document
https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02intro.html
"""
# NAG Copyright 2017-2023.
[docs]def real_gen_eigsys(crit, a, wl, wu, mest):
r"""
``real_gen_eigsys`` computes selected eigenvalues and eigenvectors of a real general matrix.
.. _f02ec-py2-py-doc:
For full information please refer to the NAG Library document for f02ec
https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02ecf.html
.. _f02ec-py2-py-parameters:
**Parameters**
**crit** : str, length 1
Indicates the criterion for selecting eigenvalues.
:math:`\mathrm{crit} = \texttt{'M'}`
Eigenvalues are selected according to their moduli: :math:`w_l\leq \left\lvert \lambda_i\right\rvert \leq w_u`.
:math:`\mathrm{crit} = \texttt{'R'}`
Eigenvalues are selected according to their real parts: :math:`w_l\leq \mathrm{Re}\left(\lambda_i\right)\leq w_u`.
**a** : float, array-like, shape :math:`\left(n, n\right)`
The :math:`n\times n` general matrix :math:`A`.
**wl** : float
:math:`w_l` and :math:`w_u`, the lower and upper bounds on the criterion for the selected eigenvalues (see :math:`\mathrm{crit}`).
**wu** : float
:math:`w_l` and :math:`w_u`, the lower and upper bounds on the criterion for the selected eigenvalues (see :math:`\mathrm{crit}`).
**mest** : int
:math:`\mathrm{mest}` must be an upper bound on :math:`m`, the number of eigenvalues and eigenvectors selected. No eigenvectors are computed if :math:`\mathrm{mest} < m`.
**Returns**
**a** : float, ndarray, shape :math:`\left(n, n\right)`
Contains the Hessenberg form of the balanced input matrix :math:`A^{\prime }` (see `Further Comments <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02ecf.html#fcomments>`__).
**m** : int
:math:`m`, the number of eigenvalues actually selected.
**wr** : float, ndarray, shape :math:`\left(n\right)`
The first :math:`\mathrm{m}` elements of :math:`\mathrm{wr}` and :math:`\mathrm{wi}` hold the real and imaginary parts, respectively, of the selected eigenvalues; elements :math:`\mathrm{m}+1` to :math:`\textit{n}` contain the other eigenvalues. Complex conjugate pairs of eigenvalues are stored in consecutive elements of the arrays, with the eigenvalue having positive imaginary part first. See also `Further Comments <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02ecf.html#fcomments>`__.
**wi** : float, ndarray, shape :math:`\left(n\right)`
The first :math:`\mathrm{m}` elements of :math:`\mathrm{wr}` and :math:`\mathrm{wi}` hold the real and imaginary parts, respectively, of the selected eigenvalues; elements :math:`\mathrm{m}+1` to :math:`\textit{n}` contain the other eigenvalues. Complex conjugate pairs of eigenvalues are stored in consecutive elements of the arrays, with the eigenvalue having positive imaginary part first. See also `Further Comments <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02ecf.html#fcomments>`__.
**vr** : float, ndarray, shape :math:`\left(n, \mathrm{mest}\right)`
Contains the real parts of the selected eigenvectors, with the :math:`i`\ th column holding the real part of the eigenvector associated with the eigenvalue :math:`\lambda_i` (stored in :math:`\mathrm{wr}[i-1]` and :math:`\mathrm{wi}[i-1]`).
**vi** : float, ndarray, shape :math:`\left(n, \mathrm{mest}\right)`
Contains the imaginary parts of the selected eigenvectors, with the :math:`i`\ th column holding the imaginary part of the eigenvector associated with the eigenvalue :math:`\lambda_i` (stored in :math:`\mathrm{wr}[i-1]` and :math:`\mathrm{wi}[i-1]`).
.. _f02ec-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{mest} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{mest}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{wu} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{wl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{wu} > \mathrm{wl}`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{crit} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{crit} = \texttt{'M'}` or :math:`\texttt{'R'}`.
(`errno` :math:`2`)
The :math:`QR` algorithm failed to converge: only :math:`\langle\mathit{\boldsymbol{value}}\rangle` eigenvalues have been computed; no eigenvectors have been computed.
(`errno` :math:`3`)
There are more than :math:`\mathrm{mest}` eigenvalues in the specified range. :math:`\mathrm{m}` (number of eigenvalues in range) :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{mest} = \langle\mathit{\boldsymbol{value}}\rangle`. No eigenvectors have been computed. Rerun with second dimension of :math:`\mathrm{vr}` and :math:`\mathrm{vi} = \mathrm{mest}\geq \mathrm{m}`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`4`)
Inverse iteration failed to compute all the specified eigenvectors. The number of eigenvectors which failed to converge is :math:`\langle\mathit{\boldsymbol{value}}\rangle`. The corresponding columns of :math:`\mathrm{vr}` and :math:`\mathrm{vi}` are set to zero.
.. _f02ec-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``real_gen_eigsys`` computes selected eigenvalues and the corresponding right eigenvectors of a real general matrix :math:`A`:
.. math::
Ax_i = \lambda_ix_i\text{.}
Eigenvalues :math:`\lambda_i` may be selected either by `modulus`, satisfying:
.. math::
w_l\leq \left\lvert \lambda_i\right\rvert \leq w_u\text{,}
or by `real part`, satisfying:
.. math::
w_l\leq \mathrm{Re}\left(\lambda_i\right)\leq w_u\text{.}
Note that even though :math:`A` is real, :math:`\lambda_i` and :math:`x_i` may be complex.
If :math:`x_i` is an eigenvector corresponding to a complex eigenvalue :math:`\lambda_i`, then the complex conjugate vector :math:`\bar{x}_i` is the eigenvector corresponding to the complex conjugate eigenvalue :math:`\bar{\lambda }_i`.
The eigenvalues in a complex conjugate pair :math:`\lambda_i` and :math:`\bar{\lambda }_i` are either both selected or both not selected.
.. _f02ec-py2-py-references:
**References**
Golub, G H and Van Loan, C F, 1996, `Matrix Computations`, (3rd Edition), Johns Hopkins University Press, Baltimore
"""
raise NotImplementedError
[docs]def real_gen_sparse_arnoldi(a, icolzp, irowix, nev, ncv, sigma, monit=None, option=None, data=None, io_manager=None):
r"""
``real_gen_sparse_arnoldi`` computes selected eigenvalues and eigenvectors of a real sparse general matrix.
Note: this function uses optional algorithmic parameters, see also: :meth:`sparseig.real_init <naginterfaces.library.sparseig.real_init>`, :meth:`sparseig.real_iter <naginterfaces.library.sparseig.real_iter>`, :meth:`sparseig.real_proc <naginterfaces.library.sparseig.real_proc>`, :meth:`sparseig.real_option <naginterfaces.library.sparseig.real_option>`, :meth:`sparseig.real_monit <naginterfaces.library.sparseig.real_monit>`.
.. _f02ek-py2-py-doc:
For full information please refer to the NAG Library document for f02ek
https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02ekf.html
.. _f02ek-py2-py-parameters:
**Parameters**
**a** : float, array-like, shape :math:`\left(\textit{nnz}\right)`
The array of nonzero elements (and diagonal elements if a nonzero inverse shift is to be applied) of the :math:`n\times n` general matrix :math:`A`.
**icolzp** : int, array-like, shape :math:`\left(n+1\right)`
:math:`\mathrm{icolzp}[i-1]` contains the index in :math:`\mathrm{a}` of the start of column :math:`\textit{i}`, for :math:`\textit{i} = 1,2,\ldots,n`; :math:`\mathrm{icolzp}[n]` must contain the value :math:`\textit{nnz}+1`. Thus the number of nonzero elements in column :math:`\textit{i}` of :math:`A` is :math:`\mathrm{icolzp}[i]-\mathrm{icolzp}[i-1]`; when shifts are applied this includes diagonal elements irrespective of value. See `the F11 Introduction <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f11/f11intro.html#background12>`__.
**irowix** : int, array-like, shape :math:`\left(\textit{nnz}\right)`
:math:`\mathrm{irowix}[i-1]` contains the row index for each entry in :math:`\mathrm{a}`. See `the F11 Introduction <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f11/f11intro.html#background12>`__.
**nev** : int
The number of eigenvalues to be computed.
**ncv** : int
The number of Arnoldi basis vectors to use during the computation.
At present there is no `a priori` analysis to guide the selection of :math:`\mathrm{ncv}` relative to :math:`\mathrm{nev}`.
However, it is recommended that :math:`\mathrm{ncv}\geq 2\times \mathrm{nev}+1`.
If many problems of the same type are to be solved, you should experiment with increasing :math:`\mathrm{ncv}` while keeping :math:`\mathrm{nev}` fixed for a given test problem.
This will usually decrease the required number of matrix-vector operations but it also increases the work and storage required to maintain the orthogonal basis vectors.
The optimal 'cross-over' with respect to CPU time is problem dependent and must be determined empirically.
**sigma** : float
If the 'Shifted Inverse Real' mode has been selected then :math:`\mathrm{sigma}` contains the real shift used; otherwise :math:`\mathrm{sigma}` is not referenced. This mode can be selected by setting the appropriate options in the user-supplied function :math:`\mathrm{option}`.
**monit** : None or callable monit(niter, nconv, w, rzest, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.
:math:`\mathrm{monit}` is used to monitor the progress of ``real_gen_sparse_arnoldi``. :math:`\mathrm{monit}` may be **None** if no monitoring is actually required. :math:`\mathrm{monit}` is called after the solution of each eigenvalue sub-problem and also just prior to return from ``real_gen_sparse_arnoldi``.
**Parameters**
**niter** : int
The number of the current Arnoldi iteration.
**nconv** : int
The number of converged eigenvalues so far.
**w** : complex, ndarray, shape :math:`\left(\textit{ncv}\right)`
The first :math:`\mathrm{nconv}` elements of :math:`\mathrm{w}` contain the converged approximate eigenvalues.
**rzest** : float, ndarray, shape :math:`\left(\textit{ncv}\right)`
The first :math:`\mathrm{nconv}` elements of :math:`\mathrm{rzest}` contain the Ritz estimates (error bounds) on the converged approximate eigenvalues.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**option** : None or callable option(comm, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.
You can supply non-default options to the Arnoldi eigensolver by repeated calls to :meth:`sparseig.real_option <naginterfaces.library.sparseig.real_option>` from within :math:`\mathrm{option}`. (Please note that it is only necessary to call :meth:`sparseig.real_option <naginterfaces.library.sparseig.real_option>`; no call to :meth:`sparseig.real_init <naginterfaces.library.sparseig.real_init>` is required from within :math:`\mathrm{option}`.) For example, you can set the mode to 'Shifted Inverse Real', you can increase the 'Iteration Limit' beyond its default and you can print varying levels of detail on the iterative process using 'Print Level'.
If only the default options (including that the eigenvalues of largest magnitude are sought) are to be used then :math:`\mathrm{option}` may be **None**.
**Parameters**
**comm** : dict, communication object, modifiable in place
Communication structure.
This argument has been initialized to hold the default option set.
If you wish to supply non-default options to the Arnoldi eigensolver, this argument may be passed as argument :math:`\textit{comm}` in a call to :meth:`sparseig.real_option <naginterfaces.library.sparseig.real_option>`.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**data** : arbitrary, optional
User-communication data for callback functions.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**a** : float, ndarray, shape :math:`\left(\textit{nnz}\right)`
If a nonzero shifted inverse is to be applied then the diagonal elements of :math:`A` have the shift value, as supplied in :math:`\mathrm{sigma}`, subtracted.
**nconv** : int
The number of converged approximations to the selected eigenvalues. On successful exit, this will normally be either :math:`\mathrm{nev}` or :math:`\mathrm{nev}+1` depending on the number of complex conjugate pairs of eigenvalues returned.
**w** : complex, ndarray, shape :math:`\left(\mathrm{ncv}\right)`
The first :math:`\mathrm{nconv}` elements contain the converged approximations to the selected eigenvalues. Since complex conjugate pairs of eigenvalues appear together, it is possible (given an odd number of converged real eigenvalues) for ``real_gen_sparse_arnoldi`` to return one more eigenvalue than requested.
**v** : float, ndarray, shape :math:`\left(n, \mathrm{ncv}\right)`
Contains the eigenvectors associated with the eigenvalue :math:`\lambda_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{nconv}` (stored in :math:`\mathrm{w}`). For a real eigenvalue, :math:`\lambda_j`, the corresponding eigenvector is real and is stored in :math:`\mathrm{v}[\textit{i}-1,j-1]`, for :math:`\textit{i} = 1,2,\ldots,n`. For complex conjugate pairs of eigenvalues, :math:`w_{{j+1}} = \bar{w_j}`, the real and imaginary parts of the corresponding eigenvectors are stored, respectively, in :math:`\mathrm{v}[\textit{i}-1,j-1]` and :math:`\mathrm{v}[\textit{i}-1,j]`, for :math:`\textit{i} = 1,2,\ldots,n`. The imaginary parts stored are for the first of the conjugate pair of eigenvectors; the other eigenvector in the pair is obtained by negating these imaginary parts.
**resid** : float, ndarray, shape :math:`\left(\mathrm{nev}+1\right)`
The residual :math:`\left\lVert Aw_{\textit{i}}-\lambda_{\textit{i}}w_{\textit{i}}\right\rVert_2` for the estimates to the eigenpair :math:`\lambda_{\textit{i}}` and :math:`w_{\textit{i}}` is returned in :math:`\mathrm{resid}[\textit{i}-1]`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{nconv}`.
.. _f02ek-py2-py-other_params:
**Other Parameters**
**'Advisory'** : int
Default :math:`\text{} = 0`
If the option 'List' is set then option specifications are listed in a 'List' file by setting the option to a file identification (unit) number associated with advisory messages (see :class:`~naginterfaces.base.utils.FileObjManager` and the ``FileObjManager`` method :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`).
**'Defaults'** : valueless
This special keyword may be used to reset all options to their default values.
**'Iteration Limit'** : int
Default :math:`\text{} = 300`
The limit on the number of Arnoldi iterations that can be performed before ``real_gen_sparse_arnoldi`` exits with :math:`\mathrm{errno}` = 24.
**'Largest Magnitude'** : valueless
Default
The Arnoldi iterative method converges on a number of eigenvalues with given properties.
The default is to compute the eigenvalues of largest magnitude using 'Largest Magnitude'.
Alternatively, eigenvalues may be chosen which have 'Largest Real' part, 'Largest Imaginary' part, 'Smallest Magnitude', 'Smallest Real' part or 'Smallest Imaginary' part.
Note that these options select the eigenvalue properties for eigenvalues of :math:`\mathrm{op}` the linear operator determined by the computational mode and problem type.
**'Largest Imaginary'** : valueless
The Arnoldi iterative method converges on a number of eigenvalues with given properties.
The default is to compute the eigenvalues of largest magnitude using 'Largest Magnitude'.
Alternatively, eigenvalues may be chosen which have 'Largest Real' part, 'Largest Imaginary' part, 'Smallest Magnitude', 'Smallest Real' part or 'Smallest Imaginary' part.
Note that these options select the eigenvalue properties for eigenvalues of :math:`\mathrm{op}` the linear operator determined by the computational mode and problem type.
**'Largest Real'** : valueless
The Arnoldi iterative method converges on a number of eigenvalues with given properties.
The default is to compute the eigenvalues of largest magnitude using 'Largest Magnitude'.
Alternatively, eigenvalues may be chosen which have 'Largest Real' part, 'Largest Imaginary' part, 'Smallest Magnitude', 'Smallest Real' part or 'Smallest Imaginary' part.
Note that these options select the eigenvalue properties for eigenvalues of :math:`\mathrm{op}` the linear operator determined by the computational mode and problem type.
**'Smallest Imaginary'** : valueless
The Arnoldi iterative method converges on a number of eigenvalues with given properties.
The default is to compute the eigenvalues of largest magnitude using 'Largest Magnitude'.
Alternatively, eigenvalues may be chosen which have 'Largest Real' part, 'Largest Imaginary' part, 'Smallest Magnitude', 'Smallest Real' part or 'Smallest Imaginary' part.
Note that these options select the eigenvalue properties for eigenvalues of :math:`\mathrm{op}` the linear operator determined by the computational mode and problem type.
**'Smallest Magnitude'** : valueless
The Arnoldi iterative method converges on a number of eigenvalues with given properties.
The default is to compute the eigenvalues of largest magnitude using 'Largest Magnitude'.
Alternatively, eigenvalues may be chosen which have 'Largest Real' part, 'Largest Imaginary' part, 'Smallest Magnitude', 'Smallest Real' part or 'Smallest Imaginary' part.
Note that these options select the eigenvalue properties for eigenvalues of :math:`\mathrm{op}` the linear operator determined by the computational mode and problem type.
**'Smallest Real'** : valueless
The Arnoldi iterative method converges on a number of eigenvalues with given properties.
The default is to compute the eigenvalues of largest magnitude using 'Largest Magnitude'.
Alternatively, eigenvalues may be chosen which have 'Largest Real' part, 'Largest Imaginary' part, 'Smallest Magnitude', 'Smallest Real' part or 'Smallest Imaginary' part.
Note that these options select the eigenvalue properties for eigenvalues of :math:`\mathrm{op}` the linear operator determined by the computational mode and problem type.
**'Nolist'** : valueless
Default
Option 'List' enables printing of each option specification as it is supplied. 'Nolist' suppresses this printing.
**'List'** : valueless
Option 'List' enables printing of each option specification as it is supplied. 'Nolist' suppresses this printing.
**'Monitoring'** : int
Default :math:`\text{} = -1`
If :math:`i > 0`, monitoring information is output to unit number (see :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`) :math:`i` during the solution of each problem; this may be the same as the 'Advisory' unit number (see :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`).
The type of information produced is dependent on the value of 'Print Level', see the description of the option 'Print Level' for details of the information produced.
Please see the ``FileObjManager`` method :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj` to associate a file with a given unit number (see :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`).
**'Print Level'** : int
Default :math:`\text{} = 0`
This controls the amount of printing produced by ``real_gen_sparse_arnoldi`` as follows.
.. rst-class:: nag-rules-none nag-align-left
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`i` |Output |
+===============+==============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================+
|:math:`= 0` |No output except error messages. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`> 0` |The set of selected options. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`= 2` |Problem and timing statistics when all calls to :meth:`sparseig.real_iter <naginterfaces.library.sparseig.real_iter>` have been completed. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 5` |A single line of summary output at each Arnoldi iteration. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 10`|If :math:`\text{‘Monitoring'} > 0`, then at each iteration, the length and additional steps of the current Arnoldi factorization and the number of converged Ritz values; during re-orthogonalization, the norm of initial/restarted starting vector. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 20`|Problem and timing statistics on final exit from :meth:`sparseig.real_iter <naginterfaces.library.sparseig.real_iter>`. If :math:`\text{‘Monitoring'} > 0`, then at each iteration, the number of shifts being applied, the eigenvalues and estimates of the Hessenberg matrix :math:`H`, the size of the Arnoldi basis, the wanted Ritz values and associated Ritz estimates and the shifts applied; vector norms prior to and following re-orthogonalization. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 30`|If :math:`\text{‘Monitoring'} > 0`, then on final iteration, the norm of the residual; when computing the Schur form, the eigenvalues and Ritz estimates both before and after sorting; for each iteration, the norm of residual for compressed factorization and the compressed upper Hessenberg matrix :math:`H`; during re-orthogonalization, the initial/restarted starting vector; during the Arnoldi iteration loop, a restart is flagged and the number of the residual requiring iterative refinement; while applying shifts, the indices of the shifts being applied.|
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 40`|If :math:`\text{‘Monitoring'} > 0`, then during the Arnoldi iteration loop, the Arnoldi vector number and norm of the current residual; while applying shifts, key measures of progress and the order of :math:`H`; while computing eigenvalues of :math:`H`, the last rows of the Schur and eigenvector matrices; when computing implicit shifts, the eigenvalues and Ritz estimates of :math:`H`. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 50`|If :math:`\text{‘Monitoring'} > 0`, then during Arnoldi iteration loop: norms of key components and the active column of :math:`H`, norms of residuals during iterative refinement, the final upper Hessenberg matrix :math:`H`; while applying shifts: number of shifts, shift values, block indices, updated matrix :math:`H`; while computing eigenvalues of :math:`H`: the matrix :math:`H`, the computed eigenvalues and Ritz estimates. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
**'Regular'** : valueless
Default
These options define the computational mode which in turn defines the form of operation :math:`\mathrm{op}\left(x\right)` to be performed.
Given a standard eigenvalue problem in the form :math:`Ax = \lambda x` then the following modes are available with the appropriate operator :math:`\mathrm{op}\left(x\right)`.
.. rst-class:: nag-rules-none nag-align-left
+----------------------+-------------------------------------------------------------------------------+
|'Regular' |:math:`\mathrm{op} = A` |
+----------------------+-------------------------------------------------------------------------------+
|'Shifted Inverse Real'|:math:`\mathrm{op} = \left(A-\sigma I\right)^{-1}` where :math:`\sigma` is real|
+----------------------+-------------------------------------------------------------------------------+
**'Shifted Inverse Real'** : valueless
These options define the computational mode which in turn defines the form of operation :math:`\mathrm{op}\left(x\right)` to be performed.
Given a standard eigenvalue problem in the form :math:`Ax = \lambda x` then the following modes are available with the appropriate operator :math:`\mathrm{op}\left(x\right)`.
.. rst-class:: nag-rules-none nag-align-left
+----------------------+-------------------------------------------------------------------------------+
|'Regular' |:math:`\mathrm{op} = A` |
+----------------------+-------------------------------------------------------------------------------+
|'Shifted Inverse Real'|:math:`\mathrm{op} = \left(A-\sigma I\right)^{-1}` where :math:`\sigma` is real|
+----------------------+-------------------------------------------------------------------------------+
**'Tolerance'** : float
Default :math:`\text{} = \epsilon`
An approximate eigenvalue has deemed to have converged when the corresponding Ritz estimate is within 'Tolerance' relative to the magnitude of the eigenvalue.
.. _f02ek-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n > 0`.
(`errno` :math:`2`)
On entry, :math:`\textit{nnz} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nnz} > 0`.
(`errno` :math:`2`)
On entry, :math:`\textit{nnz} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nnz}\leq n\times n`.
(`errno` :math:`3`)
On entry, in shifted inverse mode, the :math:`j`\ th diagonal element of :math:`A` is not defined, for :math:`j = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`4`)
On entry, :math:`\mathrm{icolzp}[0] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{icolzp}[0] = 1`.
(`errno` :math:`4`)
On entry, for :math:`i = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{icolzp}[i-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{icolzp}[i] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{icolzp}[i-1] < \mathrm{icolzp}[i]`.
(`errno` :math:`4`)
On entry, :math:`\mathrm{icolzp}[n] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{nnz} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{icolzp}[n] = \textit{nnz}+1`.
(`errno` :math:`5`)
On entry, in specification of column :math:`\langle\mathit{\boldsymbol{value}}\rangle`, and for :math:`j = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{irowix}[j-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{irowix}[j] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{irowix}[j-1] < \mathrm{irowix}[j]`.
(`errno` :math:`6`)
On entry, :math:`\mathrm{nev} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{nev} > 0`.
(`errno` :math:`7`)
On entry, :math:`\mathrm{ncv} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{nev} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ncv} > \mathrm{nev}+1`.
(`errno` :math:`7`)
On entry, :math:`\mathrm{ncv} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ncv}\leq n`.
(`errno` :math:`8`)
On entry, the matrix :math:`A-\sigma \times I` is numerically singular and could not be inverted. Try perturbing the value of :math:`\sigma`.
(`errno` :math:`8`)
On entry, the matrix :math:`A-\sigma \times I` is nearly numerically singular and could not be inverted. Try perturbing the value of :math:`\sigma`. Norm of matrix :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle`, Reciprocal condition number :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`21`)
The maximum number of iterations :math:`\text{}\leq 0`, the option 'Iteration Limit' has been set to :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`22`)
An unexpected internal error occurred when solving an eigenproblem.
This error should not occur. Please contact `NAG <https://www.nag.com>`__.
(`errno` :math:`23`)
An unexpected internal error occurred when solving an eigenproblem.
This error should not occur. Please contact `NAG <https://www.nag.com>`__.
(`errno` :math:`24`)
The maximum number of iterations has been reached.
The maximum number of iterations :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle`.
The number of converged eigenvalues :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle`.
See the function document for further details.
(`errno` :math:`25`)
No shifts could be applied during a cycle of the implicitly restarted Arnoldi iteration.
(`errno` :math:`26`)
Could not build an Arnoldi factorization. The size of the current Arnoldi factorization :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`27`)
Error in internal call to compute eigenvalues and corresponding error bounds of the current upper Hessenberg matrix.
Please contact `NAG <https://www.nag.com>`__.
(`errno` :math:`32`)
An unexpected internal error occurred when postprocessing an eigenproblem.
This error should not occur. Please contact `NAG <https://www.nag.com>`__.
(`errno` :math:`33`)
The number of eigenvalues found to sufficient accuracy is zero.
(`errno` :math:`34`)
Internal inconsistency in the number of converged Ritz values. Number counted :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle`, number expected :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`35`)
During calculation of a real Schur form, there was a failure to compute :math:`\langle\mathit{\boldsymbol{value}}\rangle` eigenvalues in a total of :math:`\langle\mathit{\boldsymbol{value}}\rangle` iterations.
(`errno` :math:`36`)
The computed Schur form could not be reordered by an internal call.
This routine returned with :math:`\textit{errno} = \langle\mathit{\boldsymbol{value}}\rangle`.
Please contact `NAG <https://www.nag.com>`__.
(`errno` :math:`37`)
In calculating eigenvectors, an internal call returned with an error.
Please contact `NAG <https://www.nag.com>`__.
**Warns**
**NagCallbackTerminateWarning**
(`errno` :math:`9`)
User requested termination in :math:`\mathrm{monit}`, :math:`\mathrm{istat} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`10`)
User requested termination in :math:`\mathrm{option}`, :math:`\mathrm{istat} = \langle\mathit{\boldsymbol{value}}\rangle`.
.. _f02ek-py2-py-notes:
**Notes**
``real_gen_sparse_arnoldi`` computes selected eigenvalues and the corresponding right eigenvectors of a real sparse general matrix :math:`A`:
.. math::
Aw_i = \lambda_iw_i\text{.}
A specified number, :math:`n_{{ev}}`, of eigenvalues :math:`\lambda_i`, or the shifted inverses :math:`\mu_i = 1/\left(\lambda_i-\sigma \right)`, may be selected either by largest or smallest modulus, largest or smallest real part, or, largest or smallest imaginary part.
Convergence is generally faster when selecting larger eigenvalues, smaller eigenvalues can always be selected by choosing a zero inverse shift (:math:`\sigma = 0.0`).
When eigenvalues closest to a given real value are required then the shifted inverses of largest magnitude should be selected with shift equal to the required real value.
Note that even though :math:`A` is real, :math:`\lambda_i` and :math:`w_i` may be complex.
If :math:`w_i` is an eigenvector corresponding to a complex eigenvalue :math:`\lambda_i`, then the complex conjugate vector :math:`\bar{w}_i` is the eigenvector corresponding to the complex conjugate eigenvalue :math:`\bar{\lambda }_i`.
The eigenvalues in a complex conjugate pair :math:`\lambda_i` and :math:`\bar{\lambda }_i` are either both selected or both not selected.
The sparse matrix :math:`A` is stored in compressed column storage (CCS) format.
See `the F11 Introduction <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f11/f11intro.html#background12>`__.
``real_gen_sparse_arnoldi`` uses an implicitly restarted Arnoldi iterative method to converge approximations to a set of required eigenvalues and corresponding eigenvectors.
Further algorithmic information is given in `Further Comments <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02ekf.html#fcomments>`__ while a fuller discussion is provided in `the F12 Introduction <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f12/f12intro.html>`__.
If shifts are to be performed then operations using shifted inverse matrices are performed using a direct sparse solver; further information on the solver used is provided in `the F11 Introduction <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f11/f11intro.html>`__.
.. _f02ek-py2-py-references:
**References**
Golub, G H and Van Loan, C F, 1996, `Matrix Computations`, (3rd Edition), Johns Hopkins University Press, Baltimore
Lehoucq, R B, Sorensen, D C and Yang, C, 1998, `ARPACK Users' Guide: Solution of Large-scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods`, SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def real_symm_sparse_eigsys(m, noits, tol, dot, image, novecs, x, monit=None, data=None):
r"""
``real_symm_sparse_eigsys`` finds eigenvalues and eigenvectors of a real sparse symmetric or generalized symmetric eigenvalue problem.
.. _f02fj-py2-py-doc:
For full information please refer to the NAG Library document for f02fj
https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02fjf.html
.. _f02fj-py2-py-parameters:
**Parameters**
**m** : int
:math:`m`, the number of eigenvalues required.
**noits** : int
The maximum number of major iterations (eigenvalue sub-problems) to be performed. If :math:`\mathrm{noits}\leq 0`, the value :math:`100` is used in place of :math:`\mathrm{noits}`.
**tol** : float
A relative tolerance to be used in accepting eigenvalues and eigenvectors. If the eigenvalues are required to about :math:`t` significant figures, :math:`\mathrm{tol}` should be set to about :math:`10^{{-t}}`. :math:`d_i` is accepted as an eigenvalue as soon as two successive approximations to :math:`d_i` differ by less than :math:`\left(\left\lvert \tilde{d}_i\right\rvert \times \mathrm{tol}\right)/10`, where :math:`\tilde{d}_i` is the latest approximation to :math:`d_i`. Once an eigenvalue has been accepted, an eigenvector is accepted as soon as :math:`\left(d_if_i\right)/\left(d_i-d_k\right) < \mathrm{tol}`, where :math:`f_i` is the normalized residual of the current approximation to the eigenvector (see `Further Comments <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02fjf.html#fcomments>`__ for further information). The values of the :math:`f_i` and :math:`d_i` can be printed from :math:`\mathrm{monit}`. If :math:`\mathrm{tol}` is supplied outside the range (:math:`\epsilon,1.0`), where :math:`\epsilon` is the machine precision, the value :math:`\epsilon` is used in place of :math:`\mathrm{tol}`.
**dot** : callable retval = dot(z, w, data=None)
:math:`\mathrm{dot}` must return the value :math:`w^\mathrm{T}Bz` for given vectors :math:`w` and :math:`z`.
For the standard eigenvalue problem, where :math:`B = I`, :math:`\mathrm{dot}` must return the dot product :math:`w^\mathrm{T}z`.
**Parameters**
**z** : float, ndarray, shape :math:`\left(n\right)`
The vector :math:`z` for which :math:`w^\mathrm{T}Bz` is required.
**w** : float, ndarray, shape :math:`\left(n\right)`
The vector :math:`w` for which :math:`w^\mathrm{T}Bz` is required.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value :math:`w^\mathrm{T}Bz`.
**image** : callable w = image(z, data=None)
:math:`\mathrm{image}` must return the vector :math:`w = Cz` for a given vector :math:`z`.
**Parameters**
**z** : float, ndarray, shape :math:`\left(n\right)`
The vector :math:`z` for which :math:`Cz` is required.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**w** : float, array-like, shape :math:`\left(n\right)`
The vector :math:`w = Cz`.
**novecs** : int
The number of approximate vectors that are being supplied in :math:`\mathrm{x}`. If :math:`\mathrm{novecs}` is outside the range :math:`\left(0, k\right)`, the value :math:`0` is used in place of :math:`\mathrm{novecs}`.
**x** : float, array-like, shape :math:`\left(n, k\right)`
If :math:`0 < \mathrm{novecs}\leq k`, the first :math:`\mathrm{novecs}` columns of :math:`\mathrm{x}` must contain approximations to the eigenvectors corresponding to the :math:`\mathrm{novecs}` eigenvalues of largest absolute value of :math:`C`. Supplying approximate eigenvectors can be useful when reasonable approximations are known, or when ``real_symm_sparse_eigsys`` is being restarted with a larger value of :math:`\textit{k}`. Otherwise it is not necessary to supply approximate vectors, as simultaneous iteration vectors will be generated randomly by ``real_symm_sparse_eigsys``.
**monit** : None or callable monit(istate, nextit, nevals, nevecs, f, d, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.
:math:`\mathrm{monit}` is used to monitor the progress of ``real_symm_sparse_eigsys``. :math:`\mathrm{monit}` may be **None** if no monitoring is actually required. :math:`\mathrm{monit}` is called after the solution of each eigenvalue sub-problem and also just prior to return from ``real_symm_sparse_eigsys``.
The arguments :math:`\mathrm{istate}` and :math:`\mathrm{nextit}` allow selective printing by :math:`\mathrm{monit}`.
**Parameters**
**istate** : int
Specifies the state of ``real_symm_sparse_eigsys``.
:math:`\mathrm{istate} = 0`
No eigenvalue or eigenvector has just been accepted.
:math:`\mathrm{istate} = 1`
One or more eigenvalues have been accepted since the last call to :math:`\mathrm{monit}`.
:math:`\mathrm{istate} = 2`
One or more eigenvectors have been accepted since the last call to :math:`\mathrm{monit}`.
:math:`\mathrm{istate} = 3`
One or more eigenvalues and eigenvectors have been accepted since the last call to :math:`\mathrm{monit}`.
:math:`\mathrm{istate} = 4`
Return from ``real_symm_sparse_eigsys`` is about to occur.
**nextit** : int
The number of the next iteration.
**nevals** : int
The number of eigenvalues accepted so far.
**nevecs** : int
The number of eigenvectors accepted so far.
**f** : float, ndarray, shape :math:`\left(k\right)`
A vector of error quantities measuring the state of convergence of the simultaneous iteration vectors. See :math:`\mathrm{tol}` and `Further Comments <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02fjf.html#fcomments>`__ for further details. Each element of :math:`\mathrm{f}` is initially set to the value :math:`4.0` and an element remains at :math:`4.0` until the corresponding vector is tested.
**d** : float, ndarray, shape :math:`\left(k\right)`
:math:`\mathrm{d}[i-1]` contains the latest approximation to the absolute value of the :math:`i`\ th eigenvalue of :math:`C`.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**data** : arbitrary, optional
User-communication data for callback functions.
**Returns**
**m** : int
:math:`m^{\prime }`, the number of eigenvalues actually found. It is equal to :math:`m` if no exception or warning is raised on exit, and is less than :math:`m` if :math:`\mathrm{errno}` = 2, 3 or 4. See :ref:`Exceptions <f02fj-py2-py-errors>` and `Further Comments <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02fjf.html#fcomments>`__ for further information.
**noits** : int
The number of iterations actually performed.
**x** : float, ndarray, shape :math:`\left(n, k\right)`
If the function exits successfully or :math:`\mathrm{errno}` = 2, 3 or 4, the first :math:`m^{\prime }` columns contain the eigenvectors corresponding to the eigenvalues returned in the first :math:`m^{\prime }` elements of :math:`\mathrm{d}`; and the next :math:`k-m^{\prime }-1` columns contain approximations to the eigenvectors corresponding to the approximate eigenvalues returned in the next :math:`k-m^{\prime }-1` elements of :math:`\mathrm{d}`. Here :math:`m^{\prime }` is the value returned in :math:`\mathrm{m}`, the number of eigenvalues actually found. The :math:`k`\ th column is used as workspace.
**d** : float, ndarray, shape :math:`\left(k\right)`
If the function exits successfully or :math:`\mathrm{errno}` = 2, 3 or 4, the first :math:`m^{\prime }` elements contain the first :math:`m^{\prime }` eigenvalues in decreasing order of magnitude; and the next :math:`k-m^{\prime }-1` elements contain approximations to the next :math:`k-m^{\prime }-1` eigenvalues. Here :math:`m^{\prime }` is the value returned in :math:`\mathrm{m}`, the number of eigenvalues actually found. :math:`\mathrm{d}[k-1]` contains the value :math:`e` where :math:`\left({-e}, e\right)` is the latest interval over which Chebyshev acceleration is performed.
.. _f02fj-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`k = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`k\leq n`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`k = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m} < k`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 1`.
(`errno` :math:`5`)
Not all requested eigenvalues and vectors have been obtained.
(`errno` :math:`5`)
Convergence of eigenvalue sub-problem occurred.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`2`)
Not all requested eigenvalues and vectors have been obtained.
(`errno` :math:`3`)
Not all requested eigenvalues and vectors have been obtained.
(`errno` :math:`4`)
Not all requested eigenvalues and vectors have been obtained.
**NagCallbackTerminateWarning**
(`errno` :math:`i < 0`)
User set :math:`\mathrm{iflag}` negative in :math:`\mathrm{dot}` or :math:`\mathrm{image}`.
.. _f02fj-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``real_symm_sparse_eigsys`` finds the :math:`m` eigenvalues of largest absolute value and the corresponding eigenvectors for the real eigenvalue problem
.. math::
Cx = \lambda x
where :math:`C` is an :math:`n\times n` matrix such that
.. math::
BC = C^\mathrm{T}B
for a given positive definite matrix :math:`B`. :math:`C` is said to be :math:`B`-symmetric.
Different specifications of :math:`C` allow for the solution of a variety of eigenvalue problems.
For example, when
.. math::
C = A\quad \text{ and }\quad B = I\quad \text{ where }\quad A = A^\mathrm{T}
the function finds the :math:`m` eigenvalues of largest absolute magnitude for the standard symmetric eigenvalue problem
.. math::
Ax = \lambda x\text{.}
The function is intended for the case where :math:`A` is sparse.
As a second example, when
.. math::
C = B^{-1}A
where
.. math::
A = A^\mathrm{T}
the function finds the :math:`m` eigenvalues of largest absolute magnitude for the generalized symmetric eigenvalue problem
.. math::
Ax = \lambda Bx\text{.}
The function is intended for the case where :math:`A` and :math:`B` are sparse.
The function does not require :math:`C` explicitly, but :math:`C` is specified via :math:`\mathrm{image}` which, given an :math:`n`-element vector :math:`z`, computes the image :math:`w` given by
.. math::
w = Cz\text{.}
For instance, in the above example, where :math:`C = B^{-1}A`, :math:`\mathrm{image}` will need to solve the positive definite system of equations :math:`Bw = Az` for :math:`w`.
To find the :math:`m` eigenvalues of smallest absolute magnitude of `(3) <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02fjf.html#eqn3>`__ we can choose :math:`C = A^{-1}` and hence find the reciprocals of the required eigenvalues, so that :math:`\mathrm{image}` will need to solve :math:`Aw = z` for :math:`w`, and correspondingly for `(4) <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02fjf.html#eqn4>`__ we can choose :math:`C = A^{-1}B` and solve :math:`Aw = Bz` for :math:`w`.
A table of examples of choice of :math:`\mathrm{image}` is given in Table [label omitted].
It should be remembered that the function also returns the corresponding eigenvectors and that :math:`B` is positive definite.
Throughout :math:`A` is assumed to be symmetric and, where necessary, nonsingularity is also assumed.
[table omitted]
The matrix :math:`B` also need not be supplied explicitly, but is specified via :math:`\mathrm{dot}` which, given :math:`n`-element vectors :math:`z` and :math:`w`, computes the generalized dot product :math:`w^\mathrm{T}Bz`.
``real_symm_sparse_eigsys`` is based upon routine SIMITZ (see Nikolai (1979)), which is itself a derivative of the Algol procedure ritzit (see Rutishauser (1970)), and uses the method of simultaneous (subspace) iteration. (See Parlett (1998) for a description, analysis and advice on the use of the method.)
The function performs simultaneous iteration on :math:`k > m` vectors.
Initial estimates to :math:`p\leq k` eigenvectors, corresponding to the :math:`p` eigenvalues of :math:`C` of largest absolute value, may be supplied to ``real_symm_sparse_eigsys``.
When possible :math:`k` should be chosen so that the :math:`k`\ th eigenvalue is not too close to the :math:`m` required eigenvalues, but if :math:`k` is initially chosen too small then ``real_symm_sparse_eigsys`` may be re-entered, supplying approximations to the :math:`k` eigenvectors found so far and with :math:`k` then increased.
At each major iteration ``real_symm_sparse_eigsys`` solves an :math:`r\times r` (:math:`r\leq k`) eigenvalue sub-problem in order to obtain an approximation to the eigenvalues for which convergence has not yet occurred.
This approximation is refined by Chebyshev acceleration.
.. _f02fj-py2-py-references:
**References**
Nikolai, P J, 1979, `Algorithm 538: Eigenvectors and eigenvalues of real generalized symmetric matrices by simultaneous iteration`, ACM Trans. Math. Software (5), 118--125
Parlett, B N, 1998, `The Symmetric Eigenvalue Problem`, SIAM, Philadelphia
Rutishauser, H, 1969, `Computational aspects of F L Bauer's simultaneous iteration method`, Numer. Math. (13), 4--13
Rutishauser, H, 1970, `Simultaneous iteration method for symmetric matrices`, Numer. Math. (16), 205--223
"""
raise NotImplementedError
[docs]def real_symm_sparse_arnoldi(n, a, irow, icol, nev, ncv, sigma, monit=None, option=None, data=None, io_manager=None):
r"""
``real_symm_sparse_arnoldi`` computes selected eigenvalues and eigenvectors of a real sparse symmetric matrix.
Note: this function uses optional algorithmic parameters, see also: :meth:`sparseig.real_symm_init <naginterfaces.library.sparseig.real_symm_init>`, :meth:`sparseig.real_symm_iter <naginterfaces.library.sparseig.real_symm_iter>`, :meth:`sparseig.real_symm_proc <naginterfaces.library.sparseig.real_symm_proc>`, :meth:`sparseig.real_symm_option <naginterfaces.library.sparseig.real_symm_option>`, :meth:`sparseig.real_symm_monit <naginterfaces.library.sparseig.real_symm_monit>`.
.. _f02fk-py2-py-doc:
For full information please refer to the NAG Library document for f02fk
https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02fkf.html
.. _f02fk-py2-py-parameters:
**Parameters**
**n** : int
:math:`n`, the order of the matrix :math:`A`.
**a** : float, array-like, shape :math:`\left(\textit{nnz}\right)`
The array of nonzero elements of the lower triangular part of the :math:`n\times n` symmetric matrix :math:`A`.
**irow** : int, array-like, shape :math:`\left(\textit{nnz}\right)`
The row and column indices of the elements supplied in array :math:`\mathrm{a}`.
If :math:`\mathrm{irow}[k-1] = i` and :math:`\mathrm{icol}[k-1] = j` then :math:`A_{{ij}}` is stored in :math:`\mathrm{a}[k-1]`. :math:`\mathrm{irow}` does not need to be ordered: an internal sort will force the correct ordering.
**icol** : int, array-like, shape :math:`\left(\textit{nnz}\right)`
The row and column indices of the elements supplied in array :math:`\mathrm{a}`.
If :math:`\mathrm{irow}[k-1] = i` and :math:`\mathrm{icol}[k-1] = j` then :math:`A_{{ij}}` is stored in :math:`\mathrm{a}[k-1]`. :math:`\mathrm{irow}` does not need to be ordered: an internal sort will force the correct ordering.
**nev** : int
The number of eigenvalues to be computed.
**ncv** : int
The number of Arnoldi basis vectors to use during the computation.
At present there is no `a priori` analysis to guide the selection of :math:`\mathrm{ncv}` relative to :math:`\mathrm{nev}`.
However, it is recommended that :math:`\mathrm{ncv}\geq 2\times \mathrm{nev}+1`.
If many problems of the same type are to be solved, you should experiment with increasing :math:`\mathrm{ncv}` while keeping :math:`\mathrm{nev}` fixed for a given test problem.
This will usually decrease the required number of matrix-vector operations but it also increases the work and storage required to maintain the orthogonal basis vectors.
The optimal 'cross-over' with respect to computation time is problem dependent and must be determined empirically.
**sigma** : float
If the 'Shifted Inverse' mode has been selected then :math:`\mathrm{sigma}` contains the real shift used; otherwise :math:`\mathrm{sigma}` is not referenced. This mode can be selected by setting the appropriate options in the user-supplied function :math:`\mathrm{option}`.
**monit** : None or callable monit(niter, nconv, w, rzest, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.
:math:`\mathrm{monit}` is used to monitor the progress of ``real_symm_sparse_arnoldi``. :math:`\mathrm{monit}` may be **None** if no monitoring is actually required. :math:`\mathrm{monit}` is called after the solution of each eigenvalue sub-problem and also just prior to return from ``real_symm_sparse_arnoldi``.
**Parameters**
**niter** : int
The number of the current Arnoldi iteration.
**nconv** : int
The number of converged eigenvalues so far.
**w** : float, ndarray, shape :math:`\left(\textit{ncv}\right)`
The first :math:`\mathrm{nconv}` elements of :math:`\mathrm{w}` contain the converged approximate eigenvalues.
**rzest** : float, ndarray, shape :math:`\left(\textit{ncv}\right)`
The first :math:`\mathrm{nconv}` elements of :math:`\mathrm{rzest}` contain the Ritz estimates (error bounds) on the converged approximate eigenvalues.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**option** : None or callable option(comm, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.
You can supply non-default options to the Arnoldi eigensolver by repeated calls to :meth:`sparseig.real_symm_option <naginterfaces.library.sparseig.real_symm_option>` from within :math:`\mathrm{option}`. (Please note that it is only necessary to call :meth:`sparseig.real_symm_option <naginterfaces.library.sparseig.real_symm_option>`; no call to :meth:`sparseig.real_symm_init <naginterfaces.library.sparseig.real_symm_init>` is required from within :math:`\mathrm{option}`.) For example, you can set the mode to 'Shifted Inverse', you can increase the 'Iteration Limit' beyond its default and you can print varying levels of detail on the iterative process using 'Print Level'.
If only the default options (including that the eigenvalues of largest magnitude are sought) are to be used then :math:`\mathrm{option}` may be **None**.
**Parameters**
**comm** : dict, communication object, modifiable in place
Communication structure.
This argument has been initialized to hold the default option set.
If you wish to supply non-default options to the Arnoldi eigensolver, this argument may be passed as argument :math:`\textit{comm}` in a call to :meth:`sparseig.real_symm_option <naginterfaces.library.sparseig.real_symm_option>`.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**data** : arbitrary, optional
User-communication data for callback functions.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**nconv** : int
The number of converged approximations to the selected eigenvalues. On successful exit, this will normally be :math:`\mathrm{nev}`.
**w** : float, ndarray, shape :math:`\left(\mathrm{ncv}\right)`
The first :math:`\mathrm{nconv}` elements contain the converged approximations to the selected eigenvalues.
**v** : float, ndarray, shape :math:`\left(\mathrm{n}, \mathrm{ncv}\right)`
Contains the eigenvectors associated with the eigenvalue :math:`\lambda_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{nconv}` (stored in :math:`\mathrm{w}`). For eigenvalue, :math:`\lambda_j`, the corresponding eigenvector is stored in :math:`\mathrm{v}[\textit{i}-1,j-1]`, for :math:`\textit{i} = 1,2,\ldots,n`.
**resid** : float, ndarray, shape :math:`\left(\mathrm{nev}\right)`
The residual :math:`\left\lVert Aw_{\textit{i}}-\lambda_{\textit{i}}w_{\textit{i}}\right\rVert_2` for the estimates to the eigenpair :math:`\lambda_{\textit{i}}` and :math:`w_{\textit{i}}` is returned in :math:`\mathrm{resid}[\textit{i}-1]`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{nconv}`.
.. _f02fk-py2-py-other_params:
**Other Parameters**
**'Advisory'** : int
Default :math:`= \text{advisory message unit number}`
If the option 'List' is set then option specifications are listed in a 'List' file by setting the option to a file identification (unit) number associated with advisory messages (see :class:`~naginterfaces.base.utils.FileObjManager` and the ``FileObjManager`` method :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`).
**'Defaults'** : valueless
This special keyword may be used to reset all options to their default values.
**'Iteration Limit'** : int
Default :math:`\text{} = 300`
The limit on the number of Lanczos iterations that can be performed before :meth:`sparseig.real_symm_iter <naginterfaces.library.sparseig.real_symm_iter>` exits.
If not all requested eigenvalues have converged to within 'Tolerance' and the number of Lanczos iterations has reached this limit then :meth:`sparseig.real_symm_iter <naginterfaces.library.sparseig.real_symm_iter>` exits with an error; :meth:`sparseig.real_symm_proc <naginterfaces.library.sparseig.real_symm_proc>` can still be called subsequently to return the number of converged eigenvalues, the converged eigenvalues and, if requested, the corresponding eigenvectors.
**'Largest Magnitude'** : valueless
Default
The Lanczos iterative method converges on a number of eigenvalues with given properties.
The default is for :meth:`sparseig.real_symm_iter <naginterfaces.library.sparseig.real_symm_iter>` to compute the eigenvalues of largest magnitude using 'Largest Magnitude'.
Alternatively, eigenvalues may be chosen which have 'Largest Algebraic' part, 'Smallest Magnitude', or 'Smallest Algebraic' part; or eigenvalues which are from 'Both Ends' of the algebraic spectrum.
**'Both Ends'** : valueless
The Lanczos iterative method converges on a number of eigenvalues with given properties.
The default is for :meth:`sparseig.real_symm_iter <naginterfaces.library.sparseig.real_symm_iter>` to compute the eigenvalues of largest magnitude using 'Largest Magnitude'.
Alternatively, eigenvalues may be chosen which have 'Largest Algebraic' part, 'Smallest Magnitude', or 'Smallest Algebraic' part; or eigenvalues which are from 'Both Ends' of the algebraic spectrum.
**'Largest Algebraic'** : valueless
The Lanczos iterative method converges on a number of eigenvalues with given properties.
The default is for :meth:`sparseig.real_symm_iter <naginterfaces.library.sparseig.real_symm_iter>` to compute the eigenvalues of largest magnitude using 'Largest Magnitude'.
Alternatively, eigenvalues may be chosen which have 'Largest Algebraic' part, 'Smallest Magnitude', or 'Smallest Algebraic' part; or eigenvalues which are from 'Both Ends' of the algebraic spectrum.
**'Smallest Algebraic'** : valueless
The Lanczos iterative method converges on a number of eigenvalues with given properties.
The default is for :meth:`sparseig.real_symm_iter <naginterfaces.library.sparseig.real_symm_iter>` to compute the eigenvalues of largest magnitude using 'Largest Magnitude'.
Alternatively, eigenvalues may be chosen which have 'Largest Algebraic' part, 'Smallest Magnitude', or 'Smallest Algebraic' part; or eigenvalues which are from 'Both Ends' of the algebraic spectrum.
**'Smallest Magnitude'** : valueless
The Lanczos iterative method converges on a number of eigenvalues with given properties.
The default is for :meth:`sparseig.real_symm_iter <naginterfaces.library.sparseig.real_symm_iter>` to compute the eigenvalues of largest magnitude using 'Largest Magnitude'.
Alternatively, eigenvalues may be chosen which have 'Largest Algebraic' part, 'Smallest Magnitude', or 'Smallest Algebraic' part; or eigenvalues which are from 'Both Ends' of the algebraic spectrum.
**'Nolist'** : valueless
Default
Option 'List' enables printing of each option specification as it is supplied. 'Nolist' suppresses this printing.
**'List'** : valueless
Option 'List' enables printing of each option specification as it is supplied. 'Nolist' suppresses this printing.
**'Monitoring'** : int
Default :math:`\text{} = -1`
If :math:`i > 0`, monitoring information is output to unit number (see :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`) :math:`i` during the solution of each problem; this may be the same as the 'Advisory' unit number (see :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`).
The type of information produced is dependent on the value of 'Print Level', see the description of the option 'Print Level' for details of the information produced.
Please see the ``FileObjManager`` method :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj` to associate a file with a given unit number (see :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`).
**'Print Level'** : int
Default :math:`\text{} = 0`
This controls the amount of printing produced by ``real_symm_sparse_arnoldi`` as follows.
.. rst-class:: nag-rules-none nag-align-left
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`i` |Output |
+===============+==================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================+
|:math:`= 0` |No output except error messages. If you want to suppress all output, set :math:`\text{‘Print Level'} = 0`. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`> 0` |The set of selected options. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`= 2` |Problem and timing statistics on final exit from :meth:`sparseig.real_symm_iter <naginterfaces.library.sparseig.real_symm_iter>`. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 5` |A single line of summary output at each Lanczos iteration. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 10`|If 'Monitoring' is set, then at each iteration, the length and additional steps of the current Lanczos factorization and the number of converged Ritz values; during re-orthogonalization, the norm of initial/restarted starting vector; on a final Lanczos iteration, the number of update iterations taken, the number of converged eigenvalues, the converged eigenvalues and their Ritz estimates. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 20`|Problem and timing statistics on final exit from :meth:`sparseig.real_symm_iter <naginterfaces.library.sparseig.real_symm_iter>`. If :math:`\text{‘Monitoring'} > 0`, 'Monitoring' is set, then at each iteration, the number of shifts being applied, the eigenvalues and estimates of the symmetric tridiagonal matrix :math:`H`, the size of the Lanczos basis, the wanted Ritz values and associated Ritz estimates and the shifts applied; vector norms prior to and following re-orthogonalization. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 30`|If :math:`\text{‘Monitoring'} > 0`, 'Monitoring' is set, then on final iteration, the norm of the residual; when computing the Schur form, the eigenvalues and Ritz estimates both before and after sorting; for each iteration, the norm of residual for compressed factorization and the symmetric tridiagonal matrix :math:`H`; during re-orthogonalization, the initial/restarted starting vector; during the Lanczos iteration loop, a restart is flagged and the number of the residual requiring iterative refinement; while applying shifts, some indices.|
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 40`|If :math:`\text{‘Monitoring'} > 0`, 'Monitoring' is set, then during the Lanczos iteration loop, the Lanczos vector number and norm of the current residual; while applying shifts, key measures of progress and the order of :math:`H`; while computing eigenvalues of :math:`H`, the last rows of the Schur and eigenvector matrices; when computing implicit shifts, the eigenvalues and Ritz estimates of :math:`H`. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\geq 50`|If 'Monitoring' is set, then during Lanczos iteration loop: norms of key components and the active column of :math:`H`, norms of residuals during iterative refinement, the final symmetric tridiagonal matrix :math:`H`; while applying shifts: number of shifts, shift values, block indices, updated tridiagonal matrix :math:`H`; while computing eigenvalues of :math:`H`: the diagonals of :math:`H`, the computed eigenvalues and Ritz estimates. |
+---------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
Note that setting :math:`\text{‘Print Level'} \geq 30` can result in very lengthy 'Monitoring' output.
**'Regular'** : valueless
Default
These options define the computational mode which in turn defines the form of operation :math:`\mathrm{op}\left(x\right)` to be performed.
.. rst-class:: nag-rules-none nag-align-left
+-----------------+-------------------------------------------------------------------------------+
|'Regular' |:math:`\mathrm{op} = A` |
+-----------------+-------------------------------------------------------------------------------+
|'Shifted Inverse'|:math:`\mathrm{op} = \left(A-\sigma I\right)^{-1}` where :math:`\sigma` is real|
+-----------------+-------------------------------------------------------------------------------+
|'Regular Inverse'|:math:`\mathrm{op} = A^{-1}` |
+-----------------+-------------------------------------------------------------------------------+
**'Regular Inverse'** : valueless
These options define the computational mode which in turn defines the form of operation :math:`\mathrm{op}\left(x\right)` to be performed.
.. rst-class:: nag-rules-none nag-align-left
+-----------------+-------------------------------------------------------------------------------+
|'Regular' |:math:`\mathrm{op} = A` |
+-----------------+-------------------------------------------------------------------------------+
|'Shifted Inverse'|:math:`\mathrm{op} = \left(A-\sigma I\right)^{-1}` where :math:`\sigma` is real|
+-----------------+-------------------------------------------------------------------------------+
|'Regular Inverse'|:math:`\mathrm{op} = A^{-1}` |
+-----------------+-------------------------------------------------------------------------------+
**'Shifted Inverse'** : valueless
These options define the computational mode which in turn defines the form of operation :math:`\mathrm{op}\left(x\right)` to be performed.
.. rst-class:: nag-rules-none nag-align-left
+-----------------+-------------------------------------------------------------------------------+
|'Regular' |:math:`\mathrm{op} = A` |
+-----------------+-------------------------------------------------------------------------------+
|'Shifted Inverse'|:math:`\mathrm{op} = \left(A-\sigma I\right)^{-1}` where :math:`\sigma` is real|
+-----------------+-------------------------------------------------------------------------------+
|'Regular Inverse'|:math:`\mathrm{op} = A^{-1}` |
+-----------------+-------------------------------------------------------------------------------+
**'Tolerance'** : float
Default :math:`\text{} = \epsilon`
An approximate eigenvalue has deemed to have converged when the corresponding Ritz estimate is within 'Tolerance' relative to the magnitude of the eigenvalue.
.. _f02fk-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n} > 0`.
(`errno` :math:`2`)
On entry, :math:`\textit{nnz} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nnz} > 0`.
(`errno` :math:`2`)
On entry, :math:`\textit{nnz} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nnz}\leq \mathrm{n}\times \left(\mathrm{n}+1\right)/2`.
(`errno` :math:`4`)
On entry, for :math:`i = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{irow}[i-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq \mathrm{irow}[i-1]\leq \mathrm{n}`.
(`errno` :math:`5`)
On entry, for :math:`i = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{icol}[i-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{irow}[i-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq \mathrm{icol}[i-1]\leq \mathrm{irow}[i-1]`.
(`errno` :math:`6`)
On entry, :math:`\mathrm{nev} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{nev} > 0`.
(`errno` :math:`6`)
On entry, :math:`\mathrm{nev} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{nev} < \left(\mathrm{n}-1\right)`.
(`errno` :math:`7`)
On entry, :math:`\mathrm{ncv} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{nev} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ncv} > \mathrm{nev}`.
(`errno` :math:`7`)
On entry, :math:`\mathrm{ncv} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ncv}\leq \mathrm{n}`.
(`errno` :math:`8`)
On entry, the matrix :math:`\left(A-\sigma I\right)` is numerically singular and could not be inverted. Try perturbing the value of :math:`\sigma`.
(`errno` :math:`20`)
The maximum number of iterations, through the option 'Iteration Limit', has been set to a non-positive value.
(`errno` :math:`21`)
The option 'Both Ends' has been set but only :math:`1` eigenvalue is requested.
(`errno` :math:`22`)
The maximum number of iterations has been reached.
The maximum number of iterations :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle`.
The number of converged eigenvalues :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle`.
See the function document for further details.
(`errno` :math:`30`)
A serious error, code :math:`\left(\langle\mathit{\boldsymbol{value}}\rangle, \langle\mathit{\boldsymbol{value}}\rangle\right)`, has occurred in an internal call. Check all function calls and array sizes. If the call is correct then please contact `NAG <https://www.nag.com>`__ for assistance.
**Warns**
**NagCallbackTerminateWarning**
(`errno` :math:`9`)
User requested termination in :math:`\mathrm{monit}`, :math:`\mathrm{istat} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`10`)
User requested termination in :math:`\mathrm{option}`, :math:`\mathrm{istat} = \langle\mathit{\boldsymbol{value}}\rangle`.
.. _f02fk-py2-py-notes:
**Notes**
``real_symm_sparse_arnoldi`` computes selected eigenvalues and the corresponding right eigenvectors of a real sparse symmetric matrix :math:`A`:
.. math::
Av_i = \lambda_iv_i\text{.}
A specified number, :math:`n_{{ev}}`, of eigenvalues :math:`\lambda_i`, or the shifted inverses :math:`\mu_i = 1/\left(\lambda_i-\sigma \right)`, may be selected either by largest or smallest modulus, largest or smallest value, or, largest and smallest values (both ends).
Convergence is generally faster when selecting larger eigenvalues, smaller eigenvalues can always be selected by choosing a zero inverse shift (:math:`\sigma = 0.0`).
When eigenvalues closest to a given value are required then the shifted inverses of largest magnitude should be selected with shift equal to the required value.
The sparse matrix :math:`A` is stored in symmetric coordinate storage (SCS) format.
See `the F11 Introduction <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f11/f11intro.html#background11>`__.
``real_symm_sparse_arnoldi`` uses an implicitly restarted Arnoldi (Lanczos) iterative method to converge approximations to a set of required eigenvalues and corresponding eigenvectors.
Further algorithmic information is given in `Further Comments <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02fkf.html#fcomments>`__ while a fuller discussion is provided in `the F12 Introduction <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f12/f12intro.html>`__.
If shifts are to be performed then operations using shifted inverse matrices are performed using a direct sparse solver.
.. _f02fk-py2-py-references:
**References**
Golub, G H and Van Loan, C F, 1996, `Matrix Computations`, (3rd Edition), Johns Hopkins University Press, Baltimore
HSL, 2011, `A collection of Fortran codes for large-scale scientific computation`, http://www.hsl.rl.ac.uk/
Lehoucq, R B, Sorensen, D C and Yang, C, 1998, `ARPACK Users' Guide: Solution of Large-scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods`, SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_gen_eigsys(crit, n, a, wl, wu, mest):
r"""
``complex_gen_eigsys`` computes selected eigenvalues and eigenvectors of a complex general matrix.
.. _f02gc-py2-py-doc:
For full information please refer to the NAG Library document for f02gc
https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02gcf.html
.. _f02gc-py2-py-parameters:
**Parameters**
**crit** : str, length 1
Indicates the criterion for selecting eigenvalues.
:math:`\mathrm{crit} = \texttt{'M'}`
Eigenvalues are selected according to their moduli: :math:`w_l\leq \left\lvert \lambda_i\right\rvert \leq w_u`.
:math:`\mathrm{crit} = \texttt{'R'}`
Eigenvalues are selected according to their real parts: :math:`w_l\leq \mathrm{Re}\left(\lambda_i\right)\leq w_u`.
**n** : int
:math:`n`, the order of the matrix :math:`A`.
**a** : complex, array-like, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
The :math:`n\times n` general matrix :math:`A`.
**wl** : float
:math:`w_l` and :math:`w_u`, the lower and upper bounds on the criterion for the selected eigenvalues (see :math:`\mathrm{crit}`).
**wu** : float
:math:`w_l` and :math:`w_u`, the lower and upper bounds on the criterion for the selected eigenvalues (see :math:`\mathrm{crit}`).
**mest** : int
:math:`\mathrm{mest}` must be an upper bound on :math:`m`, the number of eigenvalues and eigenvectors selected. No eigenvectors are computed if :math:`\mathrm{mest} < m`.
**Returns**
**a** : complex, ndarray, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
Contains the Hessenberg form of the balanced input matrix :math:`A^{\prime }` (see `Further Comments <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02gcf.html#fcomments>`__).
**m** : int
:math:`m`, the number of eigenvalues actually selected.
**w** : complex, ndarray, shape :math:`\left(\mathrm{n}\right)`
The first :math:`\mathrm{m}` elements of :math:`\mathrm{w}` hold the selected eigenvalues; elements :math:`\mathrm{m}+1` to :math:`\mathrm{n}` contain the other eigenvalues.
**v** : complex, ndarray, shape :math:`\left(\mathrm{n}, \mathrm{mest}\right)`
Contains the selected eigenvectors, with the :math:`i`\ th column holding the eigenvector associated with the eigenvalue :math:`\lambda_i` (stored in :math:`\mathrm{w}[i-1]`).
.. _f02gc-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{mest} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{mest}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{wu} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{wl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{wu} > \mathrm{wl}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{crit} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{crit} = \texttt{'M'}` or :math:`\texttt{'R'}`.
(`errno` :math:`2`)
The :math:`QR` algorithm failed to converge: only :math:`\langle\mathit{\boldsymbol{value}}\rangle` eigenvalues have been computed; no eigenvectors have been computed.
(`errno` :math:`3`)
There are more than :math:`\mathrm{mest}` eigenvalues in the specified range. :math:`\mathrm{m}` (number of eigenvalues in range) :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{mest} = \langle\mathit{\boldsymbol{value}}\rangle`. No eigenvectors have been computed. Rerun with second dimension of :math:`\mathrm{v} = \mathrm{mest}\geq \mathrm{m}`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`4`)
Inverse iteration failed to compute all the specified eigenvectors. The number of eigenvectors which failed to converge is :math:`\langle\mathit{\boldsymbol{value}}\rangle`. The corresponding columns of :math:`\mathrm{v}` are set to zero.
.. _f02gc-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``complex_gen_eigsys`` computes selected eigenvalues and the corresponding right eigenvectors of a complex general matrix :math:`A`:
.. math::
Ax_i = \lambda_ix_i\text{.}
Eigenvalues :math:`\lambda_i` may be selected either by `modulus`, satisfying
.. math::
w_l\leq \left\lvert \lambda_i\right\rvert \leq w_u\text{,}
or by `real part`, satisfying
.. math::
w_l\leq \mathrm{Re}\left(\lambda_i\right)\leq w_u\text{.}
.. _f02gc-py2-py-references:
**References**
Golub, G H and Van Loan, C F, 1996, `Matrix Computations`, (3rd Edition), Johns Hopkins University Press, Baltimore
"""
raise NotImplementedError
[docs]def real_gen_quad(scal, jobvl, jobvr, sense, tol, a, b, c):
r"""
``real_gen_quad`` solves the quadratic eigenvalue problem
.. math::
\left(\lambda^2A+\lambda B+C\right)x = 0\text{,}
where :math:`A`, :math:`B` and :math:`C` are real :math:`n\times n` matrices.
The function returns the :math:`2n` eigenvalues, :math:`\lambda_{\textit{j}}`, for :math:`\textit{j} = 1,2,\ldots,2n`, and can optionally return the corresponding right eigenvectors, :math:`x_j` and/or left eigenvectors, :math:`y_j` as well as estimates of the condition numbers of the computed eigenvalues and backward errors of the computed right and left eigenvectors.
A left eigenvector satisfies the equation
.. math::
y^\mathrm{H}\left(\lambda^2A+\lambda B+C\right) = 0\text{,}
where :math:`y^\mathrm{H}` is the complex conjugate transpose of :math:`y`.
:math:`\lambda` is represented as the pair :math:`\left(\alpha, \beta \right)`, such that :math:`\lambda = \alpha /\beta`.
Note that the computation of :math:`\alpha /\beta` may overflow and indeed :math:`\beta` may be zero.
.. _f02jc-py2-py-doc:
For full information please refer to the NAG Library document for f02jc
https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02jcf.html
.. _f02jc-py2-py-parameters:
**Parameters**
**scal** : int
Determines the form of scaling to be performed on :math:`A`, :math:`B` and :math:`C`.
:math:`\mathrm{scal} = 0`
No scaling.
:math:`\mathrm{scal} = 1` (the recommended value)
Fan, Lin and Van Dooren scaling if :math:`\frac{\left\lVert B\right\rVert }{\sqrt{\left\lVert A\right\rVert \times \left\lVert C\right\rVert }} < 10` and no scaling otherwise where :math:`\left\lVert Z\right\rVert` is the Frobenius norm of :math:`Z`.
:math:`\mathrm{scal} = 2`
Fan, Lin and Van Dooren scaling.
:math:`\mathrm{scal} = 3`
Tropical scaling with largest root.
:math:`\mathrm{scal} = 4`
Tropical scaling with smallest root.
**jobvl** : str, length 1
If :math:`\mathrm{jobvl} = \texttt{'N'}`, do not compute left eigenvectors.
If :math:`\mathrm{jobvl} = \texttt{'V'}`, compute the left eigenvectors.
If :math:`\mathrm{sense} = 1`, :math:`2`, :math:`4`, :math:`5`, :math:`6` or :math:`7`, :math:`\mathrm{jobvl}` must be set to 'V'.
**jobvr** : str, length 1
If :math:`\mathrm{jobvr} = \texttt{'N'}`, do not compute right eigenvectors.
If :math:`\mathrm{jobvr} = \texttt{'V'}`, compute the right eigenvectors.
If :math:`\mathrm{sense} = 1`, :math:`3`, :math:`4`, :math:`5`, :math:`6` or :math:`7`, :math:`\mathrm{jobvr}` must be set to 'V'.
**sense** : int
Determines whether, or not, condition numbers and backward errors are computed.
:math:`\mathrm{sense} = 0`
Do not compute condition numbers, or backward errors.
:math:`\mathrm{sense} = 1`
Just compute condition numbers for the eigenvalues.
:math:`\mathrm{sense} = 2`
Just compute backward errors for the left eigenpairs.
:math:`\mathrm{sense} = 3`
Just compute backward errors for the right eigenpairs.
:math:`\mathrm{sense} = 4`
Compute backward errors for the left and right eigenpairs.
:math:`\mathrm{sense} = 5`
Compute condition numbers for the eigenvalues and backward errors for the left eigenpairs.
:math:`\mathrm{sense} = 6`
Compute condition numbers for the eigenvalues and backward errors for the right eigenpairs.
:math:`\mathrm{sense} = 7`
Compute condition numbers for the eigenvalues and backward errors for the left and right eigenpairs.
**tol** : float
:math:`\mathrm{tol}` is used as the tolerance for making decisions on rank in the deflation procedure. If :math:`\mathrm{tol}` is zero on entry then :math:`n\times \text{machine precision}` is used in place of :math:`\mathrm{tol}`, where machine precision is as returned by function :meth:`machine.precision <naginterfaces.library.machine.precision>`. A diagonal element of a triangular matrix, :math:`R`, is regarded as zero if :math:`\left\lvert r_{{jj}}\right\rvert \leq \mathrm{tol}\times \mathrm{size}\left(X\right)`, or :math:`n\times \text{machine precision}\times \mathrm{size}\left(X\right)` when :math:`\mathrm{tol}` is zero, where :math:`\mathrm{size}\left(X\right)` is based on the size of the absolute values of the elements of the matrix :math:`X` containing the matrix :math:`R`. See Hammarling `et al.` (2013) for the motivation. If :math:`\mathrm{tol}` is :math:`-1.0` on entry then no deflation is attempted. The recommended value for :math:`\mathrm{tol}` is zero.
**a** : float, array-like, shape :math:`\left(n, n\right)`
The :math:`n\times n` matrix :math:`A`.
**b** : float, array-like, shape :math:`\left(n, n\right)`
The :math:`n\times n` matrix :math:`B`.
**c** : float, array-like, shape :math:`\left(n, n\right)`
The :math:`n\times n` matrix :math:`C`.
**Returns**
**a** : float, ndarray, shape :math:`\left(n, n\right)`
:math:`\mathrm{a}` is used as internal workspace, but if :math:`\mathrm{jobvl} = \texttt{'V'}` or :math:`\mathrm{jobvr} = \texttt{'V'}`, :math:`\mathrm{a}` is restored on exit.
**b** : float, ndarray, shape :math:`\left(n, n\right)`
:math:`\mathrm{b}` is used as internal workspace, but is restored on exit.
**c** : float, ndarray, shape :math:`\left(n, n\right)`
:math:`\mathrm{c}` is used as internal workspace, but if :math:`\mathrm{jobvl} = \texttt{'V'}` or :math:`\mathrm{jobvr} = \texttt{'V'}`, :math:`\mathrm{c}` is restored on exit.
**alphar** : float, ndarray, shape :math:`\left(2\times n\right)`
:math:`\mathrm{alphar}[\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,2n`, contains the real part of :math:`\alpha_j` for the :math:`j`\ th eigenvalue pair :math:`\left(\alpha_j, \beta_j\right)` of the quadratic eigenvalue problem.
**alphai** : float, ndarray, shape :math:`\left(2\times n\right)`
:math:`\mathrm{alphai}[\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,2n`, contains the imaginary part of :math:`\alpha_j` for the :math:`j`\ th eigenvalue pair :math:`\left(\alpha_j, \beta_j\right)` of the quadratic eigenvalue problem. If :math:`\mathrm{alphai}[j-1]` is zero then the :math:`j`\ th eigenvalue is real; if :math:`\mathrm{alphai}[j-1]` is positive then the :math:`j`\ th and :math:`\left(j+1\right)`\ th eigenvalues are a complex conjugate pair, with :math:`\mathrm{alphai}[j]` negative.
**beta** : float, ndarray, shape :math:`\left(2\times n\right)`
:math:`\mathrm{beta}[\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,2n`, contains the second part of the :math:`j`\ th eigenvalue pair :math:`\left(\alpha_j, \beta_j\right)` of the quadratic eigenvalue problem, with :math:`\beta_j\geq 0`. Infinite eigenvalues have :math:`\beta_j` set to zero.
**vl** : None or float, ndarray, shape :math:`\left(n, :\right)`
If :math:`\mathrm{jobvl} = \texttt{'V'}`, the left eigenvectors :math:`y_j` are stored one after another in the columns of :math:`\mathrm{vl}`, in the same order as the corresponding eigenvalues. If the :math:`j`\ th eigenvalue is real, then :math:`y_j = \mathrm{vl}[:-1,j-1]`, the :math:`j`\ th column of :math:`\mathrm{vl}`. If the :math:`j`\ th and :math:`\left(j+1\right)`\ th eigenvalues form a complex conjugate pair, then :math:`y_j = \mathrm{vl}[:-1,j-1]+i\times \mathrm{vl}[:-1,j]` and :math:`y_{{j+1}} = \mathrm{vl}[:-1,j-1]-i\times \mathrm{vl}[:-1,j]`. Each eigenvector will be normalized with length unity and with the element of largest modulus real and positive.
If :math:`\mathrm{jobvl} = \texttt{'N'}`, :math:`\mathrm{vl}` is not referenced.
**vr** : None or float, ndarray, shape :math:`\left(n, :\right)`
If :math:`\mathrm{jobvr} = \texttt{'V'}`, the right eigenvectors :math:`x_j` are stored one after another in the columns of :math:`\mathrm{vr}`, in the same order as the corresponding eigenvalues. If the :math:`j`\ th eigenvalue is real, then :math:`x_j = \mathrm{vr}[:-1,j-1]`, the :math:`j`\ th column of :math:`\mathrm{vr}`. If the :math:`j`\ th and :math:`\left(j+1\right)`\ th eigenvalues form a complex conjugate pair, then :math:`x_j = \mathrm{vr}[:-1,j-1]+i\times \mathrm{vr}[:-1,j]` and :math:`x_{{j+1}} = \mathrm{vr}[:-1,j-1]-i\times \mathrm{vr}[:-1,j]`. Each eigenvector will be normalized with length unity and with the element of largest modulus real and positive.
If :math:`\mathrm{jobvr} = \texttt{'N'}`, :math:`\mathrm{vr}` is not referenced.
**s** : None or float, ndarray, shape :math:`\left(:\right)`
If :math:`\mathrm{sense} = 1`, :math:`5`, :math:`6` or :math:`7`, :math:`\mathrm{s}[j-1]` contains the condition number estimate for the :math:`j`\ th eigenvalue (large condition numbers imply that the problem is near one with multiple eigenvalues). Infinite condition numbers are returned as the largest model `float` number (:meth:`machine.real_largest <naginterfaces.library.machine.real_largest>`).
If :math:`\mathrm{sense} = 0`, :math:`2`, :math:`3` or :math:`4`, :math:`\mathrm{s}` is not referenced.
**bevl** : None or float, ndarray, shape :math:`\left(:\right)`
If :math:`\mathrm{sense} = 2`, :math:`4`, :math:`5` or :math:`7`, :math:`\mathrm{bevl}[j-1]` contains the backward error estimate for the computed left eigenpair :math:`\left(\lambda_j, y_j\right)`.
If :math:`\mathrm{sense} = 0`, :math:`1`, :math:`3` or :math:`6`, :math:`\mathrm{bevl}` is not referenced.
**bevr** : None or float, ndarray, shape :math:`\left(:\right)`
If :math:`\mathrm{sense} = 3`, :math:`4`, :math:`6` or :math:`7`, :math:`\mathrm{bevr}[j-1]` contains the backward error estimate for the computed right eigenpair :math:`\left(\lambda_j, x_j\right)`.
If :math:`\mathrm{sense} = 0`, :math:`1`, :math:`2` or :math:`5`, :math:`\mathrm{bevr}` is not referenced.
**iwarn** : int
:math:`\mathrm{iwarn}` will be positive if there are warnings, otherwise :math:`\mathrm{iwarn}` will be :math:`0`.
If no exception or warning is raised then:
if :math:`\mathrm{iwarn} = 1` then one, or both, of the matrices :math:`A` and :math:`C` is zero. In this case no scaling is performed, even if :math:`\mathrm{scal} > 0`;
if :math:`\mathrm{iwarn} = 2` then the matrices :math:`A` and :math:`C` are singular, or nearly singular, so the problem is potentially ill-posed;
if :math:`\mathrm{iwarn} = 3` then both the conditions for :math:`\mathrm{iwarn} = 1` and :math:`\mathrm{iwarn} = 2` above, apply. If :math:`\mathrm{iwarn} = 4`, :math:`\left\lVert \mathrm{b}\right\rVert \geq 10\sqrt{\left\lVert A\right\rVert ·\left\lVert C\right\rVert }` and backward stability cannot be guaranteed.
If :math:`\mathrm{errno}` = 2, :math:`\mathrm{iwarn}` returns the value of :math:`\textit{info}` from :meth:`lapackeig.dgges3 <naginterfaces.library.lapackeig.dgges3>`.
If :math:`\mathrm{errno}` = 3, :math:`\mathrm{iwarn}` returns the value of :math:`\textit{info}` from :meth:`lapackeig.dggev3 <naginterfaces.library.lapackeig.dggev3>`.
.. _f02jc-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-19`)
On entry, :math:`\textit{ldvr} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{jobvr} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: when :math:`\mathrm{jobvr} = \texttt{'V'}`, :math:`\textit{ldvr}\geq n`.
(`errno` :math:`-17`)
On entry, :math:`\textit{ldvl} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{jobvl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: when :math:`\mathrm{jobvl} = \texttt{'V'}`, :math:`\textit{ldvl}\geq n`.
(`errno` :math:`-6`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`-4`)
On entry, :math:`\mathrm{sense} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{sense} = 0`, :math:`1`, :math:`2`, :math:`3`, :math:`4`, :math:`5`, :math:`6` or :math:`7`.
(`errno` :math:`-3`)
On entry, :math:`\mathrm{sense} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{jobvr} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: when :math:`\mathrm{jobvr} = \texttt{'N'}`, :math:`\mathrm{sense} = 0` or :math:`2`,
when :math:`\mathrm{jobvr} = \texttt{'V'}`, :math:`\mathrm{sense} = 1`, :math:`3`, :math:`4`, :math:`5`, :math:`6` or :math:`7`.
(`errno` :math:`-3`)
On entry, :math:`\mathrm{jobvr} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{jobvr} = \texttt{'N'}` or :math:`\texttt{'V'}`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{sense} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{jobvl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: when :math:`\mathrm{jobvl} = \texttt{'N'}`, :math:`\mathrm{sense} = 0` or :math:`3`,
when :math:`\mathrm{jobvl} = \texttt{'V'}`, :math:`\mathrm{sense} = 1`, :math:`2`, :math:`4`, :math:`5`, :math:`6` or :math:`7`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{jobvl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{jobvl} = \texttt{'N'}` or :math:`\texttt{'V'}`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{scal} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{scal} = 0`, :math:`1`, :math:`2`, :math:`3` or :math:`4`.
(`errno` :math:`1`)
The quadratic matrix polynomial is nonregular (singular).
(`errno` :math:`2`)
The :math:`QZ` iteration failed in :meth:`lapackeig.dgges3 <naginterfaces.library.lapackeig.dgges3>`.
:math:`\mathrm{iwarn}` returns the value of :math:`\textit{errno}` returned by :meth:`lapackeig.dgges3 <naginterfaces.library.lapackeig.dgges3>`. This failure is unlikely to happen, but if it does, please contact `NAG <https://www.nag.com>`__.
(`errno` :math:`3`)
The :math:`QZ` iteration failed in :meth:`lapackeig.dggev3 <naginterfaces.library.lapackeig.dggev3>`.
:math:`\mathrm{iwarn}` returns the value of :math:`\textit{errno}` returned by :meth:`lapackeig.dggev3 <naginterfaces.library.lapackeig.dggev3>`. This failure is unlikely to happen, but if it does, please contact `NAG <https://www.nag.com>`__.
.. _f02jc-py2-py-notes:
**Notes**
The quadratic eigenvalue problem is solved by linearizing the problem and solving the resulting :math:`2n\times 2n` generalized eigenvalue problem.
The linearization is chosen to have favourable conditioning and backward stability properties.
An initial preprocessing step is performed that reveals and deflates the zero and infinite eigenvalues contributed by singular leading and trailing matrices.
The algorithm is backward stable for problems that are not too heavily damped, that is :math:`\left\lVert B\right\rVert \leq 10\sqrt{\left\lVert A\right\rVert ·\left\lVert C\right\rVert }`.
Further details on the algorithm are given in Hammarling `et al.` (2013).
.. _f02jc-py2-py-references:
**References**
Fan, H -Y, Lin, W.-W and Van Dooren, P., 2004, `Normwise scaling of second order polynomial matrices`, SIAM J. Matrix Anal. Appl. (26, 1), 252--256
Gaubert, S and Sharify, M, 2009, `Tropical scaling of polynomial matrices`, Lecture Notes in Control and Information Sciences Series (389), 291--303, Springer--Verlag
Hammarling, S, Munro, C J and Tisseur, F, 2013, `An algorithm for the complete solution of quadratic eigenvalue problems`, ACM Trans. Math. Software. (39(3):18:1--18:119), http://eprints.maths.manchester.ac.uk/2061/
"""
raise NotImplementedError
[docs]def complex_gen_quad(scal, jobvl, jobvr, sense, tol, a, b, c):
r"""
``complex_gen_quad`` solves the quadratic eigenvalue problem
.. math::
\left(\lambda^2A+\lambda B+C\right)x = 0\text{,}
where :math:`A`, :math:`B` and :math:`C` are complex :math:`n\times n` matrices.
The function returns the :math:`2n` eigenvalues, :math:`\lambda_{\textit{j}}`, for :math:`\textit{j} = 1,2,\ldots,2n`, and can optionally return the corresponding right eigenvectors, :math:`x_j` and/or left eigenvectors, :math:`y_j` as well as estimates of the condition numbers of the computed eigenvalues and backward errors of the computed right and left eigenvectors.
A left eigenvector satisfies the equation
.. math::
y^\mathrm{H}\left(\lambda^2A+\lambda B+C\right) = 0\text{,}
where :math:`y^\mathrm{H}` is the complex conjugate transpose of :math:`y`.
:math:`\lambda` is represented as the pair :math:`\left(\alpha, \beta \right)`, such that :math:`\lambda = \alpha /\beta`.
Note that the computation of :math:`\alpha /\beta` may overflow and indeed :math:`\beta` may be zero.
.. _f02jq-py2-py-doc:
For full information please refer to the NAG Library document for f02jq
https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02jqf.html
.. _f02jq-py2-py-parameters:
**Parameters**
**scal** : int
Determines the form of scaling to be performed on :math:`A`, :math:`B` and :math:`C`.
:math:`\mathrm{scal} = 0`
No scaling.
:math:`\mathrm{scal} = 1` (the recommended value)
Fan, Lin and Van Dooren scaling if :math:`\frac{\left\lVert B\right\rVert }{\sqrt{\left\lVert A\right\rVert \times \left\lVert C\right\rVert }} < 10` and no scaling otherwise where :math:`\left\lVert Z\right\rVert` is the Frobenius norm of :math:`Z`.
:math:`\mathrm{scal} = 2`
Fan, Lin and Van Dooren scaling.
:math:`\mathrm{scal} = 3`
Tropical scaling with largest root.
:math:`\mathrm{scal} = 4`
Tropical scaling with smallest root.
**jobvl** : str, length 1
If :math:`\mathrm{jobvl} = \texttt{'N'}`, do not compute left eigenvectors.
If :math:`\mathrm{jobvl} = \texttt{'V'}`, compute the left eigenvectors.
If :math:`\mathrm{sense} = 1`, :math:`2`, :math:`4`, :math:`5`, :math:`6` or :math:`7`, :math:`\mathrm{jobvl}` must be set to 'V'.
**jobvr** : str, length 1
If :math:`\mathrm{jobvr} = \texttt{'N'}`, do not compute right eigenvectors.
If :math:`\mathrm{jobvr} = \texttt{'V'}`, compute the right eigenvectors.
If :math:`\mathrm{sense} = 1`, :math:`3`, :math:`4`, :math:`5`, :math:`6` or :math:`7`, :math:`\mathrm{jobvr}` must be set to 'V'.
**sense** : int
Determines whether, or not, condition numbers and backward errors are computed.
:math:`\mathrm{sense} = 0`
Do not compute condition numbers, or backward errors.
:math:`\mathrm{sense} = 1`
Just compute condition numbers for the eigenvalues.
:math:`\mathrm{sense} = 2`
Just compute backward errors for the left eigenpairs.
:math:`\mathrm{sense} = 3`
Just compute backward errors for the right eigenpairs.
:math:`\mathrm{sense} = 4`
Compute backward errors for the left and right eigenpairs.
:math:`\mathrm{sense} = 5`
Compute condition numbers for the eigenvalues and backward errors for the left eigenpairs.
:math:`\mathrm{sense} = 6`
Compute condition numbers for the eigenvalues and backward errors for the right eigenpairs.
:math:`\mathrm{sense} = 7`
Compute condition numbers for the eigenvalues and backward errors for the left and right eigenpairs.
**tol** : float
:math:`\mathrm{tol}` is used as the tolerance for making decisions on rank in the deflation procedure. If :math:`\mathrm{tol}` is zero on entry then :math:`n\times \text{machine precision}` is used in place of :math:`\mathrm{tol}`, where machine precision is as returned by function :meth:`machine.precision <naginterfaces.library.machine.precision>`. A diagonal element of a triangular matrix, :math:`R`, is regarded as zero if :math:`\left\lvert r_{{jj}}\right\rvert \leq \mathrm{tol}\times \mathrm{size}\left(X\right)`, or :math:`n\times \text{machine precision}\times \mathrm{size}\left(X\right)` when :math:`\mathrm{tol}` is zero, where :math:`\mathrm{size}\left(X\right)` is based on the size of the absolute values of the elements of the matrix :math:`X` containing the matrix :math:`R`. See Hammarling `et al.` (2013) for the motivation. If :math:`\mathrm{tol}` is :math:`-1.0` on entry then no deflation is attempted. The recommended value for :math:`\mathrm{tol}` is zero.
**a** : complex, array-like, shape :math:`\left(n, n\right)`
The :math:`n\times n` matrix :math:`A`.
**b** : complex, array-like, shape :math:`\left(n, n\right)`
The :math:`n\times n` matrix :math:`B`.
**c** : complex, array-like, shape :math:`\left(n, n\right)`
The :math:`n\times n` matrix :math:`C`.
**Returns**
**a** : complex, ndarray, shape :math:`\left(n, n\right)`
:math:`\mathrm{a}` is used as internal workspace, but if :math:`\mathrm{jobvl} = \texttt{'V'}` or :math:`\mathrm{jobvr} = \texttt{'V'}`, :math:`\mathrm{a}` is restored on exit.
**b** : complex, ndarray, shape :math:`\left(n, n\right)`
:math:`\mathrm{b}` is used as internal workspace, but is restored on exit.
**c** : complex, ndarray, shape :math:`\left(n, n\right)`
:math:`\mathrm{c}` is used as internal workspace, but if :math:`\mathrm{jobvl} = \texttt{'V'}` or :math:`\mathrm{jobvr} = \texttt{'V'}`, :math:`\mathrm{c}` is restored on exit.
**alpha** : complex, ndarray, shape :math:`\left(2\times n\right)`
:math:`\mathrm{alpha}[\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,2n`, contains the first part of the :math:`j`\ th eigenvalue pair :math:`\left(\alpha_j, \beta_j\right)` of the quadratic eigenvalue problem.
**beta** : complex, ndarray, shape :math:`\left(2\times n\right)`
:math:`\mathrm{beta}[\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,2n`, contains the second part of the :math:`j`\ th eigenvalue pair :math:`\left(\alpha_j, \beta_j\right)` of the quadratic eigenvalue problem. Although :math:`\mathrm{beta}` is declared complex, it is actually real and non-negative. Infinite eigenvalues have :math:`\beta_j` set to zero.
**vl** : None or complex, ndarray, shape :math:`\left(n, :\right)`
If :math:`\mathrm{jobvl} = \texttt{'V'}`, the left eigenvectors :math:`y_j` are stored one after another in the columns of :math:`\mathrm{vl}`, in the same order as the corresponding eigenvalues. Each eigenvector will be normalized with length unity and with the element of largest modulus real and positive.
If :math:`\mathrm{jobvl} = \texttt{'N'}`, :math:`\mathrm{vl}` is not referenced.
**vr** : None or complex, ndarray, shape :math:`\left(n, :\right)`
If :math:`\mathrm{jobvr} = \texttt{'V'}`, the right eigenvectors :math:`x_j` are stored one after another in the columns of :math:`\mathrm{vr}`, in the same order as the corresponding eigenvalues. Each eigenvector will be normalized with length unity and with the element of largest modulus real and positive.
If :math:`\mathrm{jobvr} = \texttt{'N'}`, :math:`\mathrm{vr}` is not referenced.
**s** : None or float, ndarray, shape :math:`\left(:\right)`
If :math:`\mathrm{sense} = 1`, :math:`5`, :math:`6` or :math:`7`, :math:`\mathrm{s}[j-1]` contains the condition number estimate for the :math:`j`\ th eigenvalue (large condition numbers imply that the problem is near one with multiple eigenvalues). Infinite condition numbers are returned as the largest model real number (:meth:`machine.real_largest <naginterfaces.library.machine.real_largest>`).
If :math:`\mathrm{sense} = 0`, :math:`2`, :math:`3` or :math:`4`, :math:`\mathrm{s}` is not referenced.
**bevl** : None or float, ndarray, shape :math:`\left(:\right)`
If :math:`\mathrm{sense} = 2`, :math:`4`, :math:`5` or :math:`7`, :math:`\mathrm{bevl}[j-1]` contains the backward error estimate for the computed left eigenpair :math:`\left(\lambda_j, y_j\right)`.
If :math:`\mathrm{sense} = 0`, :math:`1`, :math:`3` or :math:`6`, :math:`\mathrm{bevl}` is not referenced.
**bevr** : None or float, ndarray, shape :math:`\left(:\right)`
If :math:`\mathrm{sense} = 3`, :math:`4`, :math:`6` or :math:`7`, :math:`\mathrm{bevr}[j-1]` contains the backward error estimate for the computed right eigenpair :math:`\left(\lambda_j, x_j\right)`.
If :math:`\mathrm{sense} = 0`, :math:`1`, :math:`2` or :math:`5`, :math:`\mathrm{bevr}` is not referenced.
**iwarn** : int
:math:`\mathrm{iwarn}` will be positive if there are warnings, otherwise :math:`\mathrm{iwarn}` will be :math:`0`.
If no exception or warning is raised then:
if :math:`\mathrm{iwarn} = 1` then one, or both, of the matrices :math:`A` and :math:`C` is zero. In this case no scaling is performed, even if :math:`\mathrm{scal} > 0`;
if :math:`\mathrm{iwarn} = 2` then the matrices :math:`A` and :math:`C` are singular, or nearly singular, so the problem is potentially ill-posed;
if :math:`\mathrm{iwarn} = 3` then both the conditions for :math:`\mathrm{iwarn} = 1` and :math:`\mathrm{iwarn} = 2` above, apply. If :math:`\mathrm{iwarn} = 4`, :math:`\left\lVert \mathrm{b}\right\rVert \geq 10\sqrt{\left\lVert A\right\rVert ·\left\lVert C\right\rVert }` and backward stability cannot be guaranteed.
If :math:`\mathrm{errno}` = 2, :math:`\mathrm{iwarn}` returns the value of :math:`\textit{info}` from :meth:`lapackeig.zgges3 <naginterfaces.library.lapackeig.zgges3>`.
If :math:`\mathrm{errno}` = 3, :math:`\mathrm{iwarn}` returns the value of :math:`\textit{info}` from :meth:`lapackeig.zggev3 <naginterfaces.library.lapackeig.zggev3>`.
.. _f02jq-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-18`)
On entry, :math:`\textit{ldvr} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{jobvr} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: when :math:`\mathrm{jobvr} = \texttt{'V'}`, :math:`\textit{ldvr}\geq n`.
(`errno` :math:`-16`)
On entry, :math:`\textit{ldvl} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{jobvl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: when :math:`\mathrm{jobvl} = \texttt{'V'}`, :math:`\textit{ldvl}\geq n`.
(`errno` :math:`-6`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`-4`)
On entry, :math:`\mathrm{sense} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{sense} = 0`, :math:`1`, :math:`2`, :math:`3`, :math:`4`, :math:`5`, :math:`6` or :math:`7`.
(`errno` :math:`-3`)
On entry, :math:`\mathrm{sense} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{jobvr} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: when :math:`\mathrm{jobvr} = \texttt{'N'}`, :math:`\mathrm{sense} = 0` or :math:`2`,
when :math:`\mathrm{jobvr} = \texttt{'V'}`, :math:`\mathrm{sense} = 1`, :math:`3`, :math:`4`, :math:`5`, :math:`6` or :math:`7`.
(`errno` :math:`-3`)
On entry, :math:`\mathrm{jobvr} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{jobvr} = \texttt{'N'}` or :math:`\texttt{'V'}`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{sense} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{jobvl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: when :math:`\mathrm{jobvl} = \texttt{'N'}`, :math:`\mathrm{sense} = 0` or :math:`3`,
when :math:`\mathrm{jobvl} = \texttt{'V'}`, :math:`\mathrm{sense} = 1`, :math:`2`, :math:`4`, :math:`5`, :math:`6` or :math:`7`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{jobvl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{jobvl} = \texttt{'N'}` or :math:`\texttt{'V'}`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{scal} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{scal} = 0`, :math:`1`, :math:`2`, :math:`3` or :math:`4`.
(`errno` :math:`1`)
The quadratic matrix polynomial is nonregular (singular).
(`errno` :math:`2`)
The :math:`QZ` iteration failed in :meth:`lapackeig.zgges3 <naginterfaces.library.lapackeig.zgges3>`.
:math:`\mathrm{iwarn}` returns the value of :math:`\textit{errno}` returned by :meth:`lapackeig.zgges3 <naginterfaces.library.lapackeig.zgges3>`. This failure is unlikely to happen, but if it does, please contact `NAG <https://www.nag.com>`__.
(`errno` :math:`3`)
The :math:`QZ` iteration failed in :meth:`lapackeig.zggev3 <naginterfaces.library.lapackeig.zggev3>`.
:math:`\mathrm{iwarn}` returns the value of :math:`\textit{errno}` returned by :meth:`lapackeig.zggev3 <naginterfaces.library.lapackeig.zggev3>`. This failure is unlikely to happen, but if it does, please contact `NAG <https://www.nag.com>`__.
.. _f02jq-py2-py-notes:
**Notes**
The quadratic eigenvalue problem is solved by linearizing the problem and solving the resulting :math:`2n\times 2n` generalized eigenvalue problem.
The linearization is chosen to have favourable conditioning and backward stability properties.
An initial preprocessing step is performed that reveals and deflates the zero and infinite eigenvalues contributed by singular leading and trailing matrices.
The algorithm is backward stable for problems that are not too heavily damped, that is :math:`\left\lVert B\right\rVert \leq 10\sqrt{\left\lVert A\right\rVert ·\left\lVert C\right\rVert }`.
Further details on the algorithm are given in Hammarling `et al.` (2013).
.. _f02jq-py2-py-references:
**References**
Fan, H -Y, Lin, W.-W and Van Dooren, P., 2004, `Normwise scaling of second order polynomial matrices`, SIAM J. Matrix Anal. Appl. (26, 1), 252--256
Gaubert, S and Sharify, M, 2009, `Tropical scaling of polynomial matrices`, Lecture Notes in Control and Information Sciences Series (389), 291--303, Springer--Verlag
Hammarling, S, Munro, C J and Tisseur, F, 2013, `An algorithm for the complete solution of quadratic eigenvalue problems`, ACM Trans. Math. Software. (39(3):18:1--18:119), http://eprints.maths.manchester.ac.uk/2061/
"""
raise NotImplementedError
[docs]def real_gen_partialsvd(m, n, k, ncv, av, data=None, io_manager=None):
r"""
``real_gen_partialsvd`` returns leading terms in the singular value decomposition (SVD) of a real general matrix and computes the corresponding left and right singular vectors.
.. _f02wg-py2-py-doc:
For full information please refer to the NAG Library document for f02wg
https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02wgf.html
.. _f02wg-py2-py-parameters:
**Parameters**
**m** : int
:math:`m`, the number of rows of the matrix :math:`A`.
**n** : int
:math:`n`, the number of columns of the matrix :math:`A`.
**k** : int
:math:`k`, the number of singular values to be computed.
**ncv** : int
This is the number of Lanczos basis vectors to use during the computation of the largest eigenvalues of :math:`A^\mathrm{T}A` (:math:`m\geq n`) or :math:`AA^\mathrm{T}` (:math:`m < n`).
At present there is no `a priori` analysis to guide the selection of :math:`\mathrm{ncv}` relative to :math:`\mathrm{k}`.
However, it is recommended that :math:`\mathrm{ncv}\geq 2\times \mathrm{k}+1`.
If many problems of the same type are to be solved, you should experiment with varying :math:`\mathrm{ncv}` while keeping :math:`\mathrm{k}` fixed for a given test problem.
This will usually decrease the required number of matrix-vector operations but it also increases the internal storage required to maintain the orthogonal basis vectors.
The optimal 'cross-over' with respect to CPU time is problem dependent and must be determined empirically.
**av** : callable (iflag, ax) = av(iflag, m, n, x, data=None)
:math:`\mathrm{av}` must return the vector result of the matrix-vector product :math:`Ax` or :math:`A^\mathrm{T}x`, as indicated by the input value of :math:`\mathrm{iflag}`, for the given vector :math:`x`.
**Parameters**
**iflag** : int
If :math:`\mathrm{iflag} = 1`, :math:`\mathrm{ax}` must return the :math:`m`-vector result of the matrix-vector product :math:`Ax`.
If :math:`\mathrm{iflag} = 2`, :math:`\mathrm{ax}` must return the :math:`n`-vector result of the matrix-vector product :math:`A^\mathrm{T}x`.
**m** : int
The number of rows of the matrix :math:`A`.
**n** : int
The number of columns of the matrix :math:`A`.
**x** : float, ndarray, shape :math:`\left(:\right)`
The vector to be pre-multiplied by the matrix :math:`A` or :math:`A^\mathrm{T}`.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**iflag** : int
May be used as a flag to indicate a failure in the computation of :math:`Ax` or :math:`A^\mathrm{T}x`. If :math:`\mathrm{iflag}` is negative on exit from :math:`\mathrm{av}`, ``real_gen_partialsvd`` will exit immediately with :math:`\textit{errno}` set to :math:`\mathrm{iflag}`.
**ax** : float, array-like, shape :math:`\left(:\right)`
If :math:`\mathrm{iflag} = 1`, contains the :math:`m`-vector result of the matrix-vector product :math:`Ax`.
If :math:`\mathrm{iflag} = 2`, contains the :math:`n`-vector result of the matrix-vector product :math:`A^\mathrm{T}x`.
**data** : arbitrary, optional
User-communication data for callback functions.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**nconv** : int
The number of converged singular values found.
**sigma** : float, ndarray, shape :math:`\left(\mathrm{ncv}\right)`
The :math:`\mathrm{nconv}` converged singular values are stored in the first :math:`\mathrm{nconv}` elements of :math:`\mathrm{sigma}`.
**u** : float, ndarray, shape :math:`\left(\mathrm{m}, \mathrm{ncv}\right)`
The left singular vectors corresponding to the singular values stored in :math:`\mathrm{sigma}`.
The :math:`\textit{i}`\ th element of the :math:`\textit{j}`\ th left singular vector :math:`u_{\textit{j}}` is stored in :math:`\mathrm{u}[\textit{i}-1,\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,\mathrm{nconv}`, for :math:`\textit{i} = 1,2,\ldots,m`.
**v** : float, ndarray, shape :math:`\left(\mathrm{n}, \mathrm{ncv}\right)`
The right singular vectors corresponding to the singular values stored in :math:`\mathrm{sigma}`.
The :math:`\textit{i}`\ th element of the :math:`\textit{j}`\ th right singular vector :math:`v_{\textit{j}}` is stored in :math:`\mathrm{v}[\textit{i}-1,\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,\mathrm{nconv}`, for :math:`\textit{i} = 1,2,\ldots,n`.
**resid** : float, ndarray, shape :math:`\left(\mathrm{ncv}\right)`
The residual :math:`\left\lVert Av_j-\sigma_ju_j\right\rVert`, for :math:`m\geq n`, or :math:`\left\lVert A^\mathrm{T}u_j-\sigma_jv_j\right\rVert`, for :math:`m < n`, for each of the converged singular values :math:`\sigma_j` and corresponding left and right singular vectors :math:`u_j` and :math:`v_j`.
.. _f02wg-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`i < 0`)
On output from user-defined function :math:`\mathrm{av}`, :math:`\mathrm{iflag}` was set to a negative value, :math:`\mathrm{iflag} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m}\geq 0`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`3`)
On entry, :math:`\mathrm{k} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{k} > 0`.
(`errno` :math:`4`)
On entry, :math:`\mathrm{k} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{ncv} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{k} < \mathrm{ncv}\leq \mathrm{min}\left(\mathrm{m}, \mathrm{n}\right)`.
(`errno` :math:`8`)
The maximum number of iterations has been reached. The maximum number of iterations :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle`. The number of converged eigenvalues :math:`\text{} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`9`)
No shifts could be applied during a cycle of the implicitly restarted Lanczos iteration.
(`errno` :math:`10`)
Could not build a full Lanczos factorization.
(`errno` :math:`11`)
The number of eigenvalues found to sufficient accuracy is zero.
(`errno` :math:`20`)
An error occurred during an internal call. Consider increasing the size of :math:`\mathrm{ncv}` relative to :math:`\mathrm{k}`.
.. _f02wg-py2-py-notes:
**Notes**
``real_gen_partialsvd`` computes a few, :math:`k`, of the largest singular values and corresponding vectors of an :math:`m\times n` matrix :math:`A`.
The value of :math:`k` should be small relative to :math:`m` and :math:`n`, for example :math:`k\sim O\left(\mathrm{min}\left(m, n\right)\right)`.
The full singular value decomposition (SVD) of an :math:`m\times n` matrix :math:`A` is given by
.. math::
A = U\Sigma V^\mathrm{T}\text{,}
where :math:`U` and :math:`V` are orthogonal and :math:`\Sigma` is an :math:`m\times n` diagonal matrix with real diagonal elements, :math:`\sigma_i`, such that
.. math::
\sigma_1\geq \sigma_2\geq \cdots \geq \sigma_{{\mathrm{min}\left(m, n\right)}}\geq 0\text{.}
The :math:`\sigma_i` are the `singular values` of :math:`A` and the first :math:`\mathrm{min}\left(m, n\right)` columns of :math:`U` and :math:`V` are the `left` and `right singular vectors` of :math:`A`.
If :math:`U_k`, :math:`V_k` denote the leading :math:`k` columns of :math:`U` and :math:`V` respectively, and if :math:`Σ_k` denotes the leading principal submatrix of :math:`Σ`, then
.. math::
A_k\equiv U_kΣ_kV^\mathrm{T}_k
is the best rank-:math:`k` approximation to :math:`A` in both the :math:`2`-norm and the Frobenius norm.
The singular values and singular vectors satisfy
.. math::
Av_i = \sigma_iu_i\quad \text{ and }\quad A^\mathrm{T}u_i = \sigma_iv_i\quad \text{ so that }\quad A^\mathrm{T}A\nu_i = \sigma_i^2\nu_i\text{ and }AA^\mathrm{T}u_i = \sigma_i^2u_i\text{,}
where :math:`u_i` and :math:`v_i` are the :math:`i`\ th columns of :math:`U_k` and :math:`V_k` respectively.
Thus, for :math:`m\geq n`, the largest singular values and corresponding right singular vectors are computed by finding eigenvalues and eigenvectors for the symmetric matrix :math:`A^\mathrm{T}A`.
For :math:`m < n`, the largest singular values and corresponding left singular vectors are computed by finding eigenvalues and eigenvectors for the symmetric matrix :math:`AA^\mathrm{T}`.
These eigenvalues and eigenvectors are found using functions from submodule :mod:`~naginterfaces.library.sparseig`.
You should read `the F12 Introduction <https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f12/f12intro.html>`__ for full details of the method used here.
The real matrix :math:`A` is not explicitly supplied to ``real_gen_partialsvd``.
Instead, you are required to supply a function, :math:`\mathrm{av}`, that must calculate one of the requested matrix-vector products :math:`Ax` or :math:`A^\mathrm{T}x` for a given real vector :math:`x` (of length :math:`n` or :math:`m` respectively).
.. _f02wg-py2-py-references:
**References**
Wilkinson, J H, 1978, `Singular Value Decomposition -- Basic Aspects`, Numerical Software -- Needs and Availability, (ed D A H Jacobs), Academic Press
"""
raise NotImplementedError
[docs]def real_triang_svd(n, a, b, wantq, wantp):
r"""
``real_triang_svd`` returns all, or part, of the singular value decomposition of a real upper triangular matrix.
.. _f02wu-py2-py-doc:
For full information please refer to the NAG Library document for f02wu
https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02wuf.html
.. _f02wu-py2-py-parameters:
**Parameters**
**n** : int
:math:`n`, the order of the matrix :math:`R`.
If :math:`\mathrm{n} = 0`, an immediate return is effected.
**a** : float, array-like, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
The leading :math:`n\times n` upper triangular part of the array :math:`\mathrm{a}` must contain the upper triangular matrix :math:`R`.
**b** : float, array-like, shape :math:`\left(:, \textit{ncolb}\right)`
Note: the required extent for this argument in dimension 1 is determined as follows: if :math:`\textit{ncolb} > 0`: :math:`\mathrm{n}`; otherwise: :math:`1`.
With :math:`\textit{ncolb} > 0`, the leading :math:`n\times \textit{ncolb}` part of the array :math:`\mathrm{b}` must contain the matrix to be transformed.
**wantq** : bool
Must be :math:`\mathbf{True}` if the matrix :math:`Q` is required.
If :math:`\mathrm{wantq} = \mathbf{False}`, the array :math:`\mathrm{q}` is not referenced.
**wantp** : bool
Must be :math:`\mathbf{True}` if the matrix :math:`P^\mathrm{T}` is required, in which case :math:`P^\mathrm{T}` is overwritten on the array :math:`\mathrm{a}`, otherwise :math:`\mathrm{wantp}` must be :math:`\mathbf{False}`.
**Returns**
**a** : float, ndarray, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
If :math:`\mathrm{wantp} = \mathbf{True}`, the :math:`n\times n` part of :math:`\mathrm{a}` will contain the :math:`n\times n` orthogonal matrix :math:`P^T`, otherwise the :math:`n\times n` upper triangular part of :math:`\mathrm{a}` is used as internal workspace, but the strictly lower triangular part of :math:`\mathrm{a}` is not referenced.
**b** : float, ndarray, shape :math:`\left(:, \textit{ncolb}\right)`
The leading :math:`n\times \textit{ncolb}` part of the array :math:`\mathrm{b}` is overwritten by the matrix :math:`Q^\mathrm{T}B`.
**q** : float, ndarray, shape :math:`\left(:, :\right)`
With :math:`\mathrm{wantq} = \mathbf{True}`, the leading :math:`n\times n` part of the array :math:`\mathrm{q}` will contain the orthogonal matrix :math:`Q`. Otherwise the array :math:`\mathrm{q}` is not referenced.
**sv** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception or warning is raised the array :math:`\mathrm{sv}` will contain the diagonal elements of the matrix.
If :math:`\mathrm{errno}` = 1 the array :math:`\mathrm{sv}` will contain the diagonal elements of the bidiagonal matrix :math:`E` in the factorization :math:`R = QEP^T`; the superdiagonal elements of :math:`E` will be contained in the first :math:`\mathrm{n}-1` elements of :math:`\mathrm{work}`.
**work** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
:math:`\mathrm{work}[0:\mathrm{n}-1]` contains the super-diagonal elements of the bidiagonal matrix :math:`E` computed during the bidiagonalization stage; :math:`\mathrm{work}[\mathrm{n}-1]` contains the total number of iterations taken by the :math:`QR` algorithm.
The rest of the array is used as internal workspace.
.. _f02wu-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-1`)
On entry, :math:`\textit{ldq} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{wantq} = \mathbf{True}`, :math:`\textit{ldq}\geq \mathrm{n}`.
(`errno` :math:`-1`)
On entry, :math:`\textit{ldb} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{ncolb} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\textit{ncolb} > 0`, :math:`\textit{ldb}\geq \mathrm{n}`.
(`errno` :math:`-1`)
On entry, :math:`\textit{ncolb} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ncolb}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`1`)
The :math:`QR` algorithm has failed to converge. :math:`\langle\mathit{\boldsymbol{value}}\rangle` singular values have **not** been found.
.. _f02wu-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
The :math:`n\times n` upper triangular matrix :math:`R` is factorized as
.. math::
R = QSP^\mathrm{T}\text{,}
where :math:`Q` and :math:`P` are :math:`n\times n` orthogonal matrices and :math:`S` is an :math:`n\times n` diagonal matrix with non-negative diagonal elements, :math:`\sigma_1,\sigma_2,\ldots,\sigma_n`, ordered such that
.. math::
\sigma_1\geq \sigma_2\geq \cdots \geq \sigma_n\geq 0\text{.}
The columns of :math:`Q` are the left-hand singular vectors of :math:`R`, the diagonal elements of :math:`S` are the singular values of :math:`R` and the columns of :math:`P` are the right-hand singular vectors of :math:`R`.
Either or both of :math:`Q` and :math:`P^\mathrm{T}` may be requested and the matrix :math:`C` given by
.. math::
C = Q^\mathrm{T}B\text{,}
where :math:`B` is an :math:`n\times \textit{ncolb}` given matrix, may also be requested.
The function obtains the singular value decomposition by first reducing :math:`R` to bidiagonal form by means of Givens plane rotations and then using the :math:`QR` algorithm to obtain the singular value decomposition of the bidiagonal form.
Good background descriptions to the singular value decomposition are given in Chan (1982), Dongarra `et al.` (1979), Golub and Van Loan (1996), Hammarling (1985) and Wilkinson (1978).
Note that if :math:`K` is any orthogonal diagonal matrix so that
.. math::
KK^\mathrm{T} = I
(that is the diagonal elements of :math:`K` are :math:`+1` or :math:`-1`) then
.. math::
A = \left({QK}\right)S\left({PK}\right)^\mathrm{T}
is also a singular value decomposition of :math:`A`.
.. _f02wu-py2-py-references:
**References**
Chan, T F, 1982, `An improved algorithm for computing the singular value decomposition`, ACM Trans. Math. Software (8), 72--83
Dongarra, J J, Moler, C B, Bunch, J R and Stewart, G W, 1979, `LINPACK Users' Guide`, SIAM, Philadelphia
Golub, G H and Van Loan, C F, 1996, `Matrix Computations`, (3rd Edition), Johns Hopkins University Press, Baltimore
Hammarling, S, 1985, `The singular value decomposition in multivariate statistics`, SIGNUM Newsl. (20(3)), 2--25
Wilkinson, J H, 1978, `Singular Value Decomposition -- Basic Aspects`, Numerical Software -- Needs and Availability, (ed D A H Jacobs), Academic Press
"""
raise NotImplementedError
[docs]def complex_triang_svd(n, a, b, wantq, wantp):
r"""
``complex_triang_svd`` returns all, or part, of the singular value decomposition of a complex upper triangular matrix.
.. _f02xu-py2-py-doc:
For full information please refer to the NAG Library document for f02xu
https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/f02/f02xuf.html
.. _f02xu-py2-py-parameters:
**Parameters**
**n** : int
:math:`n`, the order of the matrix :math:`R`.
If :math:`\mathrm{n} = 0`, an immediate return is effected.
**a** : complex, array-like, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
The leading :math:`n\times n` upper triangular part of the array :math:`\mathrm{a}` must contain the upper triangular matrix :math:`R`.
**b** : complex, array-like, shape :math:`\left(:, \textit{ncolb}\right)`
Note: the required extent for this argument in dimension 1 is determined as follows: if :math:`\textit{ncolb} > 0`: :math:`\mathrm{n}`; otherwise: :math:`1`.
If :math:`\textit{ncolb} > 0`, the leading :math:`n\times \textit{ncolb}` part of the array :math:`\mathrm{b}` must contain the matrix to be transformed.
**wantq** : bool
Must be :math:`\mathbf{True}` if the matrix :math:`Q` is required.
If :math:`\mathrm{wantq} = \mathbf{False}` then the array :math:`\mathrm{q}` is not referenced.
**wantp** : bool
Must be :math:`\mathbf{True}` if the matrix :math:`P^\mathrm{H}` is required, in which case :math:`P^\mathrm{H}` is returned in the array :math:`\mathrm{a}`, otherwise :math:`\mathrm{wantp}` must be :math:`\mathbf{False}`.
**Returns**
**a** : complex, ndarray, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
If :math:`\mathrm{wantp} = \mathbf{True}`, the :math:`n\times n` part of :math:`\mathrm{a}` will contain the :math:`n\times n` unitary matrix :math:`P^H`, otherwise the :math:`n\times n` upper triangular part of :math:`\mathrm{a}` is used as internal workspace, but the strictly lower triangular part of :math:`\mathrm{a}` is not referenced.
**b** : complex, ndarray, shape :math:`\left(:, \textit{ncolb}\right)`
Is overwritten by the :math:`n\times \textit{ncolb}` matrix :math:`Q^\mathrm{H}B`.
**q** : complex, ndarray, shape :math:`\left(:, :\right)`
If :math:`\mathrm{wantq} = \mathbf{True}`, the leading :math:`n\times n` part of the array :math:`\mathrm{q}` will contain the unitary matrix :math:`Q`. Otherwise the array :math:`\mathrm{q}` is not referenced.
**sv** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
The :math:`n` diagonal elements of the matrix :math:`S`.
**rwork** : float, ndarray, shape :math:`\left(:\right)`
:math:`\mathrm{rwork}[\mathrm{n}-1]` contains the total number of iterations taken by the :math:`QR` algorithm.
The rest of the array is used as workspace.
.. _f02xu-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-1`)
On entry, :math:`\textit{ldq} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{wantq} = \mathbf{True}`, :math:`\textit{ldq}\geq \mathrm{n}`.
(`errno` :math:`-1`)
On entry, :math:`\textit{ldb} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{ncolb} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\textit{ncolb} > 0`, :math:`\textit{ldb}\geq \mathrm{n}`.
(`errno` :math:`-1`)
On entry, :math:`\textit{ncolb} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ncolb}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`1`)
The :math:`QR` algorithm has failed to converge. :math:`\langle\mathit{\boldsymbol{value}}\rangle` singular values have **not** been found.
.. _f02xu-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
The :math:`n\times n` upper triangular matrix :math:`R` is factorized as
.. math::
R = QSP^\mathrm{H}\text{,}
where :math:`Q` and :math:`P` are :math:`n\times n` unitary matrices and :math:`S` is an :math:`n\times n` diagonal matrix with real non-negative diagonal elements, :math:`sv_1,sv_2,\ldots,sv_n`, ordered such that
.. math::
sv_1\geq sv_2\geq \cdots \geq sv_n\geq 0\text{.}
The columns of :math:`Q` are the left-hand singular vectors of :math:`R`, the diagonal elements of :math:`S` are the singular values of :math:`R` and the columns of :math:`P` are the right-hand singular vectors of :math:`R`.
Either or both of :math:`Q` and :math:`P^\mathrm{H}` may be requested and the matrix :math:`C` given by
.. math::
C = Q^\mathrm{H}B\text{,}
where :math:`B` is an :math:`n\times \textit{ncolb}` given matrix, may also be requested.
``complex_triang_svd`` obtains the singular value decomposition by first reducing :math:`R` to bidiagonal form by means of Givens plane rotations and then using the :math:`QR` algorithm to obtain the singular value decomposition of the bidiagonal form.
Good background descriptions to the singular value decomposition are given in Dongarra `et al.` (1979), Hammarling (1985) and Wilkinson (1978).
Note that if :math:`K` is any unitary diagonal matrix so that
.. math::
KK^\mathrm{H} = I\text{,}
then
.. math::
A = \left({QK}\right)S\left({PK}\right)^H
is also a singular value decomposition of :math:`A`.
.. _f02xu-py2-py-references:
**References**
Dongarra, J J, Moler, C B, Bunch, J R and Stewart, G W, 1979, `LINPACK Users' Guide`, SIAM, Philadelphia
Hammarling, S, 1985, `The singular value decomposition in multivariate statistics`, SIGNUM Newsl. (20(3)), 2--25
Wilkinson, J H, 1978, `Singular Value Decomposition -- Basic Aspects`, Numerical Software -- Needs and Availability, (ed D A H Jacobs), Academic Press
"""
raise NotImplementedError
