# Source code for naginterfaces.library.eigen

# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 29.0 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/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/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/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

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

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/flhtml/f02/f02intro.html
"""

[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/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/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/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/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/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/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/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**
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/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/flhtml/f02/f02ekf.html#fcomments>__ while a fuller discussion is provided in the F12 Introduction <https://www.nag.com/numeric/nl/nagdoc_29/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/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/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/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/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/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::

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/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/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/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**
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/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/flhtml/f02/f02fkf.html#fcomments>__ while a fuller discussion is provided in the F12 Introduction <https://www.nag.com/numeric/nl/nagdoc_29/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/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/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/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/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/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::

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/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/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/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