# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 30.0 `linsys` Chapter.
``linsys`` - Simultaneous Linear Equations
This module is concerned with the solution of the matrix equation :math:`AX = B`, where :math:`B` may be a single vector or a matrix of multiple right-hand sides.
The matrix :math:`A` may be real, complex, symmetric, Hermitian, positive definite, positive definite Toeplitz or banded.
It may also be rectangular, in which case a least squares solution is obtained.
Much of the functionality of this module has been superseded by functions from submodule :mod:`~naginterfaces.library.lapacklin` and submodule :mod:`~naginterfaces.library.lapackeig` (LAPACK routines) as those modules have grown and have included driver and expert driver functions.
For a general introduction to sparse systems of equations, see `the F11 Introduction <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f11/f11intro.html>`__, which provides functions for large sparse systems. Some functions for sparse problems are also included in this module; they are described in `Sparse Matrix Functions <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04intro.html#recomm_4>`__.
Functionality Index
-------------------
**Black Box functions,** :math:`Ax = b`
complex general band matrix: :meth:`complex_band_solve`
complex general matrix: :meth:`complex_square_solve`
complex general tridiagonal matrix: :meth:`complex_tridiag_solve`
complex Hermitian matrix
packed matrix format: :meth:`complex_herm_packed_solve`
standard matrix format: :meth:`complex_herm_solve`
complex Hermitian positive definite band matrix: :meth:`complex_posdef_band_solve`
complex Hermitian positive definite matrix
packed matrix format: :meth:`complex_posdef_packed_solve`
standard matrix format: :meth:`complex_posdef_solve`
complex Hermitian positive definite tridiagonal matrix: :meth:`complex_posdef_tridiag_solve`
complex symmetric matrix
packed matrix format: :meth:`complex_symm_packed_solve`
standard matrix format: :meth:`complex_symm_solve`
real general band matrix: :meth:`real_band_solve`
real general matrix
multiple right-hand sides, standard precision: :meth:`real_square_solve`
real general tridiagonal matrix: :meth:`real_tridiag_solve`
real symmetric matrix
packed matrix format: :meth:`real_symm_packed_solve`
standard matrix format: :meth:`real_symm_solve`
real symmetric positive definite band matrix: :meth:`real_posdef_band_solve`
real symmetric positive definite matrix
multiple right-hand sides, standard precision: :meth:`real_posdef_solve`
packed matrix format: :meth:`real_posdef_packed_solve`
real symmetric positive definite Toeplitz matrix
general right-hand side: :meth:`real_toeplitz_solve`
Yule--Walker equations: :meth:`real_toeplitz_yule`
real symmetric positive definite tridiagonal matrix: :meth:`real_posdef_tridiag_solve`
**General Purpose functions,** :math:`Ax = b`
real almost block-diagonal matrix: :meth:`real_blkdiag_fac_solve`
real band symmetric positive definite matrix, variable bandwidth: :meth:`real_posdef_vband_solve`
real sparse matrix
direct method: :meth:`real_sparse_fac_solve`
iterative method: :meth:`real_gen_sparse_lsqsol`
real symmetric positive definite Toeplitz matrix
general right-hand side, update solution: :meth:`real_toeplitz_update`
Yule--Walker equations, update solution: :meth:`real_toeplitz_yule_update`
real tridiagonal matrix: :meth:`real_tridiag_fac_solve`
**Least squares and Homogeneous Equations**
real :math:`m\times n` matrix
:math:`m\geq n`, rank :math:`\text{} = n` or minimal solution: :meth:`real_gen_solve`
rank :math:`\text{} = n`, iterative refinement: :meth:`real_gen_lsqsol`
real sparse matrix: :meth:`real_gen_sparse_lsqsol`
**Service Functions**
complex rectangular matrix
norm and condition number estimation: :meth:`complex_gen_norm_rcomm`
real matrix
covariance matrix for linear least squares problems: :meth:`real_gen_lsq_covmat`
real rectangular matrix
norm and condition number estimation: :meth:`real_gen_norm_rcomm`
For full information please refer to the NAG Library document
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04intro.html
"""
# NAG Copyright 2017-2024.
[docs]def real_gen_lsqsol(a, b, eps):
r"""
``real_gen_lsqsol`` calculates the accurate least squares solution of a set of :math:`m` linear equations in :math:`n` unknowns, :math:`m\geq n` and rank :math:`\text{} = n`, with multiple right-hand sides, :math:`AX = B`, using a :math:`QR` factorization and iterative refinement.
.. _f04am-py2-py-doc:
For full information please refer to the NAG Library document for f04am
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04amf.html
.. _f04am-py2-py-parameters:
**Parameters**
**a** : float, array-like, shape :math:`\left(m, n\right)`
The :math:`m\times n` matrix :math:`A`.
**b** : float, array-like, shape :math:`\left(m, \textit{ir}\right)`
The :math:`m\times r` right-hand side matrix :math:`B`.
**eps** : float
Must be set to the value of the machine precision.
**Returns**
**x** : float, ndarray, shape :math:`\left(n, \textit{ir}\right)`
The :math:`n\times r` solution matrix :math:`X`.
**qr** : float, ndarray, shape :math:`\left(m, n\right)`
Details of the :math:`QR` factorization.
**alpha** : float, ndarray, shape :math:`\left(n\right)`
The diagonal elements of the upper triangular matrix :math:`R`.
**ipiv** : int, ndarray, shape :math:`\left(n\right)`
Details of the column interchanges.
.. _f04am-py2-py-errors:
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`1`)
The rank of :math:`\mathrm{a}` is less than :math:`\textit{n}`. The problem does not have a unique solution.
(`errno` :math:`2`)
The iterative refinement fails to converge. The matrix :math:`A` is too ill-conditioned.
.. _f04am-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
To compute the least squares solution to a set of :math:`m` linear equations in :math:`n` unknowns :math:`\left(m\geq n\right)` :math:`AX = B`, ``real_gen_lsqsol`` first computes a :math:`QR` factorization of :math:`A` with column pivoting, :math:`AP = QR`, where :math:`R` is upper triangular, :math:`Q` is an :math:`m\times m` orthogonal matrix, and :math:`P` is a permutation matrix. :math:`Q^\mathrm{T}` is applied to the :math:`m\times r` right-hand side matrix :math:`B` to give :math:`C = Q^\mathrm{T}B`, and the :math:`n\times r` solution matrix :math:`X` is calculated, to a first approximation, by back-substitution in :math:`RX = C`.
The residual matrix :math:`S = B-AX` is calculated using additional precision, and a correction :math:`D` to :math:`X` is computed as the least squares solution to :math:`AD = S`. :math:`X` is replaced by :math:`X+D` and this iterative refinement of the solution is repeated until full machine accuracy has been obtained.
.. _f04am-py2-py-references:
**References**
Wilkinson, J H and Reinsch, C, 1971, `Handbook for Automatic Computation II, Linear Algebra`, Springer--Verlag
"""
raise NotImplementedError
[docs]def real_sparse_fac_solve(a, comm, rhs, mtype):
r"""
``real_sparse_fac_solve`` calculates the approximate solution of a set of real sparse linear equations with a single right-hand side, :math:`Ax = b` or :math:`A^\mathrm{T}x = b`, where :math:`A` has been factorized by :meth:`matop.real_gen_sparse_lu <naginterfaces.library.matop.real_gen_sparse_lu>` or :meth:`matop.real_gen_sparse_lu_reuse <naginterfaces.library.matop.real_gen_sparse_lu_reuse>`.
.. _f04ax-py2-py-doc:
For full information please refer to the NAG Library document for f04ax
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04axf.html
.. _f04ax-py2-py-parameters:
**Parameters**
**a** : float, array-like, shape :math:`\left(\textit{licn}\right)`
The nonzero elements in the factorization of the matrix :math:`A`, as returned by :meth:`matop.real_gen_sparse_lu <naginterfaces.library.matop.real_gen_sparse_lu>` or :meth:`matop.real_gen_sparse_lu_reuse <naginterfaces.library.matop.real_gen_sparse_lu_reuse>`.
**comm** : dict, communication object
Communication structure.
This argument must have been initialized by a prior call to :meth:`matop.real_gen_sparse_lu <naginterfaces.library.matop.real_gen_sparse_lu>`.
**rhs** : float, array-like, shape :math:`\left(n\right)`
The right-hand side vector :math:`b`.
**mtype** : int
:math:`\mathrm{mtype}` specifies the task to be performed.
:math:`\mathrm{mtype} = 1`
Solve :math:`Ax = b`.
:math:`\mathrm{mtype}\neq 1`
Solve :math:`A^\mathrm{T}x = b`.
**Returns**
**rhs** : float, ndarray, shape :math:`\left(n\right)`
:math:`\mathrm{rhs}` is overwritten by the solution vector :math:`x`.
**resid** : float
The value of the maximum residual, :math:`\mathrm{max}\left(\left\lvert b_i-\sum_j^{{}}a_{{ij}}x_j\right\rvert \right)`, over all the unsatisfied equations, in case :meth:`matop.real_gen_sparse_lu <naginterfaces.library.matop.real_gen_sparse_lu>` or :meth:`matop.real_gen_sparse_lu_reuse <naginterfaces.library.matop.real_gen_sparse_lu_reuse>` has been used to factorize a singular or rectangular matrix.
.. _f04ax-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
.. _f04ax-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
To solve a system of real linear equations :math:`Ax = b` or :math:`A^\mathrm{T}x = b`, where :math:`A` is a general sparse matrix, :math:`A` must first be factorized by :meth:`matop.real_gen_sparse_lu <naginterfaces.library.matop.real_gen_sparse_lu>` or :meth:`matop.real_gen_sparse_lu_reuse <naginterfaces.library.matop.real_gen_sparse_lu_reuse>`. ``real_sparse_fac_solve`` then computes :math:`x` by block forward or backward substitution using simple forward and backward substitution within each diagonal block.
The method is fully described in Duff (1977).
A more recent method is available through solver function :meth:`sparse.direct_real_gen_solve <naginterfaces.library.sparse.direct_real_gen_solve>` and factorization function :meth:`sparse.direct_real_gen_lu <naginterfaces.library.sparse.direct_real_gen_lu>`.
.. _f04ax-py2-py-references:
**References**
Duff, I S, 1977, `MA28 -- a set of Fortran subroutines for sparse unsymmetric linear equations`, AERE Report R8730, HMSO
"""
raise NotImplementedError
[docs]def real_square_solve(n, a, b):
r"""
``real_square_solve`` computes the solution to a real system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` matrix and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04ba-py2-py-doc:
For full information please refer to the NAG Library document for f04ba
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04baf.html
.. _f04ba-py2-py-parameters:
**Parameters**
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**a** : float, array-like, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
The :math:`n\times n` coefficient matrix :math:`A`.
**b** : float, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**a** : float, ndarray, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
If no exception is raised, the factors :math:`L` and :math:`U` from the factorization :math:`A = PLU`. The unit diagonal elements of :math:`L` are not stored.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, the pivot indices that define the permutation matrix :math:`P`; at the :math:`i`\ th step row :math:`i` of the matrix was interchanged with row :math:`\mathrm{ipiv}[i-1]`. :math:`\mathrm{ipiv}[i-1] = i` indicates a row interchange was not required.
**b** : float, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no constraints are violated, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04ba-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-2`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\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:`\left(i > 0\right)\text{ and }\left(i \leq \mathrm{n}\right)`)
Diagonal element :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the upper triangular factor is zero. The factorization has been completed, but the solution could not be computed.
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04ba-py2-py-notes:
**Notes**
The :math:`LU` decomposition with partial pivoting and row interchanges is used to factor :math:`A` as :math:`A = PLU`, where :math:`P` is a permutation matrix, :math:`L` is unit lower triangular, and :math:`U` is upper triangular.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04ba-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def real_band_solve(n, kl, ku, ab, b):
r"""
``real_band_solve`` computes the solution to a real system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` band matrix, with :math:`k_l` subdiagonals and :math:`k_u` superdiagonals, and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04bb-py2-py-doc:
For full information please refer to the NAG Library document for f04bb
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04bbf.html
.. _f04bb-py2-py-parameters:
**Parameters**
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**kl** : int
The number of subdiagonals :math:`k_l`, within the band of :math:`A`.
**ku** : int
The number of superdiagonals :math:`k_u`, within the band of :math:`A`.
**ab** : float, array-like, shape :math:`\left(2\times \mathrm{kl}+\mathrm{ku}+1, \mathrm{n}\right)`
The :math:`n\times n` matrix :math:`A`.
See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04bbf.html#fcomments>`__ for further details.
**b** : float, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**ab** : float, ndarray, shape :math:`\left(2\times \mathrm{kl}+\mathrm{ku}+1, \mathrm{n}\right)`
If no exception is raised, :math:`\mathrm{ab}` is overwritten by details of the factorization.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, the pivot indices that define the permutation matrix :math:`P`; at the :math:`i`\ th step row :math:`i` of the matrix was interchanged with row :math:`\mathrm{ipiv}[i-1]`. :math:`\mathrm{ipiv}[i-1] = i` indicates a row interchange was not required.
**b** : float, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no constraints are violated, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04bb-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-6`)
On entry, :math:`\textit{ldab} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{kl} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ku} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ldab}\geq 2\times \mathrm{kl}+\mathrm{ku}+1`.
(`errno` :math:`-4`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-3`)
On entry, :math:`\mathrm{ku} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ku}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{kl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{kl}\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:`\left(i > 0\right)\text{ and }\left(i \leq \mathrm{n}\right)`)
Diagonal element :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the upper triangular factor is zero. The factorization has been completed, but the solution could not be computed.
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04bb-py2-py-notes:
**Notes**
The :math:`LU` decomposition with partial pivoting and row interchanges is used to factor :math:`A` as :math:`A = PLU`, where :math:`P` is a permutation matrix, :math:`L` is the product of permutation matrices and unit lower triangular matrices with :math:`k_l` subdiagonals, and :math:`U` is upper triangular with :math:`\left(k_l+k_u\right)` superdiagonals.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04bb-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def real_tridiag_solve(n, dl, d, du, b):
r"""
``real_tridiag_solve`` computes the solution to a real system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` tridiagonal matrix and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04bc-py2-py-doc:
For full information please refer to the NAG Library document for f04bc
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04bcf.html
.. _f04bc-py2-py-parameters:
**Parameters**
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**dl** : float, array-like, shape :math:`\left(\mathrm{n}-1\right)`
Must contain the :math:`\left(n-1\right)` subdiagonal elements of the matrix :math:`A`.
**d** : float, array-like, shape :math:`\left(\mathrm{n}\right)`
Must contain the :math:`n` diagonal elements of the matrix :math:`A`.
**du** : float, array-like, shape :math:`\left(\mathrm{n}-1\right)`
Must contain the :math:`\left(n-1\right)` superdiagonal elements of the matrix :math:`A`
**b** : float, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**dl** : float, ndarray, shape :math:`\left(\mathrm{n}-1\right)`
If no exception is raised, :math:`\mathrm{dl}` is overwritten by the :math:`\left(n-1\right)` multipliers that define the matrix :math:`L` from the :math:`LU` factorization of :math:`A`.
**d** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, :math:`\mathrm{d}` is overwritten by the :math:`n` diagonal elements of the upper triangular matrix :math:`U` from the :math:`LU` factorization of :math:`A`.
**du** : float, ndarray, shape :math:`\left(\mathrm{n}-1\right)`
If no exception is raised, :math:`\mathrm{du}` is overwritten by the :math:`\left(n-1\right)` elements of the first superdiagonal of :math:`U`.
**du2** : float, ndarray, shape :math:`\left(\mathrm{n}-2\right)`
If no exception is raised, :math:`\mathrm{du2}` returns the :math:`\left(n-2\right)` elements of the second superdiagonal of :math:`U`.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, the pivot indices that define the permutation matrix :math:`P`; at the :math:`i`\ th step row :math:`i` of the matrix was interchanged with row :math:`\mathrm{ipiv}[i-1]`. :math:`\mathrm{ipiv}[i-1]` will always be either :math:`i` or :math:`\left(i+1\right)`; :math:`\mathrm{ipiv}[i-1] = i` indicates a row interchange was not required.
**b** : float, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no constraints are violated, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04bc-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-2`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\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:`\left(i > 0\right)\text{ and }(\mathrm{n})`)
Diagonal element :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the upper triangular factor is zero. The factorization has been completed, but the solution could not be computed.
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04bc-py2-py-notes:
**Notes**
The :math:`LU` decomposition with partial pivoting and row interchanges is used to factor :math:`A` as :math:`A = PLU`, where :math:`P` is a permutation matrix, :math:`L` is unit lower triangular with at most one nonzero subdiagonal element, and :math:`U` is an upper triangular band matrix with two superdiagonals.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
Note that the equations :math:`A^\mathrm{T}X = B` may be solved by interchanging the order of the arguments :math:`\mathrm{du}` and :math:`\mathrm{dl}`.
.. _f04bc-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def real_posdef_solve(uplo, a, b):
r"""
``real_posdef_solve`` computes the solution to a real system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` symmetric positive definite matrix and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04bd-py2-py-doc:
For full information please refer to the NAG Library document for f04bd
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04bdf.html
.. _f04bd-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**a** : float, array-like, shape :math:`\left(n, n\right)`
The :math:`n\times n` symmetric matrix :math:`A`.
If :math:`\mathrm{uplo} = \texttt{'U'}`, the leading :math:`\textit{n}` by :math:`\textit{n}` upper triangular part of :math:`\mathrm{a}` contains the upper triangular part of the matrix :math:`A`, and the strictly lower triangular part of :math:`\mathrm{a}` is not referenced.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the leading :math:`\textit{n}` by :math:`\textit{n}` lower triangular part of :math:`\mathrm{a}` contains the lower triangular part of the matrix :math:`A`, and the strictly upper triangular part of :math:`\mathrm{a}` is not referenced.
**b** : float, array-like, shape :math:`\left(n, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**a** : float, ndarray, shape :math:`\left(n, n\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the factor :math:`U` or :math:`L` from the Cholesky factorization :math:`A = U^\mathrm{T}U` or :math:`A = LL^\mathrm{T}`.
**b** : float, ndarray, shape :math:`\left(n, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04bd-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-3`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo} \neq \texttt{'U'}` or :math:`\texttt{'L'}`: :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq n\right)`)
The principal minor of order :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the matrix :math:`A` is not positive definite. The factorization has not been completed and the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`n+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04bd-py2-py-notes:
**Notes**
The Cholesky factorization is used to factor :math:`A` as :math:`A = U^\mathrm{T}U`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LL^\mathrm{T}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` is an upper triangular matrix and :math:`L` is a lower triangular matrix.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04bd-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def real_posdef_packed_solve(uplo, ap, b):
r"""
``real_posdef_packed_solve`` computes the solution to a real system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` symmetric positive definite matrix, stored in packed format, and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04be-py2-py-doc:
For full information please refer to the NAG Library document for f04be
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04bef.html
.. _f04be-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**ap** : float, array-like, shape :math:`\left(n\times \left(n+1\right)/2\right)`
The :math:`n\times n` symmetric matrix :math:`A`. The upper or lower triangular part of the symmetric matrix is packed column-wise in a linear array. The :math:`j`\ th column of :math:`A` is stored in the array :math:`\mathrm{ap}` as follows:
**b** : float, array-like, shape :math:`\left(n, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**ap** : float, ndarray, shape :math:`\left(n\times \left(n+1\right)/2\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the factor :math:`U` or :math:`L` from the Cholesky factorization :math:`A = U^\mathrm{T}U` or :math:`A = LL^\mathrm{T}`, in the same storage format as :math:`A`.
**b** : float, ndarray, shape :math:`\left(n, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04be-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-3`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo}` not one of 'U' or 'u' or 'L' or 'l': :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq n\right)`)
The principal minor of order :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the matrix :math:`A` is not positive definite. The factorization has not been completed and the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`n+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04be-py2-py-notes:
**Notes**
The Cholesky factorization is used to factor :math:`A` as :math:`A = U^\mathrm{T}U`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LL^\mathrm{T}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` is an upper triangular matrix and :math:`L` is a lower triangular matrix.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04be-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def real_posdef_band_solve(uplo, kd, ab, b):
r"""
``real_posdef_band_solve`` computes the solution to a real system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` symmetric positive definite band matrix of band width :math:`2k+1`, and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04bf-py2-py-doc:
For full information please refer to the NAG Library document for f04bf
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04bff.html
.. _f04bf-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**kd** : int
The number of superdiagonals :math:`k` (and the number of subdiagonals) of the band matrix :math:`A`.
**ab** : float, array-like, shape :math:`\left(\mathrm{kd}+1, n\right)`
The :math:`n\times n` symmetric band matrix :math:`A`. The upper or lower triangular part of the symmetric matrix is stored in the first :math:`\mathrm{kd}+1` rows of the array. The :math:`j`\ th column of :math:`A` is stored in the :math:`j`\ th column of the array :math:`\mathrm{ab}` as follows:
See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04bff.html#fcomments>`__ below for further details.
**b** : float, array-like, shape :math:`\left(n, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**ab** : float, ndarray, shape :math:`\left(\mathrm{kd}+1, n\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the factor :math:`U` or :math:`L` from the Cholesky factorization :math:`A = U^\mathrm{T}U` or :math:`A = LL^\mathrm{T}`, in the same storage format as :math:`A`.
**b** : float, ndarray, shape :math:`\left(n, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04bf-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-4`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-3`)
On entry, :math:`\mathrm{kd} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{kd}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo}` not one of 'U' or 'u' or 'L' or 'l': :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq n\right)`)
The principal minor of order :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the matrix :math:`A` is not positive definite. The factorization has not been completed and the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`n+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04bf-py2-py-notes:
**Notes**
The Cholesky factorization is used to factor :math:`A` as :math:`A = U^\mathrm{T}U`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LL^\mathrm{T}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` is an upper triangular band matrix with :math:`k` superdiagonals, and :math:`L` is a lower triangular band matrix with :math:`k` subdiagonals.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04bf-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def real_posdef_tridiag_solve(d, e, b):
r"""
``real_posdef_tridiag_solve`` computes the solution to a real system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` symmetric positive definite tridiagonal matrix and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04bg-py2-py-doc:
For full information please refer to the NAG Library document for f04bg
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04bgf.html
.. _f04bg-py2-py-parameters:
**Parameters**
**d** : float, array-like, shape :math:`\left(n\right)`
Must contain the :math:`n` diagonal elements of the tridiagonal matrix :math:`A`.
**e** : float, array-like, shape :math:`\left(n-1\right)`
Must contain the :math:`\left(n-1\right)` subdiagonal elements of the tridiagonal matrix :math:`A`.
**b** : float, array-like, shape :math:`\left(n, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**d** : float, ndarray, shape :math:`\left(n\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, :math:`\mathrm{d}` is overwritten by the :math:`n` diagonal elements of the diagonal matrix :math:`D` from the :math:`LDL^\mathrm{T}` factorization of :math:`A`.
**e** : float, ndarray, shape :math:`\left(n-1\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, :math:`\mathrm{e}` is overwritten by the :math:`\left(n-1\right)` subdiagonal elements of the unit lower bidiagonal matrix :math:`L` from the :math:`LDL^\mathrm{T}` factorization of :math:`A`. (:math:`\mathrm{e}` can also be regarded as the superdiagonal of the unit upper bidiagonal factor :math:`U` from the :math:`U^\mathrm{T}DU` factorization of :math:`A`.)
**b** : float, ndarray, shape :math:`\left(n, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04bg-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-2`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq n\right)`)
The principal minor of order :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the matrix :math:`A` is not positive definite. The factorization has not been completed and the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`n+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04bg-py2-py-notes:
**Notes**
:math:`A` is factorized as :math:`A = LDL^\mathrm{T}`, where :math:`L` is a unit lower bidiagonal matrix and :math:`D` is diagonal, and the factored form of :math:`A` is then used to solve the system of equations.
.. _f04bg-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def real_symm_solve(uplo, n, a, b):
r"""
``real_symm_solve`` computes the solution to a real system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` symmetric matrix and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04bh-py2-py-doc:
For full information please refer to the NAG Library document for f04bh
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04bhf.html
.. _f04bh-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**a** : float, array-like, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
The :math:`n\times n` symmetric matrix :math:`A`.
If :math:`\mathrm{uplo} = \texttt{'U'}`, the leading :math:`\mathrm{n}` by :math:`\mathrm{n}` upper triangular part of the array :math:`\mathrm{a}` contains the upper triangular part of the matrix :math:`A`, and the strictly lower triangular part of :math:`\mathrm{a}` is not referenced.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the leading :math:`\mathrm{n}` by :math:`\mathrm{n}` lower triangular part of the array :math:`\mathrm{a}` contains the lower triangular part of the matrix :math:`A`, and the strictly upper triangular part of :math:`\mathrm{a}` is not referenced.
**b** : float, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**a** : float, ndarray, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
If no exception is raised, the block diagonal matrix :math:`D` and the multipliers used to obtain the factor :math:`U` or :math:`L` from the factorization :math:`A = UDU^\mathrm{T}` or :math:`A = LDL^\mathrm{T}` as computed by :meth:`lapacklin.dsytrf <naginterfaces.library.lapacklin.dsytrf>`.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, details of the interchanges and the block structure of :math:`D`, as determined by :meth:`lapacklin.dsytrf <naginterfaces.library.lapacklin.dsytrf>`.
:math:`\mathrm{ipiv}[k-1] > 0`
Rows and columns :math:`k` and :math:`\mathrm{ipiv}[k-1]` were interchanged, and :math:`d_{{kk}}` is a :math:`1\times 1` diagonal block.
:math:`\mathrm{uplo} = \texttt{'U'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k-2] < 0`
Rows and columns :math:`k-1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k-1:k,k-1:k}}` is a :math:`2\times 2` diagonal block.
:math:`\mathrm{uplo} = \texttt{'L'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k] < 0`
Rows and columns :math:`k+1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k:k+1,k:k+1}}` is a :math:`2\times 2` diagonal block.
**b** : float, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no constraints are violated, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04bh-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-3`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo}` not one of 'U' or 'u' or 'L' or 'l': :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq \mathrm{n}\right)`)
Diagonal block :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the block diagonal matrix is zero. The factorization has been completed, but the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04bh-py2-py-notes:
**Notes**
The diagonal pivoting method is used to factor :math:`A` as :math:`A = UDU^\mathrm{T}`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LDL^\mathrm{T}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` (or :math:`L`) is a product of permutation and unit upper (lower) triangular matrices, and :math:`D` is symmetric and block diagonal with :math:`1\times 1` and :math:`2\times 2` diagonal blocks.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04bh-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def real_symm_packed_solve(uplo, n, ap, b):
r"""
``real_symm_packed_solve`` computes the solution to a real system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` symmetric matrix, stored in packed format and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04bj-py2-py-doc:
For full information please refer to the NAG Library document for f04bj
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04bjf.html
.. _f04bj-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**ap** : float, array-like, shape :math:`\left(\mathrm{n}\times \left(\mathrm{n}+1\right)/2\right)`
The :math:`n\times n` symmetric matrix :math:`A`, packed column-wise in a linear array. The :math:`j`\ th column of the matrix :math:`A` is stored in the array :math:`\mathrm{ap}` as follows:
**b** : float, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**ap** : float, ndarray, shape :math:`\left(\mathrm{n}\times \left(\mathrm{n}+1\right)/2\right)`
If no exception is raised, the block diagonal matrix :math:`D` and the multipliers used to obtain the factor :math:`U` or :math:`L` from the factorization :math:`A = UDU^\mathrm{T}` or :math:`A = LDL^\mathrm{T}` as computed by :meth:`lapacklin.dsptrf <naginterfaces.library.lapacklin.dsptrf>`, stored as a packed triangular matrix in the same storage format as :math:`A`.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, details of the interchanges and the block structure of :math:`D`, as determined by :meth:`lapacklin.dsptrf <naginterfaces.library.lapacklin.dsptrf>`.
If :math:`\mathrm{ipiv}[k-1] > 0`, then rows and columns :math:`k` and :math:`\mathrm{ipiv}[k-1]` were interchanged, and :math:`d_{{kk}}` is a :math:`1\times 1` diagonal block;
if :math:`\mathrm{uplo} = \texttt{'U'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k-2] < 0`, then rows and columns :math:`k-1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k-1:k,k-1:k}}` is a :math:`2\times 2` diagonal block;
if :math:`\mathrm{uplo} = \texttt{'L'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k] < 0`, then rows and columns :math:`k+1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k:k+1,k:k+1}}` is a :math:`2\times 2` diagonal block.
**b** : float, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no constraints are violated, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04bj-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-3`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo}` not one of 'U' or 'u' or 'L' or 'l': :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq \mathrm{n}\right)`)
Diagonal block :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the block diagonal matrix is zero. The factorization has been completed, but the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04bj-py2-py-notes:
**Notes**
The diagonal pivoting method is used to factor :math:`A` as :math:`A = UDU^\mathrm{T}`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LDL^\mathrm{T}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` (or :math:`L`) is a product of permutation and unit upper (lower) triangular matrices, and :math:`D` is symmetric and block diagonal with :math:`1\times 1` and :math:`2\times 2` diagonal blocks.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04bj-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_square_solve(n, a, b):
r"""
``complex_square_solve`` computes the solution to a complex system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` matrix and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04ca-py2-py-doc:
For full information please refer to the NAG Library document for f04ca
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04caf.html
.. _f04ca-py2-py-parameters:
**Parameters**
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**a** : complex, array-like, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
The :math:`n\times n` coefficient matrix :math:`A`.
**b** : complex, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**a** : complex, ndarray, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
If no exception is raised, the factors :math:`L` and :math:`U` from the factorization :math:`A = PLU`. The unit diagonal elements of :math:`L` are not stored.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, the pivot indices that define the permutation matrix :math:`P`; at the :math:`i`\ th step row :math:`i` of the matrix was interchanged with row :math:`\mathrm{ipiv}[i-1]`. :math:`\mathrm{ipiv}[i-1] = i` indicates a row interchange was not required.
**b** : complex, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no constraints are violated, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04ca-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-2`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq \mathrm{n}\right)`)
Diagonal element :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the upper triangular factor is zero. The factorization has been completed, but the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04ca-py2-py-notes:
**Notes**
The :math:`LU` decomposition with partial pivoting and row interchanges is used to factor :math:`A` as :math:`A = PLU`, where :math:`P` is a permutation matrix, :math:`L` is unit lower triangular, and :math:`U` is upper triangular.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04ca-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_band_solve(n, kl, ku, ab, b):
r"""
``complex_band_solve`` computes the solution to a complex system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` band matrix, with :math:`k_l` subdiagonals and :math:`k_u` superdiagonals, and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04cb-py2-py-doc:
For full information please refer to the NAG Library document for f04cb
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04cbf.html
.. _f04cb-py2-py-parameters:
**Parameters**
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**kl** : int
The number of subdiagonals :math:`k_l`, within the band of :math:`A`.
**ku** : int
The number of superdiagonals :math:`k_u`, within the band of :math:`A`.
**ab** : complex, array-like, shape :math:`\left(2\times \mathrm{kl}+\mathrm{ku}+1, \mathrm{n}\right)`
The :math:`n\times n` matrix :math:`A`.
See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04cbf.html#fcomments>`__ for further details.
**b** : complex, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**ab** : complex, ndarray, shape :math:`\left(2\times \mathrm{kl}+\mathrm{ku}+1, \mathrm{n}\right)`
If no exception is raised, :math:`\mathrm{ab}` is overwritten by details of the factorization.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, the pivot indices that define the permutation matrix :math:`P`; at the :math:`i`\ th step row :math:`i` of the matrix was interchanged with row :math:`\mathrm{ipiv}[i-1]`. :math:`\mathrm{ipiv}[i-1] = i` indicates a row interchange was not required.
**b** : complex, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no exception is raised, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = \left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04cb-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-6`)
On entry, :math:`\textit{ldab} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{kl} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ku} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ldab}\geq 2\times \mathrm{kl}+\mathrm{ku}+1`.
(`errno` :math:`-4`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-3`)
On entry, :math:`\mathrm{ku} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ku}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{kl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{kl}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq \mathrm{n}\right)`)
Diagonal element :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the upper triangular factor is zero. The factorization has been completed, but the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04cb-py2-py-notes:
**Notes**
The :math:`LU` decomposition with partial pivoting and row interchanges is used to factor :math:`A` as :math:`A = PLU`, where :math:`P` is a permutation matrix, :math:`L` is the product of permutation matrices and unit lower triangular matrices with :math:`k_l` subdiagonals, and :math:`U` is upper triangular with :math:`\left(k_l+k_u\right)` superdiagonals.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04cb-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_tridiag_solve(n, dl, d, du, b):
r"""
``complex_tridiag_solve`` computes the solution to a complex system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` tridiagonal matrix and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04cc-py2-py-doc:
For full information please refer to the NAG Library document for f04cc
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04ccf.html
.. _f04cc-py2-py-parameters:
**Parameters**
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**dl** : complex, array-like, shape :math:`\left(\mathrm{n}-1\right)`
Must contain the :math:`\left(n-1\right)` subdiagonal elements of the matrix :math:`A`.
**d** : complex, array-like, shape :math:`\left(\mathrm{n}\right)`
Must contain the :math:`n` diagonal elements of the matrix :math:`A`.
**du** : complex, array-like, shape :math:`\left(\mathrm{n}-1\right)`
Must contain the :math:`\left(n-1\right)` superdiagonal elements of the matrix :math:`A`
**b** : complex, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**dl** : complex, ndarray, shape :math:`\left(\mathrm{n}-1\right)`
If no exception is raised, :math:`\mathrm{dl}` is overwritten by the :math:`\left(n-1\right)` multipliers that define the matrix :math:`L` from the :math:`LU` factorization of :math:`A`.
**d** : complex, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, :math:`\mathrm{d}` is overwritten by the :math:`n` diagonal elements of the upper triangular matrix :math:`U` from the :math:`LU` factorization of :math:`A`.
**du** : complex, ndarray, shape :math:`\left(\mathrm{n}-1\right)`
If no exception is raised, :math:`\mathrm{du}` is overwritten by the :math:`\left(n-1\right)` elements of the first superdiagonal of :math:`U`.
**du2** : complex, ndarray, shape :math:`\left(\mathrm{n}-2\right)`
If no exception is raised, :math:`\mathrm{du2}` returns the :math:`\left(n-2\right)` elements of the second superdiagonal of :math:`U`.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, the pivot indices that define the permutation matrix :math:`P`; at the :math:`i`\ th step row :math:`i` of the matrix was interchanged with row :math:`\mathrm{ipiv}[i-1]`. :math:`\mathrm{ipiv}[i-1]` will always be either :math:`i` or :math:`\left(i+1\right)`; :math:`\mathrm{ipiv}[i-1] = i` indicates a row interchange was not required.
**b** : complex, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no constraints are violated, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04cc-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-2`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq \mathrm{n}\right)`)
Diagonal element :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the upper triangular factor is zero. The factorization has been completed, but the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04cc-py2-py-notes:
**Notes**
The :math:`LU` decomposition with partial pivoting and row interchanges is used to factor :math:`A` as :math:`A = PLU`, where :math:`P` is a permutation matrix, :math:`L` is unit lower triangular with at most one nonzero subdiagonal element, and :math:`U` is an upper triangular band matrix with two superdiagonals.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
Note that the equations :math:`A^\mathrm{T}X = B` may be solved by interchanging the order of the arguments :math:`\mathrm{du}` and :math:`\mathrm{dl}`.
.. _f04cc-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_posdef_solve(uplo, a, b):
r"""
``complex_posdef_solve`` computes the solution to a complex system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` Hermitian positive definite matrix and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04cd-py2-py-doc:
For full information please refer to the NAG Library document for f04cd
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04cdf.html
.. _f04cd-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**a** : complex, array-like, shape :math:`\left(n, n\right)`
The :math:`n\times n` Hermitian matrix :math:`A`.
If :math:`\mathrm{uplo} = \texttt{'U'}`, the leading :math:`\textit{n}` by :math:`\textit{n}` upper triangular part of :math:`\mathrm{a}` contains the upper triangular part of the matrix :math:`A`, and the strictly lower triangular part of :math:`\mathrm{a}` is not referenced.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the leading :math:`\textit{n}` by :math:`\textit{n}` lower triangular part of :math:`\mathrm{a}` contains the lower triangular part of the matrix :math:`A`, and the strictly upper triangular part of :math:`\mathrm{a}` is not referenced.
**b** : complex, array-like, shape :math:`\left(n, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**a** : complex, ndarray, shape :math:`\left(n, n\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the factor :math:`U` or :math:`L` from the Cholesky factorization :math:`A = U^\mathrm{H}U` or :math:`A = LL^\mathrm{H}`.
**b** : complex, ndarray, shape :math:`\left(n, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04cd-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-3`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo}` not one of 'U' or 'u' or 'L' or 'l': :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq n\right)`)
The principal minor of order :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the matrix :math:`A` is not positive definite. The factorization has not been completed and the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`n+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04cd-py2-py-notes:
**Notes**
The Cholesky factorization is used to factor :math:`A` as :math:`A = U^\mathrm{H}U`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LL^\mathrm{H}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` is an upper triangular matrix and :math:`L` is a lower triangular matrix.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04cd-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_posdef_packed_solve(uplo, ap, b):
r"""
``complex_posdef_packed_solve`` computes the solution to a complex system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` Hermitian positive definite matrix, stored in packed format, and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04ce-py2-py-doc:
For full information please refer to the NAG Library document for f04ce
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04cef.html
.. _f04ce-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**ap** : complex, array-like, shape :math:`\left(n\times \left(n+1\right)/2\right)`
The :math:`n\times n` Hermitian matrix :math:`A`. The upper or lower triangular part of the Hermitian matrix is packed column-wise in a linear array. The :math:`j`\ th column of :math:`A` is stored in the array :math:`\mathrm{ap}` as follows:
**b** : complex, array-like, shape :math:`\left(n, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**ap** : complex, ndarray, shape :math:`\left(n\times \left(n+1\right)/2\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the factor :math:`U` or :math:`L` from the Cholesky factorization :math:`A = U^\mathrm{H}U` or :math:`A = LL^\mathrm{H}`, in the same storage format as :math:`A`.
**b** : complex, ndarray, shape :math:`\left(n, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04ce-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-3`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo} \neq \texttt{'U'}` or :math:`\texttt{'L'}`: :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq n\right)`)
The principal minor of order :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the matrix :math:`A` is not positive definite. The factorization has not been completed and the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`n+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04ce-py2-py-notes:
**Notes**
The Cholesky factorization is used to factor :math:`A` as :math:`A = U^\mathrm{H}U`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LL^\mathrm{H}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` is an upper triangular matrix and :math:`L` is a lower triangular matrix.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04ce-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_posdef_band_solve(uplo, kd, ab, b):
r"""
``complex_posdef_band_solve`` computes the solution to a complex system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` Hermitian positive definite band matrix of band width :math:`2k+1`, and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04cf-py2-py-doc:
For full information please refer to the NAG Library document for f04cf
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04cff.html
.. _f04cf-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**kd** : int
The number of superdiagonals :math:`k` (and the number of subdiagonals) of the band matrix :math:`A`.
**ab** : complex, array-like, shape :math:`\left(\mathrm{kd}+1, n\right)`
The :math:`n\times n` Hermitian band matrix :math:`A`. The upper or lower triangular part of the Hermitian matrix is stored in the first :math:`\mathrm{kd}+1` rows of the array. The :math:`j`\ th column of :math:`A` is stored in the :math:`j`\ th column of the array :math:`\mathrm{ab}` as follows:
See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04cff.html#fcomments>`__ below for further details.
**b** : complex, array-like, shape :math:`\left(n, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**ab** : complex, ndarray, shape :math:`\left(\mathrm{kd}+1, n\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the factor :math:`U` or :math:`L` from the Cholesky factorization :math:`A = U^\mathrm{H}U` or :math:`A = LL^\mathrm{H}`, in the same storage format as :math:`A`.
**b** : complex, ndarray, shape :math:`\left(n, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04cf-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-4`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-3`)
On entry, :math:`\mathrm{kd} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{kd}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo}` not one of 'U' or 'u' or 'L' or 'l': :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq n\right)`)
The principal minor of order :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the matrix :math:`A` is not positive definite. The factorization has not been completed and the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`n+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04cf-py2-py-notes:
**Notes**
The Cholesky factorization is used to factor :math:`A` as :math:`A = U^\mathrm{H}U`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LL^\mathrm{H}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` is an upper triangular band matrix with :math:`k` superdiagonals, and :math:`L` is a lower triangular band matrix with :math:`k` subdiagonals.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04cf-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_posdef_tridiag_solve(d, e, b):
r"""
``complex_posdef_tridiag_solve`` computes the solution to a complex system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` Hermitian positive definite tridiagonal matrix and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04cg-py2-py-doc:
For full information please refer to the NAG Library document for f04cg
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04cgf.html
.. _f04cg-py2-py-parameters:
**Parameters**
**d** : float, array-like, shape :math:`\left(n\right)`
Must contain the :math:`n` diagonal elements of the tridiagonal matrix :math:`A`.
**e** : complex, array-like, shape :math:`\left(n-1\right)`
Must contain the :math:`\left(n-1\right)` subdiagonal elements of the tridiagonal matrix :math:`A`.
**b** : complex, array-like, shape :math:`\left(n, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**d** : float, ndarray, shape :math:`\left(n\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, :math:`\mathrm{d}` is overwritten by the :math:`n` diagonal elements of the diagonal matrix :math:`D` from the :math:`LDL^\mathrm{H}` factorization of :math:`A`.
**e** : complex, ndarray, shape :math:`\left(n-1\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, :math:`\mathrm{e}` is overwritten by the :math:`\left(n-1\right)` subdiagonal elements of the unit lower bidiagonal matrix :math:`L` from the :math:`LDL^\mathrm{H}` factorization of :math:`A`. (:math:`\mathrm{e}` can also be regarded as the conjugate of the superdiagonal of the unit upper bidiagonal factor :math:`U` from the :math:`U^\mathrm{H}DU` factorization of :math:`A`.)
**b** : complex, ndarray, shape :math:`\left(n, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1, \left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04cg-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-2`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq n\right)`)
The principal minor of order :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the matrix :math:`A` is not positive definite. The factorization has not been completed and the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`n+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04cg-py2-py-notes:
**Notes**
:math:`A` is factorized as :math:`A = LDL^\mathrm{H}`, where :math:`L` is a unit lower bidiagonal matrix and :math:`D` is a real diagonal matrix, and the factored form of :math:`A` is then used to solve the system of equations.
.. _f04cg-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_herm_solve(uplo, n, a, b):
r"""
``complex_herm_solve`` computes the solution to a complex system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` Hermitian matrix and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04ch-py2-py-doc:
For full information please refer to the NAG Library document for f04ch
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04chf.html
.. _f04ch-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**a** : complex, array-like, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
The :math:`n\times n` Hermitian matrix :math:`A`.
If :math:`\mathrm{uplo} = \texttt{'U'}`, the leading :math:`\mathrm{n}` by :math:`\mathrm{n}` upper triangular part of the array :math:`\mathrm{a}` contains the upper triangular part of the matrix :math:`A`, and the strictly lower triangular part of :math:`\mathrm{a}` is not referenced.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the leading :math:`\mathrm{n}` by :math:`\mathrm{n}` lower triangular part of the array :math:`\mathrm{a}` contains the lower triangular part of the matrix :math:`A`, and the strictly upper triangular part of :math:`\mathrm{a}` is not referenced.
**b** : complex, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**a** : complex, ndarray, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
If no exception is raised, the block diagonal matrix :math:`D` and the multipliers used to obtain the factor :math:`U` or :math:`L` from the factorization :math:`A = UDU^\mathrm{H}` or :math:`A = LDL^\mathrm{H}` as computed by :meth:`lapacklin.zhetrf <naginterfaces.library.lapacklin.zhetrf>`.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, details of the interchanges and the block structure of :math:`D`, as determined by :meth:`lapacklin.zhetrf <naginterfaces.library.lapacklin.zhetrf>`.
If :math:`\mathrm{ipiv}[k-1] > 0`, then rows and columns :math:`k` and :math:`\mathrm{ipiv}[k-1]` were interchanged, and :math:`d_{{kk}}` is a :math:`1\times 1` diagonal block;
if :math:`\mathrm{uplo} = \texttt{'U'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k-2] < 0`, then rows and columns :math:`k-1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k-1:k,k-1:k}}` is a :math:`2\times 2` diagonal block;
if :math:`\mathrm{uplo} = \texttt{'L'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k] < 0`, then rows and columns :math:`k+1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k:k+1,k:k+1}}` is a :math:`2\times 2` diagonal block.
**b** : complex, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no constraints are violated, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04ch-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-3`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo}` not one of 'U' or 'u' or 'L' or 'l': :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq \mathrm{n}\right)`)
Diagonal block :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the block diagonal matrix is zero. The factorization has been completed, but the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04ch-py2-py-notes:
**Notes**
The diagonal pivoting method is used to factor :math:`A` as :math:`A = UDU^\mathrm{H}`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LDL^\mathrm{H}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` (or :math:`L`) is a product of permutation and unit upper (lower) triangular matrices, and :math:`D` is Hermitian and block diagonal with :math:`1\times 1` and :math:`2\times 2` diagonal blocks.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04ch-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_herm_packed_solve(uplo, n, ap, b):
r"""
``complex_herm_packed_solve`` computes the solution to a complex system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` complex Hermitian matrix, stored in packed format and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04cj-py2-py-doc:
For full information please refer to the NAG Library document for f04cj
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04cjf.html
.. _f04cj-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**ap** : complex, array-like, shape :math:`\left(\mathrm{n}\times \left(\mathrm{n}+1\right)/2\right)`
The :math:`n\times n` Hermitian matrix :math:`A`, packed column-wise in a linear array. The :math:`j`\ th column of the matrix :math:`A` is stored in the array :math:`\mathrm{ap}` as follows:
**b** : complex, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**ap** : complex, ndarray, shape :math:`\left(\mathrm{n}\times \left(\mathrm{n}+1\right)/2\right)`
If no exception is raised, the block diagonal matrix :math:`D` and the multipliers used to obtain the factor :math:`U` or :math:`L` from the factorization :math:`A = UDU^\mathrm{H}` or :math:`A = LDL^\mathrm{H}` as computed by :meth:`lapacklin.zhptrf <naginterfaces.library.lapacklin.zhptrf>`, stored as a packed triangular matrix in the same storage format as :math:`A`.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, details of the interchanges and the block structure of :math:`D`, as determined by :meth:`lapacklin.zhptrf <naginterfaces.library.lapacklin.zhptrf>`.
If :math:`\mathrm{ipiv}[k-1] > 0`, then rows and columns :math:`k` and :math:`\mathrm{ipiv}[k-1]` were interchanged, and :math:`d_{{kk}}` is a :math:`1\times 1` diagonal block;
if :math:`\mathrm{uplo} = \texttt{'U'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k-2] < 0`, then rows and columns :math:`k-1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k-1:k,k-1:k}}` is a :math:`2\times 2` diagonal block;
if :math:`\mathrm{uplo} = \texttt{'L'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k] < 0`, then rows and columns :math:`k+1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k:k+1,k:k+1}}` is a :math:`2\times 2` diagonal block.
**b** : complex, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no constraints are violated, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04cj-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-3`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo} \neq \texttt{'U'}` or :math:`\texttt{'L'}`: :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq \mathrm{n}\right)`)
Diagonal block :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the block diagonal matrix is zero. The factorization has been completed, but the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04cj-py2-py-notes:
**Notes**
The diagonal pivoting method is used to factor :math:`A` as :math:`A = UDU^\mathrm{H}`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LDL^\mathrm{H}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` (or :math:`L`) is a product of permutation and unit upper (lower) triangular matrices, and :math:`D` is Hermitian and block diagonal with :math:`1\times 1` and :math:`2\times 2` diagonal blocks.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04cj-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_symm_solve(uplo, n, a, b):
r"""
``complex_symm_solve`` computes the solution to a complex system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` complex symmetric matrix and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04dh-py2-py-doc:
For full information please refer to the NAG Library document for f04dh
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04dhf.html
.. _f04dh-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**a** : complex, array-like, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
The :math:`n\times n` complex symmetric matrix :math:`A`.
If :math:`\mathrm{uplo} = \texttt{'U'}`, the leading :math:`\mathrm{n}` by :math:`\mathrm{n}` upper triangular part of the array :math:`\mathrm{a}` contains the upper triangular part of the matrix :math:`A`, and the strictly lower triangular part of :math:`\mathrm{a}` is not referenced.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the leading :math:`\mathrm{n}` by :math:`\mathrm{n}` lower triangular part of the array :math:`\mathrm{a}` contains the lower triangular part of the matrix :math:`A`, and the strictly upper triangular part of :math:`\mathrm{a}` is not referenced.
**b** : complex, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**a** : complex, ndarray, shape :math:`\left(\mathrm{n}, \mathrm{n}\right)`
If no exception is raised, the block diagonal matrix :math:`D` and the multipliers used to obtain the factor :math:`U` or :math:`L` from the factorization :math:`A = UDU^\mathrm{T}` or :math:`A = LDL^\mathrm{T}` as computed by :meth:`lapacklin.zsytrf <naginterfaces.library.lapacklin.zsytrf>`.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no exception is raised, details of the interchanges and the block structure of :math:`D`, as determined by :meth:`lapacklin.zsytrf <naginterfaces.library.lapacklin.zsytrf>`.
If :math:`\mathrm{ipiv}[k-1] > 0`, then rows and columns :math:`k` and :math:`\mathrm{ipiv}[k-1]` were interchanged, and :math:`d_{{kk}}` is a :math:`1\times 1` diagonal block;
if :math:`\mathrm{uplo} = \texttt{'U'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k-2] < 0`, then rows and columns :math:`k-1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k-1:k,k-1:k}}` is a :math:`2\times 2` diagonal block;
if :math:`\mathrm{uplo} = \texttt{'L'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k] < 0`, then rows and columns :math:`k+1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k:k+1,k:k+1}}` is a :math:`2\times 2` diagonal block.
**b** : complex, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no constraints are violated, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04dh-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-3`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo}` not one of 'U' or 'u' or 'L' or 'l': :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq \mathrm{n}\right)`)
Diagonal block :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the block diagonal matrix is zero. The factorization has been completed, but the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04dh-py2-py-notes:
**Notes**
The diagonal pivoting method is used to factor :math:`A` as :math:`A = UDU^\mathrm{T}`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LDL^\mathrm{T}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` (or :math:`L`) is a product of permutation and unit upper (lower) triangular matrices, and :math:`D` is symmetric and block diagonal with :math:`1\times 1` and :math:`2\times 2` diagonal blocks.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04dh-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def complex_symm_packed_solve(uplo, n, ap, b):
r"""
``complex_symm_packed_solve`` computes the solution to a complex system of linear equations :math:`AX = B`, where :math:`A` is an :math:`n\times n` complex symmetric matrix, stored in packed format and :math:`X` and :math:`B` are :math:`n\times r` matrices.
An estimate of the condition number of :math:`A` and an error bound for the computed solution are also returned.
.. _f04dj-py2-py-doc:
For full information please refer to the NAG Library document for f04dj
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04djf.html
.. _f04dj-py2-py-parameters:
**Parameters**
**uplo** : str, length 1
If :math:`\mathrm{uplo} = \texttt{'U'}`, the upper triangle of the matrix :math:`A` is stored.
If :math:`\mathrm{uplo} = \texttt{'L'}`, the lower triangle of the matrix :math:`A` is stored.
**n** : int
The number of linear equations :math:`n`, i.e., the order of the matrix :math:`A`.
**ap** : complex, array-like, shape :math:`\left(\mathrm{n}\times \left(\mathrm{n}+1\right)/2\right)`
The :math:`n\times n` symmetric matrix :math:`A`, packed column-wise in a linear array. The :math:`j`\ th column of the matrix :math:`A` is stored in the array :math:`\mathrm{ap}` as follows:
**b** : complex, array-like, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
The :math:`n\times r` matrix of right-hand sides :math:`B`.
**Returns**
**ap** : complex, ndarray, shape :math:`\left(\mathrm{n}\times \left(\mathrm{n}+1\right)/2\right)`
If no exception is raised, the block diagonal matrix :math:`D` and the multipliers used to obtain the factor :math:`U` or :math:`L` from the factorization :math:`A = UDU^\mathrm{H}` or :math:`A = LDL^\mathrm{H}` as computed by :meth:`lapacklin.zsptrf <naginterfaces.library.lapacklin.zsptrf>`, stored as a packed triangular matrix in the same storage format as :math:`A`.
**ipiv** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
If no constraints are violated, details of the interchanges and the block structure of :math:`D`, as determined by :meth:`lapacklin.zsptrf <naginterfaces.library.lapacklin.zsptrf>`.
If :math:`\mathrm{ipiv}[k-1] > 0`, then rows and columns :math:`k` and :math:`\mathrm{ipiv}[k-1]` were interchanged, and :math:`d_{{kk}}` is a :math:`1\times 1` diagonal block;
if :math:`\mathrm{uplo} = \texttt{'U'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k-2] < 0`, then rows and columns :math:`k-1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k-1:k,k-1:k}}` is a :math:`2\times 2` diagonal block;
if :math:`\mathrm{uplo} = \texttt{'L'}` and :math:`\mathrm{ipiv}[k-1] = \mathrm{ipiv}[k] < 0`, then rows and columns :math:`k+1` and :math:`-\mathrm{ipiv}[k-1]` were interchanged and :math:`d_{{k:k+1,k:k+1}}` is a :math:`2\times 2` diagonal block.
**b** : complex, ndarray, shape :math:`\left(\mathrm{n}, \textit{nrhs}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, the :math:`n\times r` solution matrix :math:`X`.
**rcond** : float
If no constraints are violated, an estimate of the reciprocal of the condition number of the matrix :math:`A`, computed as :math:`\mathrm{rcond} = 1/\left(\left\lVert A\right\rVert_1\left\lVert A^{-1}\right\rVert_1\right)`.
**errbnd** : float
If the function exits successfully or :math:`\mathrm{errno}` = :math:`\mathrm{n}` + 1, an estimate of the forward error bound for a computed solution :math:`\hat{x}`, such that :math:`\left\lVert {\hat{x}-x}\right\rVert_1/\left\lVert x\right\rVert_1\leq \mathrm{errbnd}`, where :math:`\hat{x}` is a column of the computed solution returned in the array :math:`\mathrm{b}` and :math:`x` is the corresponding column of the exact solution :math:`X`. If :math:`\mathrm{rcond}` is less than machine precision, :math:`\mathrm{errbnd}` is returned as unity.
.. _f04dj-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-3`)
On entry, :math:`\textit{nrhs} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nrhs}\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{uplo} \neq \texttt{'U'}` or :math:`\texttt{'L'}`: :math:`\mathrm{uplo} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`\left(i > 0\right)\text{ and }\left(i \leq \mathrm{n}\right)`)
Diagonal block :math:`\langle\mathit{\boldsymbol{value}}\rangle` of the block diagonal matrix is zero. The factorization has been completed, but the solution could not be computed.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`\mathrm{n}+1`)
A solution has been computed, but :math:`\mathrm{rcond}` is less than machine precision so that the matrix :math:`A` is numerically singular.
.. _f04dj-py2-py-notes:
**Notes**
The diagonal pivoting method is used to factor :math:`A` as :math:`A = UDU^\mathrm{T}`, if :math:`\mathrm{uplo} = \texttt{'U'}`, or :math:`A = LDL^\mathrm{T}`, if :math:`\mathrm{uplo} = \texttt{'L'}`, where :math:`U` (or :math:`L`) is a product of permutation and unit upper (lower) triangular matrices, and :math:`D` is symmetric and block diagonal with :math:`1\times 1` and :math:`2\times 2` diagonal blocks.
The factored form of :math:`A` is then used to solve the system of equations :math:`AX = B`.
.. _f04dj-py2-py-references:
**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, `LAPACK Users' Guide`, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug
Higham, N J, 2002, `Accuracy and Stability of Numerical Algorithms`, (2nd Edition), SIAM, Philadelphia
"""
raise NotImplementedError
[docs]def real_toeplitz_yule(t, wantp, wantv):
r"""
``real_toeplitz_yule`` solves the Yule--Walker equations for a real symmetric positive definite Toeplitz system.
.. _f04fe-py2-py-doc:
For full information please refer to the NAG Library document for f04fe
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04fef.html
.. _f04fe-py2-py-parameters:
**Parameters**
**t** : float, array-like, shape :math:`\left(n+1\right)`
:math:`\mathrm{t}[0]` must contain the value :math:`\tau_0` of the diagonal elements of :math:`T`, and the remaining :math:`\textit{n}` elements of :math:`\mathrm{t}` must contain the elements of the vector :math:`t`.
**wantp** : bool
Must be set to :math:`\mathbf{True}` if the partial (auto)correlation coefficients are required, and must be set to :math:`\mathbf{False}` otherwise.
**wantv** : bool
Must be set to :math:`\mathbf{True}` if the mean square prediction errors are required, and must be set to :math:`\mathbf{False}` otherwise.
**Returns**
**x** : float, ndarray, shape :math:`\left(n\right)`
The solution vector :math:`x`.
**p** : float, ndarray, shape :math:`\left(:\right)`
With :math:`\mathrm{wantp}` as :math:`\mathbf{True}`, the :math:`i`\ th element of :math:`\mathrm{p}` contains the partial (auto)correlation coefficient, or reflection coefficient, :math:`p_i` for the :math:`i`\ th step. (See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04fef.html#fcomments>`__ and submodule :mod:`~naginterfaces.library.tsa`.) If :math:`\mathrm{wantp}` is :math:`\mathbf{False}`, :math:`\mathrm{p}` is not referenced. Note that in any case, :math:`x_n = p_n`.
**v** : float, ndarray, shape :math:`\left(:\right)`
With :math:`\mathrm{wantv}` as :math:`\mathbf{True}`, the :math:`i`\ th element of :math:`\mathrm{v}` contains the mean square prediction error, or predictor error variance ratio, :math:`v_i`, for the :math:`i`\ th step. (See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04fef.html#fcomments>`__ and submodule :mod:`~naginterfaces.library.tsa`.) If :math:`\mathrm{wantv}` is :math:`\mathbf{False}`, :math:`\mathrm{v}` is not referenced.
**vlast** : float
The value of :math:`v_n`, the mean square prediction error for the final step.
.. _f04fe-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-1`)
On entry, :math:`\mathrm{t}[0] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{t}[0] > 0.0`.
(`errno` :math:`-1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`i > 0`)
Principal minor :math:`\langle\mathit{\boldsymbol{value}}\rangle` is not positive definite. Value of the reflection coefficient is :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
.. _f04fe-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``real_toeplitz_yule`` solves the equations
.. math::
Tx = -t\text{,}
where :math:`T` is the :math:`n\times n` symmetric positive definite Toeplitz matrix
.. math::
T = \begin{pmatrix}\tau_0&\tau_1&\tau_2&\ldots &\tau_{{n-1}}\\\tau_1&\tau_0&\tau_1&\ldots &\tau_{{n-2}}\\\tau_2&\tau_1&\tau_0&\ldots &\tau_{{n-3}}\\.&.&.&&.\\\tau_{{n-1}}&\tau_{{n-2}}&\tau_{{n-3}}&\ldots &\tau_0\end{pmatrix}
and :math:`t` is the vector
.. math::
t^\mathrm{T} = \left(\tau_1, {\tau_2\ldots \tau_n}\right)\text{.}
The function uses the method of Durbin (see Durbin (1960) and Golub and Van Loan (1996)).
Optionally the mean square prediction errors and/or the partial correlation coefficients for each step can be returned.
.. _f04fe-py2-py-references:
**References**
Bunch, J R, 1985, `Stability of methods for solving Toeplitz systems of equations`, SIAM J. Sci. Statist. Comput. (6), 349--364
Bunch, J R, 1987, `The weak and strong stability of algorithms in numerical linear algebra`, Linear Algebra Appl. (88/89), 49--66
Cybenko, G, 1980, `The numerical stability of the Levinson--Durbin algorithm for Toeplitz systems of equations`, SIAM J. Sci. Statist. Comput. (1), 303--319
Durbin, J, 1960, `The fitting of time series models`, Rev. Inst. Internat. Stat. (28), 233
Golub, G H and Van Loan, C F, 1996, `Matrix Computations`, (3rd Edition), Johns Hopkins University Press, Baltimore
"""
raise NotImplementedError
[docs]def real_toeplitz_solve(n, t, b, wantp):
r"""
``real_toeplitz_solve`` solves the equations :math:`Tx = b`, where :math:`T` is a real symmetric positive definite Toeplitz matrix.
.. _f04ff-py2-py-doc:
For full information please refer to the NAG Library document for f04ff
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04fff.html
.. _f04ff-py2-py-parameters:
**Parameters**
**n** : int
The order of the Toeplitz matrix :math:`T`.
**t** : float, array-like, shape :math:`\left(\mathrm{n}\right)`
:math:`\mathrm{t}[\textit{i}]` must contain the value :math:`\tau_{\textit{i}}`, for :math:`\textit{i} = 0,1,\ldots,\mathrm{n}-1`.
**b** : float, array-like, shape :math:`\left(\mathrm{n}\right)`
The right-hand side vector :math:`b`.
**wantp** : bool
Must be set to :math:`\mathbf{True}` if the reflection coefficients are required, and must be set to :math:`\mathbf{False}` otherwise.
**Returns**
**x** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
The solution vector :math:`x`.
**p** : float, ndarray, shape :math:`\left(:\right)`
With :math:`\mathrm{wantp}` as :math:`\mathbf{True}`, the :math:`i`\ th element of :math:`\mathrm{p}` contains the reflection coefficient, :math:`p_{\textit{i}}`, for the :math:`\textit{i}`\ th step, for :math:`\textit{i} = 1,2,\ldots,\mathrm{n}-1`. (See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04fff.html#fcomments>`__.) If :math:`\mathrm{wantp}` is :math:`\mathbf{False}`, :math:`\mathrm{p}` is not referenced.
.. _f04ff-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-1`)
On entry, :math:`\mathrm{t}[0] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{t}[0] > 0.0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`i > 0`)
Principal minor :math:`\langle\mathit{\boldsymbol{value}}\rangle` is not positive definite. Value of the reflection coefficient is :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
.. _f04ff-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``real_toeplitz_solve`` solves the equations
.. math::
Tx = b\text{,}
where :math:`T` is the :math:`n\times n` symmetric positive definite Toeplitz matrix
.. math::
T = \begin{pmatrix}\tau_0&\tau_1&\tau_2&\ldots &\tau_{{n-1}}\\\tau_1&\tau_0&\tau_1&\ldots &\tau_{{n-2}}\\\tau_2&\tau_1&\tau_0&\ldots &\tau_{{n-3}}\\.&.&.&&.\\\tau_{{n-1}}&\tau_{{n-2}}&\tau_{{n-3}}&\ldots &\tau_0\end{pmatrix}
and :math:`b` is an :math:`n`-element vector.
The function uses the method of Levinson (see Levinson (1947) and Golub and Van Loan (1996)).
Optionally, the reflection coefficients for each step may also be returned.
.. _f04ff-py2-py-references:
**References**
Bunch, J R, 1985, `Stability of methods for solving Toeplitz systems of equations`, SIAM J. Sci. Statist. Comput. (6), 349--364
Bunch, J R, 1987, `The weak and strong stability of algorithms in numerical linear algebra`, Linear Algebra Appl. (88/89), 49--66
Cybenko, G, 1980, `The numerical stability of the Levinson--Durbin algorithm for Toeplitz systems of equations`, SIAM J. Sci. Statist. Comput. (1), 303--319
Golub, G H and Van Loan, C F, 1996, `Matrix Computations`, (3rd Edition), Johns Hopkins University Press, Baltimore
Levinson, N, 1947, `The Weiner RMS error criterion in filter design and prediction`, J. Math. Phys. (25), 261--278
"""
raise NotImplementedError
[docs]def real_gen_solve(a, b, tol, lwork):
r"""
``real_gen_solve`` finds the solution of a linear least squares problem, :math:`Ax = b`, where :math:`A` is a real :math:`m\times n\left(m\geq n\right)` matrix and :math:`b` is an :math:`m` element vector.
If the matrix of observations is not of full rank, then the minimal least squares solution is returned.
.. _f04jg-py2-py-doc:
For full information please refer to the NAG Library document for f04jg
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04jgf.html
.. _f04jg-py2-py-parameters:
**Parameters**
**a** : float, array-like, shape :math:`\left(m, n\right)`
The :math:`m\times n` matrix :math:`A`.
**b** : float, array-like, shape :math:`\left(m\right)`
The right-hand side vector :math:`b`.
**tol** : float
A relative tolerance to be used to determine the rank of :math:`A`. :math:`\mathrm{tol}` should be chosen as approximately the largest relative error in the elements of :math:`A`. For example, if the elements of :math:`A` are correct to about :math:`4` significant figures then :math:`\mathrm{tol}` should be set to about :math:`5\times 10^{-4}`. See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04jgf.html#fcomments>`__ for a description of how :math:`\mathrm{tol}` is used to determine rank. If :math:`\mathrm{tol}` is outside the range :math:`\left(\epsilon, 1.0\right)`, where :math:`\epsilon` is the machine precision, the value :math:`\epsilon` is used in place of :math:`\mathrm{tol}`. For most problems this is unreasonably small.
**lwork** : int
The dimension of the array :math:`\mathrm{work}`.
**Returns**
**a** : float, ndarray, shape :math:`\left(m, n\right)`
If :math:`\mathrm{svd}` is returned as :math:`\mathbf{False}`, :math:`\mathrm{a}` is overwritten by details of the :math:`QR` factorization of :math:`A`.
If :math:`\mathrm{svd}` is returned as :math:`\mathbf{True}`, the first :math:`n` rows of :math:`\mathrm{a}` are overwritten by the right-hand singular vectors, stored by rows; and the remaining rows of the array are used as workspace.
**b** : float, ndarray, shape :math:`\left(m\right)`
The first :math:`n` elements of :math:`\mathrm{b}` contain the minimal least squares solution vector :math:`x`. The remaining :math:`m-n` elements are used for workspace.
**svd** : bool
Is returned as :math:`\mathbf{False}` if the least squares solution has been obtained from the :math:`QR` factorization of :math:`A`. In this case :math:`A` is of full rank. :math:`\mathrm{svd}` is returned as :math:`\mathbf{True}` if the least squares solution has been obtained from the singular value decomposition of :math:`A`.
**sigma** : float
The standard error, i.e., the value :math:`\sqrt{r^\mathrm{T}r/\left(m-k\right)}` when :math:`m > k`, and the value zero when :math:`m = k`. Here :math:`r` is the residual vector :math:`b-Ax` and :math:`k` is the rank of :math:`A`.
**irank** : int
:math:`k`, the rank of the matrix :math:`A`. It should be noted that it is possible for :math:`\mathrm{irank}` to be returned as :math:`n` and :math:`\mathrm{svd}` to be returned as :math:`\mathbf{True}`. This means that the matrix :math:`R` only just failed the test for nonsingularity.
**work** : float, ndarray, shape :math:`\left(\mathrm{lwork}\right)`
If :math:`\mathrm{svd}` is returned as :math:`\mathbf{False}`, then the first :math:`n` elements of :math:`\mathrm{work}` contain information on the :math:`QR` factorization of :math:`A` (see argument :math:`\mathrm{a}` above), and :math:`\mathrm{work}[n]` contains the condition number :math:`\left\lVert R\right\rVert_E\left\lVert R^{-1}\right\rVert_E` of the upper triangular matrix :math:`R`.
If :math:`\mathrm{svd}` is returned as :math:`\mathbf{True}`, then the first :math:`n` elements of :math:`\mathrm{work}` contain the singular values of :math:`A` arranged in descending order and :math:`\mathrm{work}[n]` contains the total number of iterations taken by the :math:`QR` algorithm.
The rest of :math:`\mathrm{work}` is used as workspace.
.. _f04jg-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{lwork}` is too small. Minimum size required: :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 1`.
(`errno` :math:`1`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m\geq n`.
(`errno` :math:`2`)
The :math:`QR` algorithm has failed to converge to the singular values in :math:`50\times n` iterations. This failure can only happen when the singular value decomposition is employed, but even then it is not likely to occur.
.. _f04jg-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
The minimal least squares solution of the problem :math:`Ax = b` is the vector :math:`x` of minimum (Euclidean) length which minimizes the length of the residual vector :math:`r = b-Ax`.
The real :math:`m\times n\left(m\geq n\right)` matrix :math:`A` is factorized as
.. math::
A = Q\begin{pmatrix}R\\0\end{pmatrix}
where :math:`Q` is an :math:`m\times m` orthogonal matrix and :math:`R` is an :math:`n\times n` upper triangular matrix.
If :math:`R` is of full rank, then the least squares solution is given by
.. math::
x = \begin{pmatrix}R^{-1}&0\end{pmatrix}Q^\mathrm{T}b\text{.}
If :math:`R` is not of full rank, then the singular value decomposition of :math:`R` is obtained so that :math:`R` is factorized as
.. math::
R = UDV^\mathrm{T}\text{,}
where :math:`U` and :math:`V` are :math:`n\times n` orthogonal matrices and :math:`D` is the :math:`n\times n` diagonal matrix
.. math::
D = \mathrm{diag}\left(\sigma_1,\sigma_2,\ldots,\sigma_n\right)\text{,}
with :math:`\sigma_1\geq \sigma_2\geq \cdots \geq \sigma_n\geq 0`, these being the singular values of :math:`A`.
If the singular values :math:`\sigma_{{k+1}},\ldots,\sigma_n` are negligible, but :math:`\sigma_k` is not negligible, relative to the data errors in :math:`A`, then the rank of :math:`A` is taken to be :math:`k` and the minimal least squares solution is given by
.. math::
x = V\begin{pmatrix}S^{-1}&0\\0&0\end{pmatrix}\begin{pmatrix}U^\mathrm{T}&0\\0&I\end{pmatrix}Q^\mathrm{T}b\text{,}
where :math:`S = \mathrm{diag}\left(\sigma_1,\sigma_2,\ldots,\sigma_k\right)`.
The function also returns the value of the standard error
.. math::
\begin{array}{ccll}\sigma & = &\sqrt{\frac{{r^\mathrm{T}r}}{{m-k}}}\text{,}&\text{if }m > k\text{,}\\&&&\\& = &0\text{,}&\text{if }m = k\text{,}\end{array}
:math:`r^\mathrm{T}r` being the residual sum of squares and :math:`k` the rank of :math:`A`.
.. _f04jg-py2-py-references:
**References**
Lawson, C L and Hanson, R J, 1974, `Solving Least Squares Problems`, Prentice--Hall
"""
raise NotImplementedError
[docs]def real_tridiag_fac_solve(job, a, b, c, d, ipiv, y, tol):
r"""
``real_tridiag_fac_solve`` solves a system of tridiagonal equations following the factorization by :meth:`matop.real_gen_tridiag_lu <naginterfaces.library.matop.real_gen_tridiag_lu>`.
This function is intended for applications such as inverse iteration as well as straightforward linear equation applications.
.. _f04le-py2-py-doc:
For full information please refer to the NAG Library document for f04le
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04lef.html
.. _f04le-py2-py-parameters:
**Parameters**
**job** : int
Must specify the equations to be solved.
:math:`\mathrm{job} = 1`
The equations :math:`\left(T-\lambda I\right)x = y` are to be solved, but diagonal elements of :math:`U` are not to be perturbed.
:math:`\mathrm{job} = -1`
The equations :math:`\left(T-\lambda I\right)x = y` are to be solved and, if overflow would otherwise occur, diagonal elements of :math:`U` are to be perturbed. See argument :math:`\mathrm{tol}`.
:math:`\mathrm{job} = 2`
The equations :math:`\left(T-\lambda I\right)^\mathrm{T}x = y` are to be solved, but diagonal elements of :math:`U` are not to be perturbed.
:math:`\mathrm{job} = -2`
The equations :math:`\left(T-\lambda I\right)^\mathrm{T}x = y` are to be solved and, if overflow would otherwise occur, diagonal elements of :math:`U` are to be perturbed. See argument :math:`\mathrm{tol}`.
:math:`\mathrm{job} = 3`
The equations :math:`Ux = y` are to be solved, but diagonal elements of :math:`U` are not to be perturbed.
:math:`\mathrm{job} = -3`
The equations :math:`Ux = y` are to be solved and, if overflow would otherwise occur, diagonal elements of :math:`U` are to be perturbed. See argument :math:`\mathrm{tol}`.
**a** : float, array-like, shape :math:`\left(n\right)`
The diagonal elements of :math:`U` as returned by ``real_tridiag_fac_solve``.
**b** : float, array-like, shape :math:`\left(n\right)`
The elements of the first superdiagonal of :math:`U` as returned by ``real_tridiag_fac_solve``.
**c** : float, array-like, shape :math:`\left(n\right)`
The subdiagonal elements of :math:`L` as returned by ``real_tridiag_fac_solve``.
**d** : float, array-like, shape :math:`\left(n\right)`
The elements of the second superdiagonal of :math:`U` as returned by ``real_tridiag_fac_solve``.
**ipiv** : int, array-like, shape :math:`\left(n\right)`
Details of the matrix :math:`P` as returned by :meth:`matop.real_gen_tridiag_lu <naginterfaces.library.matop.real_gen_tridiag_lu>`.
**y** : float, array-like, shape :math:`\left(n\right)`
The right-hand side vector :math:`y`.
**tol** : float
The minimum perturbation to be made to very small diagonal elements of :math:`U`. :math:`\mathrm{tol}` is only referenced when :math:`\mathrm{job}` is negative. :math:`\mathrm{tol}` should normally be chosen as about :math:`\epsilon \left\lVert U\right\rVert`, where :math:`\epsilon` is the machine precision, but if :math:`\mathrm{tol}` is supplied as non-positive, it is reset to :math:`\epsilon \mathrm{max}\left(\left\lvert u_{{ij}}\right\rvert \right)`.
**Returns**
**y** : float, ndarray, shape :math:`\left(n\right)`
:math:`\mathrm{y}` is overwritten by the solution vector :math:`x`.
**tol** : float
If on entry :math:`\mathrm{tol}` is non-positive, it is reset as just described. Otherwise :math:`\mathrm{tol}` is unchanged.
.. _f04le-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`i > 1`)
Overflow would occur computing the element :math:`I` in the solution. :math:`I = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{job} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{job}\neq 0` and :math:`\mathrm{job}\geq -3` and :math:`\mathrm{job}\leq 3`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 1`.
.. _f04le-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
Following the factorization of the :math:`n\times n` tridiagonal matrix :math:`\left(T-\lambda I\right)` as
.. math::
T-\lambda I = PLU
by :meth:`matop.real_gen_tridiag_lu <naginterfaces.library.matop.real_gen_tridiag_lu>`, ``real_tridiag_fac_solve`` may be used to solve any of the equations
.. math::
\left(T-\lambda I\right)x = y\text{, }\quad \left(T-\lambda I\right)^\mathrm{T}x = y\text{, }\quad Ux = y
for :math:`x`, the choice of equation being controlled by the argument :math:`\mathrm{job}`.
In each case there is an option to perturb zero or very small diagonal elements of :math:`U`, this option being intended for use in applications such as inverse iteration.
.. _f04le-py2-py-references:
**References**
Wilkinson, J H, 1965, `The Algebraic Eigenvalue Problem`, Oxford University Press, Oxford
Wilkinson, J H and Reinsch, C, 1971, `Handbook for Automatic Computation II, Linear Algebra`, Springer--Verlag
"""
raise NotImplementedError
[docs]def real_blkdiag_fac_solve(trans, blkstr, a, pivot, b):
r"""
``real_blkdiag_fac_solve`` calculates the approximate solution of a set of real linear equations with multiple right-hand sides, :math:`AX = B` or :math:`A^\mathrm{T}X = B`, where :math:`A` is an almost block-diagonal matrix which has been factorized by :meth:`matop.real_gen_blkdiag_lu <naginterfaces.library.matop.real_gen_blkdiag_lu>`.
.. _f04lh-py2-py-doc:
For full information please refer to the NAG Library document for f04lh
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04lhf.html
.. _f04lh-py2-py-parameters:
**Parameters**
**trans** : str, length 1
Specifies the equations to be solved.
:math:`\mathrm{trans} = \texttt{'N'}`
Solve :math:`AX = B`.
:math:`\mathrm{trans} = \texttt{'T'}`
Solve :math:`A^\mathrm{T}X = B`.
**blkstr** : int, array-like, shape :math:`\left(3, \textit{nbloks}\right)`
Information which describes the block structure of :math:`A`, as supplied to ``real_blkdiag_fac_solve``.
**a** : float, array-like, shape :math:`\left(\textit{lena}\right)`
The elements in the factorization of :math:`A`, as returned by ``real_blkdiag_fac_solve``.
**pivot** : int, array-like, shape :math:`\left(n\right)`
Details of the interchanges in the factorization, as returned by ``real_blkdiag_fac_solve``.
**b** : float, array-like, shape :math:`\left(n, \textit{ir}\right)`
The :math:`n\times r` right-hand side matrix :math:`B`.
**Returns**
**b** : float, ndarray, shape :math:`\left(n, \textit{ir}\right)`
:math:`\mathrm{b}` is overwritten by the :math:`n\times r` solution matrix :math:`X`.
.. _f04lh-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{lena}` is too small. :math:`\textit{lena} = \langle\mathit{\boldsymbol{value}}\rangle`. Minimum possible dimension: :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, the following equality does not hold: :math:`\mathrm{blkstr}[1,0]+\mathrm{sum}\left(\mathrm{blkstr}[1,k-1]-\mathrm{blkstr}[2,k-2]:k = 2,\textit{nbloks}\right) = n`.
(`errno` :math:`1`)
On entry, the following equality does not hold: :math:`\mathrm{sum}\left({\mathrm{blkstr}[0,k-1]:k = 1}, \textit{nbloks}\right) = n`.
(`errno` :math:`1`)
On entry, :math:`K = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{blkstr}[2,K-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{blkstr}[2,K-2] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{blkstr}[1,K-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{blkstr}[2,K-1]+\mathrm{blkstr}[2,K-2]\geq \mathrm{blkstr}[1,K-1]`.
(`errno` :math:`1`)
On entry, the following inequality was not satisfied for: :math:`J = \langle\mathit{\boldsymbol{value}}\rangle`. :math:`\mathrm{sum}\left({\mathrm{blkstr}[1,k-1]-\mathrm{blkstr}[2,k-1]:k = 1}, J\right)\leq` :math:`\mathrm{sum}\left({\mathrm{blkstr}[0,k-1]:k = 1}, J\right)\leq` :math:`\mathrm{blkstr}[1,0]+\mathrm{sum}\left({\mathrm{blkstr}[1,k-1]-\mathrm{blkstr}[2,k-2]:k = 2}, J\right)`.
(`errno` :math:`1`)
On entry, :math:`K = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{blkstr}[1,K-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{blkstr}[2,K-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{blkstr}[0,K-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{blkstr}[1,K-1]-\mathrm{blkstr}[2,K-1]\leq \mathrm{blkstr}[0,K-1]`.
(`errno` :math:`1`)
On entry, :math:`K = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{blkstr}[1,K-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{blkstr}[0,K-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{blkstr}[1,K-1]\geq \mathrm{blkstr}[0,K-1]`.
(`errno` :math:`1`)
On entry, :math:`K = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{blkstr}[2,K-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{blkstr}[2,K-1]\geq 0`.
(`errno` :math:`1`)
On entry, :math:`K = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{blkstr}[1,K-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{blkstr}[1,K-1]\geq 1`.
(`errno` :math:`1`)
On entry, :math:`K = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{blkstr}[0,K-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{blkstr}[0,K-1]\geq 1`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{nbloks} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq \textit{nbloks}`.
(`errno` :math:`1`)
On entry, :math:`\textit{nbloks} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nbloks}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\textit{ir} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ir}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{trans} \neq \texttt{'N'}` or :math:`\texttt{'T'}`: :math:`\mathrm{trans} = \langle\mathit{\boldsymbol{value}}\rangle`.
.. _f04lh-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``real_blkdiag_fac_solve`` solves a set of real linear equations :math:`AX = B` or :math:`A^\mathrm{T}X = B`, where :math:`A` is almost block-diagonal. :math:`A` must first be factorized by :meth:`matop.real_gen_blkdiag_lu <naginterfaces.library.matop.real_gen_blkdiag_lu>`. ``real_blkdiag_fac_solve`` then computes :math:`X` by forward and backward substitution over the blocks.
.. _f04lh-py2-py-references:
**References**
Diaz, J C, Fairweather, G and Keast, P, 1983, `Fortran packages for solving certain almost block diagonal linear systems by modified alternate row and column elimination`, ACM Trans. Math. Software (9), 358--375
"""
raise NotImplementedError
[docs]def real_posdef_vband_solve(al, d, nrow, b, iselct):
r"""
``real_posdef_vband_solve`` computes the approximate solution of a system of real linear equations with multiple right-hand sides, :math:`AX = B`, where :math:`A` is a symmetric positive definite variable-bandwidth matrix, which has previously been factorized by :meth:`matop.real_vband_posdef_fac <naginterfaces.library.matop.real_vband_posdef_fac>`.
Related systems may also be solved.
.. _f04mc-py2-py-doc:
For full information please refer to the NAG Library document for f04mc
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04mcf.html
.. _f04mc-py2-py-parameters:
**Parameters**
**al** : float, array-like, shape :math:`\left(\textit{lal}\right)`
The elements within the envelope of the lower triangular matrix :math:`L`, taken in row by row order, as returned by :meth:`matop.real_vband_posdef_fac <naginterfaces.library.matop.real_vband_posdef_fac>`. The unit diagonal elements of :math:`L` must be stored explicitly.
**d** : float, array-like, shape :math:`\left(:\right)`
Note: the required length for this argument is determined as follows: if :math:`\mathrm{iselct}\geq 4`: :math:`1`; otherwise: :math:`n`.
The diagonal elements of the diagonal matrix :math:`D`. :math:`\mathrm{d}` is not referenced if :math:`\mathrm{iselct}\geq 4`.
**nrow** : int, array-like, shape :math:`\left(n\right)`
:math:`\mathrm{nrow}[i-1]` must contain the width of row :math:`i` of :math:`L`, i.e., the number of elements between the first (leftmost) nonzero element and the element on the diagonal, inclusive.
**b** : float, array-like, shape :math:`\left(n, \textit{ir}\right)`
The :math:`n\times r` right-hand side matrix :math:`B`. See also `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04mcf.html#fcomments>`__.
**iselct** : int
Must specify the type of system to be solved, as follows:
:math:`\mathrm{iselct} = 1`
Solve :math:`LDL^\mathrm{T}X = B`.
:math:`\mathrm{iselct} = 2`
Solve :math:`LDX = B`.
:math:`\mathrm{iselct} = 3`
Solve :math:`DL^\mathrm{T}X = B`.
:math:`\mathrm{iselct} = 4`
Solve :math:`LL^\mathrm{T}X = B`.
:math:`\mathrm{iselct} = 5`
Solve :math:`LX = B`.
:math:`\mathrm{iselct} = 6`
Solve :math:`L^\mathrm{T}X = B`.
**Returns**
**x** : float, ndarray, shape :math:`\left(n, \textit{ir}\right)`
The :math:`n\times r` solution matrix :math:`X`. See also `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04mcf.html#fcomments>`__.
.. _f04mc-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{lal} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{nrow}[0] + \cdots +\mathrm{nrow}[n-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{lal}\geq \mathrm{nrow}[0] + \cdots +\mathrm{nrow}[n-1]`.
(`errno` :math:`1`)
On entry, :math:`I = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{nrow}[I-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{nrow}[I-1]\geq 1` and :math:`\mathrm{nrow}[I-1]\leq I`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 1`.
(`errno` :math:`2`)
On entry, :math:`\textit{ir} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ir}\geq 1`.
(`errno` :math:`3`)
On entry, :math:`\mathrm{iselct} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{iselct}\geq 1` and :math:`\mathrm{iselct}\leq 6`.
(`errno` :math:`4`)
On entry, :math:`I = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{d}[I-1]\neq 0.0`.
(`errno` :math:`5`)
At least one diagonal entry of :math:`\mathrm{al}` is not unit.
.. _f04mc-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.`
The normal use of this function is the solution of the systems :math:`AX = B`, following a call of :meth:`matop.real_vband_posdef_fac <naginterfaces.library.matop.real_vband_posdef_fac>` to determine the Cholesky factorization :math:`A = LDL^\mathrm{T}` of the symmetric positive definite variable-bandwidth matrix :math:`A`.
However, the function may be used to solve any one of the following systems of linear algebraic equations:
(1) :math:`LDL^\mathrm{T}X = B` (usual system),
(#) :math:`LDX = B` (lower triangular system),
(#) :math:`DL^\mathrm{T}X = B` (upper triangular system),
(#) :math:`LL^\mathrm{T}X = B`
(#) :math:`LX = B` (unit lower triangular system),
(#) :math:`L^\mathrm{T}X = B` (unit upper triangular system).
:math:`L` denotes a unit lower triangular variable-bandwidth matrix of order :math:`n`, :math:`D` a diagonal matrix of order :math:`n`, and :math:`B` a set of right-hand sides.
The matrix :math:`L` is represented by the elements lying within its **envelope**, i.e., between the first nonzero of each row and the diagonal.
The width :math:`\mathrm{nrow}[i-1]` of the :math:`i`\ th row is the number of elements between the first nonzero element and the element on the diagonal inclusive.
.. _f04mc-py2-py-references:
**References**
Wilkinson, J H and Reinsch, C, 1971, `Handbook for Automatic Computation II, Linear Algebra`, Springer--Verlag
"""
raise NotImplementedError
[docs]def real_toeplitz_yule_update(t, x, v):
r"""
``real_toeplitz_yule_update`` updates the solution to the Yule--Walker equations for a real symmetric positive definite Toeplitz system.
.. _f04me-py2-py-doc:
For full information please refer to the NAG Library document for f04me
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04mef.html
.. _f04me-py2-py-parameters:
**Parameters**
**t** : float, array-like, shape :math:`\left(n+1\right)`
:math:`\mathrm{t}[0]` must contain the value :math:`\tau_0` of the diagonal elements of :math:`T`, and the remaining :math:`\textit{n}` elements of :math:`\mathrm{t}` must contain the elements of the vector :math:`t_n`.
**x** : float, array-like, shape :math:`\left(n\right)`
With :math:`n > 1` the (:math:`n-1`) elements of the solution vector :math:`x_{{n-1}}` as returned by a previous call to ``real_toeplitz_yule_update``. The element :math:`\mathrm{x}[n-1]` need not be specified.
**v** : float
With :math:`n > 1` the mean square prediction error for the (:math:`n-1`)th step, as returned by a previous call to ``real_toeplitz_yule_update``.
**Returns**
**x** : float, ndarray, shape :math:`\left(n\right)`
The solution vector :math:`x_n`. The element :math:`\mathrm{x}[n-1]` returns the partial (auto)correlation coefficient, or reflection coefficient, for the :math:`n`\ th step. If :math:`\left\lvert \mathrm{x}[n-1]\right\rvert \geq 1.0`, the matrix :math:`T_{{n+1}}` will not be positive definite to working accuracy.
**v** : float
The mean square prediction error, or predictor error variance ratio, :math:`\nu_n`, for the :math:`n`\ th step. (See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04mef.html#fcomments>`__ and `the G13 Introduction <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/g13/g13intro.html>`__.)
.. _f04me-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{x}[n-2] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`n > 1`, :math:`\left\lvert \mathrm{x}[n-2]\right\rvert < 1.0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{t}[0] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{t}[0] > 0.0`.
(`errno` :math:`-1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`1`)
Matrix of order :math:`\langle\mathit{\boldsymbol{value}}\rangle` would not be positive definite. Value of the reflection coefficient is :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
.. _f04me-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``real_toeplitz_yule_update`` solves the equations
.. math::
T_nx_n = -t_n\text{,}
where :math:`T_n` is the :math:`n\times n` symmetric positive definite Toeplitz matrix
.. math::
T_n = \begin{pmatrix}\tau_0&\tau_1&\tau_2&\ldots &\tau_{{n-1}}\\\tau_1&\tau_0&\tau_1&\ldots &\tau_{{n-2}}\\\tau_2&\tau_1&\tau_0&\ldots &\tau_{{n-3}}\\.&.&.&&.\\\tau_{{n-1}}&\tau_{{n-2}}&\tau_{{n-3}}&\ldots &\tau_0\end{pmatrix}
and :math:`t_n` is the vector
.. math::
t_n^\mathrm{T} = \left(\tau_1\tau_2\ldots \tau_n\right)\text{,}
given the solution of the equations
.. math::
T_{{n-1}}x_{{n-1}} = -t_{{n-1}}\text{.}
The function will normally be used to successively solve the equations
.. math::
T_kx_k = -t_k\text{, }\quad k = 1,2,\ldots,n\text{.}
If it is desired to solve the equations for a single value of :math:`n`, then function :meth:`real_toeplitz_yule` may be called.
This function uses the method of Durbin (see Durbin (1960) and Golub and Van Loan (1996)).
.. _f04me-py2-py-references:
**References**
Bunch, J R, 1985, `Stability of methods for solving Toeplitz systems of equations`, SIAM J. Sci. Statist. Comput. (6), 349--364
Bunch, J R, 1987, `The weak and strong stability of algorithms in numerical linear algebra`, Linear Algebra Appl. (88/89), 49--66
Cybenko, G, 1980, `The numerical stability of the Levinson--Durbin algorithm for Toeplitz systems of equations`, SIAM J. Sci. Statist. Comput. (1), 303--319
Durbin, J, 1960, `The fitting of time series models`, Rev. Inst. Internat. Stat. (28), 233
Golub, G H and Van Loan, C F, 1996, `Matrix Computations`, (3rd Edition), Johns Hopkins University Press, Baltimore
"""
raise NotImplementedError
[docs]def real_toeplitz_update(t, b, x, work):
r"""
``real_toeplitz_update`` updates the solution of the equations :math:`Tx = b`, where :math:`T` is a real symmetric positive definite Toeplitz matrix.
.. _f04mf-py2-py-doc:
For full information please refer to the NAG Library document for f04mf
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04mff.html
.. _f04mf-py2-py-parameters:
**Parameters**
**t** : float, array-like, shape :math:`\left(n\right)`
:math:`\mathrm{t}[\textit{i}]` must contain the value :math:`\tau_{\textit{i}}`, for :math:`\textit{i} = 0,1,\ldots,n-1`.
**b** : float, array-like, shape :math:`\left(n\right)`
The right-hand side vector :math:`b_n`.
**x** : float, array-like, shape :math:`\left(n\right)`
With :math:`n > 1` the (:math:`n-1`) elements of the solution vector :math:`x_{{n-1}}` as returned by a previous call to ``real_toeplitz_update``. The element :math:`\mathrm{x}[n-1]` need not be specified.
**work** : float, array-like, shape :math:`\left(2\times n-1\right)`
With :math:`n > 2` the elements of :math:`\mathrm{work}` should be as returned from a previous call to ``real_toeplitz_update`` with (:math:`n-1`) as the argument :math:`\textit{n}`.
**Returns**
**x** : float, ndarray, shape :math:`\left(n\right)`
The solution vector :math:`x_n`.
**p** : float
The reflection coefficient :math:`p_{{n-1}}`. (See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04mff.html#fcomments>`__.)
**work** : float, ndarray, shape :math:`\left(2\times n-1\right)`
The first (:math:`n-1`) elements of :math:`\mathrm{work}` contain the solution to the Yule--Walker equations
.. math::
T_{{n-1}}y_{{n-1}} = -t_{{n-1}}\text{,}
where :math:`t_{{n-1}} = {\left(\tau_1\tau_2\ldots \tau_{{n-1}}\right)}^\mathrm{t}`.
.. _f04mf-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-1`)
On entry, :math:`\mathrm{t}[0] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{t}[0] > 0.0`.
(`errno` :math:`-1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`1`)
The Toeplitz Matrix is not positive definite Value of the reflection coefficient is :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
.. _f04mf-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``real_toeplitz_update`` solves the equations
.. math::
T_nx_n = b_n\text{,}
where :math:`T_n` is the :math:`n\times n` symmetric positive definite Toeplitz matrix
.. math::
T_n = \begin{pmatrix}\tau_0&\tau_1&\tau_2&\ldots &\tau_{{n-1}}\\\tau_1&\tau_0&\tau_1&\ldots &\tau_{{n-2}}\\\tau_2&\tau_1&\tau_0&\ldots &\tau_{{n-3}}\\.&.&.&&.\\\tau_{{n-1}}&\tau_{{n-2}}&\tau_{{n-3}}&\ldots &\tau_0\end{pmatrix}
and :math:`b_n` is the :math:`n`-element vector :math:`b_n = \left(\beta_1\beta_2\ldots \beta_n\right)^\mathrm{T}`, given the solution of the equations
.. math::
T_{{n-1}}x_{{n-1}} = b_{{n-1}}\text{.}
This function will normally be used to successively solve the equations
.. math::
T_kx_k = b_k\text{, }\quad k = 1,2,\ldots,n\text{.}
If it is desired to solve the equations for a single value of :math:`n`, then function :meth:`real_toeplitz_solve` may be called.
This function uses the method of Levinson (see Levinson (1947) and Golub and Van Loan (1996)).
.. _f04mf-py2-py-references:
**References**
Bunch, J R, 1985, `Stability of methods for solving Toeplitz systems of equations`, SIAM J. Sci. Statist. Comput. (6), 349--364
Bunch, J R, 1987, `The weak and strong stability of algorithms in numerical linear algebra`, Linear Algebra Appl. (88/89), 49--66
Cybenko, G, 1980, `The numerical stability of the Levinson--Durbin algorithm for Toeplitz systems of equations`, SIAM J. Sci. Statist. Comput. (1), 303--319
Golub, G H and Van Loan, C F, 1996, `Matrix Computations`, (3rd Edition), Johns Hopkins University Press, Baltimore
Levinson, N, 1947, `The Weiner RMS error criterion in filter design and prediction`, J. Math. Phys. (25), 261--278
"""
raise NotImplementedError
[docs]def real_gen_sparse_lsqsol(n, b, aprod, damp=0.0, atol=-1.0, btol=-1.0, conlim=0.0, itnlim=0, msglvl=0, data=None, io_manager=None):
r"""
``real_gen_sparse_lsqsol`` solves sparse nonsymmetric equations, sparse linear least squares problems and sparse damped linear least squares problems, using a Lanczos algorithm.
.. _f04qa-py2-py-doc:
For full information please refer to the NAG Library document for f04qa
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html
.. _f04qa-py2-py-parameters:
**Parameters**
**n** : int
:math:`n`, the number of columns of the matrix :math:`A`.
**b** : float, array-like, shape :math:`\left(m\right)`
The right-hand side vector :math:`b`.
**aprod** : callable (mode, x, y) = aprod(mode, x, y, data=None)
:math:`\mathrm{aprod}` must perform the operations :math:`y \mathrel{:=} y+Ax` and :math:`x \mathrel{:=} x+A^\mathrm{T}y` for given vectors :math:`x` and :math:`y`.
**Parameters**
**mode** : int
Specifies which operation is to be performed.
:math:`\mathrm{mode} = 1`
:math:`\mathrm{aprod}` must compute :math:`y+Ax`.
:math:`\mathrm{mode} = 2`
:math:`\mathrm{aprod}` must compute :math:`x+A^\mathrm{T}y`.
**x** : float, ndarray, shape :math:`\left(n\right)`
The vector :math:`x`.
**y** : float, ndarray, shape :math:`\left(m\right)`
The vector :math:`y`.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**mode** : int
May be used as a flag to indicate a failure in the computation of :math:`y+Ax` or :math:`x+A^\mathrm{T}y`. If :math:`\mathrm{mode}` is negative on exit from :math:`\mathrm{aprod}`, ``real_gen_sparse_lsqsol`` will exit immediately with :math:`\textit{errno}` set to :math:`\mathrm{mode}`.
**x** : float, array-like, shape :math:`\left(n\right)`
If :math:`\mathrm{mode} = 1`, :math:`\mathrm{x}` must be unchanged.
If :math:`\mathrm{mode} = 2`, :math:`\mathrm{x}` must contain :math:`x+A^\mathrm{T}y`.
**y** : float, array-like, shape :math:`\left(m\right)`
If :math:`\mathrm{mode} = 1`, :math:`\mathrm{y}` must contain :math:`y+Ax`.
If :math:`\mathrm{mode} = 2`, :math:`\mathrm{y}` must be unchanged.
**damp** : float, optional
The value :math:`\lambda`. If either problem `(1) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn1>`__ or problem `(2) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn2>`__ is to be solved, :math:`\mathrm{damp}` must be supplied as zero.
**atol** : float, optional
The tolerance, :math:`\textit{tol}_1`, of the convergence criteria `(6) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn6>`__ and `(7) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn6>`__; it should be an estimate of the largest relative error in the elements of :math:`A`. For example, if the elements of :math:`A` are correct to about :math:`4` significant figures, then :math:`\mathrm{atol}` should be set to about :math:`5\times 10^{-4}`. If :math:`\mathrm{atol}` is supplied as less than :math:`\epsilon`, where :math:`\epsilon` is the machine precision, the value :math:`\epsilon` is used instead of :math:`\mathrm{atol}`.
**btol** : float, optional
The tolerance, :math:`\textit{tol}_2`, of the convergence criterion `(6) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn6>`__; it should be an estimate of the largest relative error in the elements of :math:`B`. For example, if the elements of :math:`B` are correct to about :math:`4` significant figures, then :math:`\mathrm{btol}` should be set to about :math:`5\times 10^{-4}`. If :math:`\mathrm{btol}` is supplied as less than :math:`\epsilon`, the value :math:`\epsilon` is used instead of :math:`\mathrm{btol}`.
**conlim** : float, optional
The value :math:`c_{\mathrm{lim}}` of equation `(8) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn8>`__; it should be an upper limit on the condition number of :math:`\bar{A}`. :math:`\mathrm{conlim}` should not normally be chosen much larger than :math:`1.0/\mathrm{atol}`. If :math:`\mathrm{conlim}` is supplied as zero, the value :math:`1.0/\epsilon` is used instead of :math:`\mathrm{conlim}`.
**itnlim** : int, optional
An upper limit on the number of iterations. If :math:`\mathrm{itnlim}\leq 0`, the value :math:`\mathrm{n}` is used in place of :math:`\mathrm{itnlim}`, but for ill-conditioned problems a higher value of :math:`\mathrm{itnlim}` is likely to be necessary.
**msglvl** : int, optional
The level of printing from ``real_gen_sparse_lsqsol``. If :math:`\mathrm{msglvl}\leq 0`, then no printing occurs, but otherwise messages will be output on the file object associated with the advisory I/O unit (see :class:`~naginterfaces.base.utils.FileObjManager`). A description of the printed output is given in `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#fc-printedoutput>`__. The level of printing is determined as follows:
:math:`\mathrm{msglvl}\leq 0`
No printing.
:math:`\mathrm{msglvl} = 1`
A brief summary is printed just prior to return from ``real_gen_sparse_lsqsol``.
:math:`\mathrm{msglvl}\geq 2`
A summary line is printed periodically to monitor the progress of ``real_gen_sparse_lsqsol``, together with a brief summary just prior to return from ``real_gen_sparse_lsqsol``.
**data** : arbitrary, optional
User-communication data for callback functions.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**b** : float, ndarray, shape :math:`\left(m\right)`
:math:`\mathrm{b}` is overwritten.
**x** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
The solution vector :math:`x`.
**se** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
The estimates of the standard errors of the components of :math:`x`. Thus :math:`\mathrm{se}[i-1]` contains an estimate of :math:`\sqrt{\nu_{{ii}}}`, where :math:`\nu_{{ii}}` is the :math:`i`\ th diagonal element of the estimated variance-covariance matrix :math:`V`. The estimates returned in :math:`\mathrm{se}` will be the lower bounds on the actual estimated standard errors, but will usually be correct to at least one significant figure.
**itnlim** : int
Unchanged unless :math:`\mathrm{itnlim}\leq 0` on entry, in which case it is set to :math:`\mathrm{n}`.
**itn** : int
The number of iterations performed.
**anorm** : float
An estimate of :math:`\left\lVert \bar{A}\right\rVert` for the matrix :math:`\bar{A}` of `(4) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn4>`__.
**acond** : float
An estimate of :math:`\mathrm{cond}\left(\bar{A}\right)` which is a lower bound.
**rnorm** : float
An estimate of :math:`\left\lVert \bar{r}\right\rVert` for the residual, :math:`\bar{r}`, of `(5) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn5>`__ corresponding to the solution :math:`x` returned in :math:`\mathrm{x}`. Note that :math:`\left\lVert \bar{r}\right\rVert` is the function being minimized.
**arnorm** : float
An estimate of the :math:`\left\lVert \bar{A}^\mathrm{T}\bar{r}\right\rVert` corresponding to the solution :math:`x` returned in :math:`\mathrm{x}`.
**xnorm** : float
An estimate of :math:`\left\lVert x\right\rVert` for the solution :math:`x` returned in :math:`\mathrm{x}`.
**inform** : int
The reason for termination of ``real_gen_sparse_lsqsol``.
:math:`\mathrm{inform} = 0`
The exact solution is :math:`x = 0`. No iterations are performed in this case.
:math:`\mathrm{inform} = 1`
The termination criterion of `(6) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn6>`__ has been satisfied with :math:`\textit{tol}_1` and :math:`\textit{tol}_2` as the values supplied in :math:`\mathrm{atol}` and :math:`\mathrm{btol}` respectively.
:math:`\mathrm{inform} = 2`
The termination criterion of `(7) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn7>`__ has been satisfied with :math:`\textit{tol}_1` as the value supplied in :math:`\mathrm{atol}`.
:math:`\mathrm{inform} = 3`
The termination criterion of `(6) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn6>`__ has been satisfied with :math:`\textit{tol}_1` and/or :math:`\textit{tol}_2` as the value :math:`\epsilon`, where :math:`\epsilon` is the machine precision. One or both of the values supplied in :math:`\mathrm{atol}` and :math:`\mathrm{btol}` must have been less than :math:`\epsilon` and was too small for this machine.
:math:`\mathrm{inform} = 4`
The termination criterion of `(7) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn7>`__ has been satisfied with :math:`\textit{tol}_1` as the value :math:`\epsilon`. The value supplied in :math:`\mathrm{atol}` must have been less than :math:`\epsilon` and was too small for this machine.
The values :math:`\mathrm{inform} = 5`, :math:`6` or :math:`7` correspond to failure with :math:`\mathrm{errno}` = 2, 3 and 4 respectively (see :ref:`Exceptions <f04qa-py2-py-errors>`) and when :math:`\textit{errno}` is negative :math:`\mathrm{inform}` will be set to the same negative value.
.. _f04qa-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m\geq 1`.
(`errno` :math:`2`)
Termination criteria not satisfied, but the condition number bound supplied in :math:`\mathrm{conlim}` has been met.
(`errno` :math:`3`)
Termination criteria not satisfied, but the condition number is bounded by :math:`1`/:meth:`machine.precision <naginterfaces.library.machine.precision>`.
(`errno` :math:`4`)
Maximum number of iterations reached.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`i < 0`)
:math:`\mathrm{mode}` in :math:`\mathrm{aprod}` has been set to: :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
.. _f04qa-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``real_gen_sparse_lsqsol`` can be used to solve a system of linear equations
.. math::
Ax = b
where :math:`A` is an :math:`n\times n` sparse nonsymmetric matrix, or can be used to solve linear least squares problems, so that ``real_gen_sparse_lsqsol`` minimizes the value :math:`\rho` given by
.. math::
\rho = \left\lVert r\right\rVert \text{, }\quad r = b-Ax
where :math:`A` is an :math:`m\times n` sparse matrix and :math:`\left\lVert r\right\rVert` denotes the Euclidean length of :math:`r` so that :math:`\left\lVert r\right\rVert^2 = r^\mathrm{T}r`.
A damping argument, :math:`\lambda`, may be included in the least squares problem in which case ``real_gen_sparse_lsqsol`` minimizes the value :math:`\rho` given by
.. math::
\rho^2 = \left\lVert r\right\rVert^2+\lambda^2\left\lVert x\right\rVert^2\text{.}
:math:`\lambda` is supplied as the argument :math:`\mathrm{damp}` and should of course be zero if the solution to problems `(1) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn1>`__ and `(2) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn1>`__ is required.
Minimizing :math:`\rho` in `(3) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn3>`__ is often called ridge regression.
``real_gen_sparse_lsqsol`` is based upon algorithm LSQR (see Paige and Saunders (1982a) and Paige and Saunders (1982b)) and solves the problems by an algorithm based upon the Lanczos process.
The function does not require :math:`A` explicitly, but :math:`A` is specified via :math:`\mathrm{aprod}` which must perform the operations :math:`\left(y+Ax\right)` and :math:`\left(x+A^\mathrm{T}y\right)` for a given :math:`n`-element vector :math:`x` and :math:`m` element vector :math:`y`. An argument to :math:`\mathrm{aprod}` specifies which of the two operations is required on a given entry.
The function also returns estimates of the standard errors of the sample regression coefficients ( :math:`x_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,n`) given by the diagonal elements of the estimated variance-covariance matrix :math:`V`.
When problem `(2) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn2>`__ is being solved and :math:`A` is of full rank, then :math:`V` is given by
.. math::
V = s^2\left(A^\mathrm{T}A\right)^{-1}\text{, }\quad s^2 = \rho^2/\left(m-n\right)\text{, }\quad m > n
and when problem `(3) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn3>`__ is being solved then :math:`V` is given by
.. math::
V = s^2{\left(A^\mathrm{T}A+\lambda^2I\right)}^{-1}\text{, }\quad s^2 = \rho^2/m\text{, }\quad \lambda \neq 0\text{.}
Let :math:`\bar{A}` denote the matrix
.. math::
\bar{A} = A\text{, }\quad \lambda = 0\text{; }\quad \bar{A} = \begin{pmatrix}A\\\lambda I\end{pmatrix}\text{, }\quad \lambda \neq 0\text{,}
let :math:`\bar{r}` denote the residual vector
.. math::
\bar{r} = r\text{, }\quad \lambda = 0\text{; }\quad \bar{r} = \begin{pmatrix}b\\0\end{pmatrix}-\bar{A}x\text{, }\quad \lambda \neq 0
corresponding to an iterate :math:`x`, so that :math:`\rho = \left\lVert \bar{r}\right\rVert` is the function being minimized, and let :math:`\left\lVert A\right\rVert` denote the Frobenius (Euclidean) norm of :math:`A`.
Then the function accepts :math:`x` as a solution if it is estimated that one of the following two conditions is satisfied:
.. math::
\rho \leq \textit{tol}_1\left\lVert \bar{A}\right\rVert.\left\lVert x\right\rVert +\textit{tol}_2\left\lVert b\right\rVert
.. math::
\left\lVert \bar{A}^\mathrm{T}\bar{r}\right\rVert \leq \textit{tol}_1\left\lVert \bar{A}\right\rVert \rho
where :math:`\textit{tol}_1` and :math:`\textit{tol}_2` are user-supplied tolerances which estimate the relative errors in :math:`A` and :math:`b` respectively.
Condition `(6) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn6>`__ is appropriate for compatible problems where, in theory, we expect the residual to be zero and will be satisfied by an acceptable solution :math:`x` to a compatible problem.
Condition `(7) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn7>`__ is appropriate for incompatible systems where we do not expect the residual to be zero and is based on the observation that, in theory,
.. math::
\bar{A}^\mathrm{T}\bar{r} = 0
when :math:`x` is a solution to the least squares problem, and so `(7) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn7>`__ will be satisfied by an acceptable solution :math:`x` to a linear least squares problem.
The function also includes a test to prevent convergence to solutions, :math:`x`, with unacceptably large elements.
This can happen if :math:`A` is nearly singular or is nearly rank deficient.
If we let the singular values of :math:`\bar{A}` be
.. math::
\sigma_1\geq \sigma_2\geq \cdots \geq \sigma_n\geq 0
then the condition number of :math:`\bar{A}` is defined as
.. math::
\mathrm{cond}\left(\bar{A}\right) = \sigma_1/\sigma_k
where :math:`\sigma_k` is the smallest nonzero singular value of :math:`\bar{A}` and hence :math:`k` is the rank of :math:`\bar{A}`.
When :math:`k < n`, then :math:`\bar{A}` is rank deficient, the least squares solution is not unique and ``real_gen_sparse_lsqsol`` will normally converge to the minimal length solution.
In practice :math:`\bar{A}` will not have exactly zero singular values, but may instead have small singular values that we wish to regard as zero.
The function provides for this possibility by terminating if
.. math::
\mathrm{cond}\left(\bar{A}\right)\geq c_{\mathrm{lim}}
where :math:`c_{\mathrm{lim}}` is a user-supplied limit on the condition number of :math:`\bar{A}`.
For problem `(1) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn1>`__ termination with this condition indicates that :math:`A` is nearly singular and for problem `(2) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn2>`__ indicates that :math:`A` is nearly rank deficient and so has near linear dependencies in its columns.
In this case inspection of :math:`\left\lVert r\right\rVert`, :math:`\left\lVert A^\mathrm{T}r\right\rVert` and :math:`\left\lVert x\right\rVert`, which are all returned by the function, will indicate whether or not an acceptable solution has been found.
Condition `(8) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04qaf.html#eqn8>`__, perhaps in conjunction with :math:`\lambda \neq 0`, can be used to try and 'regularize' least squares solutions.
A full discussion of the stopping criteria is given in Section 6 of Paige and Saunders (1982a).
Introduction of a nonzero damping argument :math:`\lambda` tends to reduce the size of the computed solution and to make its components less sensitive to changes in the data, and ``real_gen_sparse_lsqsol`` is applicable when a value of :math:`\lambda` is known `a priori`.
To have an effect, :math:`\lambda` should normally be at least :math:`\sqrt{\epsilon }\left\lVert A\right\rVert` where :math:`\epsilon` is the machine precision.
For further discussion see Paige and Saunders (1982b) and the references given there.
Whenever possible the matrix :math:`A` should be scaled so that the relative errors in the elements of :math:`A` are all of comparable size.
Such a scaling helps to prevent the least squares problem from being unnecessarily sensitive to data errors and will normally reduce the number of iterations required.
At the very least, in the absence of better information, the columns of :math:`A` should be scaled to have roughly equal column length.
.. _f04qa-py2-py-references:
**References**
Paige, C C and Saunders, M A, 1982, `LSQR: An algorithm for sparse linear equations and sparse least squares`, ACM Trans. Math. Software (8), 43--71
Paige, C C and Saunders, M A, 1982, `Algorithm 583 LSQR: Sparse linear equations and least squares problems`, ACM Trans. Math. Software (8), 195--209
"""
raise NotImplementedError
[docs]def real_gen_lsq_covmat(job, sigma, a, svd, irank, sv):
r"""
``real_gen_lsq_covmat`` returns elements of the estimated variance-covariance matrix of the sample regression coefficients for the solution of a linear least squares problem.
The function can be used to find the estimated variances of the sample regression coefficients.
.. _f04ya-py2-py-doc:
For full information please refer to the NAG Library document for f04ya
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04yaf.html
.. _f04ya-py2-py-parameters:
**Parameters**
**job** : int
Specifies which elements of :math:`C` are required.
:math:`\mathrm{job} = -1`
The upper triangular part of :math:`C` is required.
:math:`\mathrm{job} = 0`
The diagonal elements of :math:`C` are required.
:math:`\mathrm{job} > 0`
The elements of column :math:`\mathrm{job}` of :math:`C` are required.
**sigma** : float
The standard error of the residual vector.
**a** : float, array-like, shape :math:`\left(:, p\right)`
Note: the required extent for this argument in dimension 1 is determined as follows: if :math:`\mathrm{not}\left(\mathrm{svd}\right)\text{ or }\mathrm{job}=-1`: :math:`p`; if :math:`\mathrm{svd}= \mathbf{True}\text{ and }\mathrm{job}\geq 0`: :math:`\mathrm{irank}`; otherwise: :math:`0`.
If :math:`\mathrm{svd} = \mathbf{False}`, :math:`\mathrm{a}` must contain the upper triangular matrix :math:`U` of the :math:`QU` factorization of :math:`X`, or of the Cholesky factorization of :math:`X^\mathrm{T}X`; elements of the array below the diagonal need not be set.
If :math:`\mathrm{svd} = \mathbf{True}`, :math:`A` must contain the first :math:`k` rows of the matrix :math:`V^\mathrm{T}`, where :math:`k` is the rank of :math:`X` and :math:`V` is the right-hand orthogonal matrix of the singular value decomposition of :math:`X`.
Thus the :math:`i`\ th row must contain the :math:`i`\ th right-hand singular vector of :math:`X`.
**svd** : bool
Must be :math:`\mathbf{True}` if the least squares solution was obtained from a singular value decomposition of :math:`X`. :math:`\mathrm{svd}` must be :math:`\mathbf{False}` if the least squares solution was obtained from either a :math:`QU` factorization of :math:`X` or a Cholesky factorization of :math:`X^\mathrm{T}X`. In the latter case the rank of :math:`X` is assumed to be :math:`p` and so is applicable only to full rank problems with :math:`n\geq p`.
**irank** : int
If :math:`\mathrm{svd} = \mathbf{True}`, :math:`\mathrm{irank}` must specify the rank :math:`k` of the matrix :math:`X`.
If :math:`\mathrm{svd} = \mathbf{False}`, :math:`\mathrm{irank}` is not referenced and the rank of :math:`X` is assumed to be :math:`p`.
**sv** : float, array-like, shape :math:`\left(p\right)`
If :math:`\mathrm{svd} = \mathbf{True}`, :math:`\mathrm{sv}` must contain the first :math:`\mathrm{irank}` singular values of :math:`X`.
If :math:`\mathrm{svd} = \mathbf{False}`, :math:`\mathrm{sv}` is not referenced.
**Returns**
**a** : float, ndarray, shape :math:`\left(:, p\right)`
If :math:`\mathrm{job}\geq 0`, :math:`A` is unchanged.
If :math:`\mathrm{job} = -1`, :math:`\mathrm{a}` contains the upper triangle of the symmetric matrix :math:`C`.
If :math:`\mathrm{svd} = \mathbf{True}`, elements of the array below the diagonal are used as workspace.
If :math:`\mathrm{svd} = \mathbf{False}`, they are unchanged.
**cj** : float, ndarray, shape :math:`\left(p\right)`
If :math:`\mathrm{job} = 0`, :math:`\mathrm{cj}` returns the diagonal elements of :math:`C`.
If :math:`\mathrm{job} = j > 0`, :math:`\mathrm{cj}` returns the :math:`j`\ th column of :math:`C`.
If :math:`\mathrm{job} = -1`, :math:`\mathrm{cj}` is not referenced.
.. _f04ya-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{lda} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`p = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{svd} = \mathbf{False}`, :math:`\textit{lda}\geq p`.
(`errno` :math:`1`)
On entry, :math:`\textit{lda} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`p = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{svd} = \mathbf{True}` and :math:`\mathrm{job} = -1`, :math:`\textit{lda}\geq p`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{job} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{lda} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{irank} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{svd} = \mathbf{True}` and :math:`\mathrm{job}\geq 0`, :math:`\textit{lda}\geq \mathrm{max}\left(1, \mathrm{irank}\right)`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{irank} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`p = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{svd} = \mathbf{True}`, :math:`\mathrm{irank}\geq 0` and :math:`\mathrm{irank}\leq p`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{job} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`p = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{job}\geq -1` and :math:`\mathrm{job}\leq p`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{sigma} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{sigma}\geq 0.0`.
(`errno` :math:`1`)
On entry, :math:`p = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`p\geq 1`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{irank} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{svd} = \mathbf{True}`, :math:`\mathrm{irank} > 0`.
(`errno` :math:`3`)
On entry, :math:`\mathrm{svd} = \mathbf{False}` and overflow would occur in calculating an element of :math:`C`. The upper triangular matrix :math:`U` must be very nearly singular.
(`errno` :math:`4`)
On entry, :math:`\mathrm{svd} = \mathbf{False}` and one of the first :math:`\mathrm{irank}` singular values is zero. Either the first :math:`\mathrm{irank}` singular values or :math:`\mathrm{irank}` must be incorrect: :math:`\mathrm{irank} = \langle\mathit{\boldsymbol{value}}\rangle`.
.. _f04ya-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
The estimated variance-covariance matrix :math:`C` of the sample regression coefficients is given by
.. math::
C = \sigma^2{\left(X^\mathrm{T}X\right)}^{{-1}},\quad \text{ }\quad X^\mathrm{T}X\quad \text{ nonsingular,}
where :math:`X^\mathrm{T}X` is the normal matrix for the linear least squares regression problem
.. math::
\mathrm{min}:\left\lVert y-Xb\right\rVert_2\text{,}
:math:`\sigma^2` is the estimated variance of the residual vector :math:`r = y-Xb`, and :math:`X` is an :math:`n\times p` observation matrix.
When :math:`X^\mathrm{T}X` is singular, :math:`C` is taken to be
.. math::
C = \sigma^2\left(X^\mathrm{T}X\right)^†\text{,}
where :math:`\left(X^\mathrm{T}X\right)^†` is the pseudo-inverse of :math:`X^\mathrm{T}X`; this assumes that the minimal least squares solution of `(1) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04yaf.html#eqn1>`__ has been found.
The diagonal elements of :math:`C` are the estimated variances of the sample regression coefficients, :math:`b`.
The function can be used to find either the diagonal elements of :math:`C`, or the elements of the :math:`j`\ th column of :math:`C`, or the upper triangular part of :math:`C`.
This function must be preceded by a function that returns either the upper triangular matrix :math:`U` of the :math:`QU` factorization of :math:`X` or of the Cholesky factorization of :math:`X^\mathrm{T}X`, or the singular values and right singular vectors of :math:`X`.
In particular this function can be preceded by one of the functions :meth:`lapackeig.dgelss <naginterfaces.library.lapackeig.dgelss>` or :meth:`real_gen_solve`, which return the arguments :math:`\mathrm{irank}`, :math:`\mathrm{sigma}`, :math:`\mathrm{a}` and :math:`\mathrm{sv}` in the required form. :meth:`real_gen_solve` returns the argument :math:`\mathrm{svd}`, but when this function is used following function :meth:`lapackeig.dgelss <naginterfaces.library.lapackeig.dgelss>` the argument :math:`\mathrm{svd}` should be set to :math:`\mathbf{True}`.
The argument :math:`\textit{p}` of this function corresponds to the argument :math:`\textit{n}` in functions :meth:`lapackeig.dgelss <naginterfaces.library.lapackeig.dgelss>` and :meth:`real_gen_solve`.
.. _f04ya-py2-py-references:
**References**
Anderson, T W, 1958, `An Introduction to Multivariate Statistical Analysis`, Wiley
Lawson, C L and Hanson, R J, 1974, `Solving Least Squares Problems`, Prentice--Hall
"""
raise NotImplementedError
[docs]def real_gen_norm_rcomm(irevcm, x, y, estnrm, seed, comm):
r"""
``real_gen_norm_rcomm`` estimates the :math:`1`-norm of a real rectangular matrix without accessing the matrix explicitly.
It uses reverse communication for evaluating matrix products.
The function may be used for estimating condition numbers of square matrices.
.. _f04yd-py2-py-doc:
For full information please refer to the NAG Library document for f04yd
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04ydf.html
.. _f04yd-py2-py-parameters:
**Parameters**
**irevcm** : int
`On initial entry`: must be set to :math:`0`.
`On intermediate entry`: :math:`\mathrm{irevcm}` must be unchanged.
**x** : float, ndarray, shape :math:`\left(n, t\right)`, modified in place
`On initial entry`: need not be set.
`On intermediate exit`: if :math:`\mathrm{irevcm} = 1`, contains the current matrix :math:`X`.
`On intermediate entry`: if :math:`\mathrm{irevcm} = 2`, must contain :math:`A^\mathrm{T}Y`.
`On final exit`: the array is undefined.
**y** : float, ndarray, shape :math:`\left(m, t\right)`, modified in place
`On initial entry`: need not be set.
`On intermediate exit`: if :math:`\mathrm{irevcm} = 2`, contains the current matrix :math:`Y`.
`On intermediate entry`: if :math:`\mathrm{irevcm} = 1`, must contain :math:`AX`.
`On final exit`: the array is undefined.
**estnrm** : float
`On initial entry`: need not be set.
`On intermediate entry`: must not be changed.
**seed** : int
The seed used for random number generation.
If :math:`t = 1`, :math:`\mathrm{seed}` is not used.
**comm** : dict, communication object, modified in place
Communication structure.
`On initial entry`: need not be set.
**Returns**
**irevcm** : int
`On intermediate exit`: :math:`\mathrm{irevcm} = 1` or :math:`2`, and :math:`\mathrm{x}` contains the :math:`n\times t` matrix :math:`X` and :math:`\mathrm{y}` contains the :math:`m\times t` matrix :math:`Y`. The calling program must
(a) if :math:`\mathrm{irevcm} = 1`, evaluate :math:`AX` and store the result in :math:`\mathrm{y}`
or
if :math:`\mathrm{irevcm} = 2`, evaluate :math:`A^\mathrm{T}Y` and store the result in :math:`\mathrm{x}`,
(#) call ``real_gen_norm_rcomm`` once again, with all the other arguments unchanged.
`On final exit`: :math:`\mathrm{irevcm} = 0`.
**estnrm** : float
`On final exit`: an estimate (a lower bound) for :math:`\left\lVert A\right\rVert_1`.
.. _f04yd-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-10`)
On entry, :math:`t = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{seed} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`t > 1`, :math:`\mathrm{seed}\geq 1`.
(`errno` :math:`-9`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`t = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq t\leq m`.
(`errno` :math:`-3`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{irevcm} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{irevcm} = 0`, :math:`1` or :math:`2`.
(`errno` :math:`-1`)
On initial entry, :math:`\mathrm{irevcm} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{irevcm} = 0`.
(`errno` :math:`1`)
Internal error; please contact `NAG <https://www.nag.com>`__.
.. _f04yd-py2-py-notes:
**Notes**
``real_gen_norm_rcomm`` computes an estimate (a lower bound) for the :math:`1`-norm
.. math::
\left\lVert A\right\rVert_1 = \mathrm{max}_{{1\leq j\leq n}}\sum_{{i = 1}}^m\left\lvert a_{{ij}}\right\rvert
of an :math:`m\times n` real matrix :math:`A = \left(a_{{ij}}\right)`.
The function regards the matrix :math:`A` as being defined by a user-supplied 'Black Box' which, given an :math:`n\times t` matrix :math:`X` (with :math:`t≪n`) or an :math:`m\times t` matrix :math:`Y`, can return :math:`AX` or :math:`A^\mathrm{T}Y`.
A reverse communication interface is used; thus control is returned to the calling program whenever a matrix product is required.
**Note:** this function is **not recommended** for use when the elements of :math:`A` are known explicitly; it is then more efficient to compute the :math:`1`-norm directly from formula `(1) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04ydf.html#eqn1>`__ above.
The **main use** of the function is for estimating :math:`\left\lVert B^{-1}\right\rVert_1` for a square matrix, :math:`B`, and hence the **condition number** :math:`\kappa_1\left(B\right) = \left\lVert B\right\rVert_1\left\lVert B^{-1}\right\rVert_1`, without forming :math:`B^{-1}` explicitly (:math:`A = B^{-1}` above).
If, for example, an :math:`LU` factorization of :math:`B` is available, the matrix products :math:`B^{-1}X` and :math:`B^{-\mathrm{T}}Y` required by ``real_gen_norm_rcomm`` may be computed by back - and forward-substitutions, without computing :math:`B^{-1}`.
The function can also be used to estimate :math:`1`-norms of matrix products such as :math:`A^{-1}B` and :math:`ABC`, without forming the products explicitly.
Further applications are described by Higham (1988).
Since :math:`\left\lVert A\right\rVert_\infty = \left\lVert A^\mathrm{T}\right\rVert_1`, ``real_gen_norm_rcomm`` can be used to estimate the :math:`\infty`-norm of :math:`A` by working with :math:`A^\mathrm{T}` instead of :math:`A`.
The algorithm used is described in Higham and Tisseur (2000).
.. _f04yd-py2-py-references:
**References**
Higham, N J, 1988, `FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation`, ACM Trans. Math. Software (14), 381--396
Higham, N J and Tisseur, F, 2000, `A block algorithm for matrix` :math:`1` `-norm estimation, with an application to` :math:`1` `-norm pseudospectra`, SIAM J. Matrix. Anal. Appl. (21), 1185--1201
"""
raise NotImplementedError
[docs]def complex_gen_norm_rcomm(irevcm, x, y, estnrm, seed, comm):
r"""
``complex_gen_norm_rcomm`` estimates the :math:`1`-norm of a complex rectangular matrix without accessing the matrix explicitly.
It uses reverse communication for evaluating matrix products.
The function may be used for estimating condition numbers of square matrices.
.. _f04zd-py2-py-doc:
For full information please refer to the NAG Library document for f04zd
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04zdf.html
.. _f04zd-py2-py-parameters:
**Parameters**
**irevcm** : int
`On initial entry`: must be set to :math:`0`.
`On intermediate entry`: :math:`\mathrm{irevcm}` must be unchanged.
**x** : complex, ndarray, shape :math:`\left(n, t\right)`, modified in place
`On initial entry`: need not be set.
`On intermediate exit`: if :math:`\mathrm{irevcm} = 1`, contains the current matrix :math:`X`.
`On intermediate entry`: if :math:`\mathrm{irevcm} = 2`, must contain :math:`A^\mathrm{H}Y`.
`On final exit`: the array is undefined.
**y** : complex, ndarray, shape :math:`\left(m, t\right)`, modified in place
`On initial entry`: need not be set.
`On intermediate exit`: if :math:`\mathrm{irevcm} = 2`, contains the current matrix :math:`Y`.
`On intermediate entry`: if :math:`\mathrm{irevcm} = 1`, must contain :math:`AX`.
`On final exit`: the array is undefined.
**estnrm** : float
`On initial entry`: need not be set.
`On intermediate entry`: must not be changed.
**seed** : int
The seed used for random number generation.
If :math:`t = 1`, :math:`\mathrm{seed}` is not used.
**comm** : dict, communication object, modified in place
Communication structure.
`On initial entry`: need not be set.
**Returns**
**irevcm** : int
`On intermediate exit`: :math:`\mathrm{irevcm} = 1` or :math:`2`, and :math:`\mathrm{x}` contains the :math:`n\times t` matrix :math:`X` and :math:`\mathrm{y}` contains the :math:`m\times t` matrix :math:`Y`. The calling program must
(a) if :math:`\mathrm{irevcm} = 1`, evaluate :math:`AX` and store the result in :math:`\mathrm{y}`
or
if :math:`\mathrm{irevcm} = 2`, evaluate :math:`A^\mathrm{H}Y` and store the result in :math:`\mathrm{x}`, where :math:`A^\mathrm{H}` is the complex conjugate transpose;
(#) call ``complex_gen_norm_rcomm`` once again, with all the arguments unchanged.
`On final exit`: :math:`\mathrm{irevcm} = 0`.
**estnrm** : float
`On final exit`: an estimate (a lower bound) for :math:`\left\lVert A\right\rVert_1`.
.. _f04zd-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`-10`)
On entry, :math:`t = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{seed} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`t > 1`, :math:`\mathrm{seed}\geq 1`.
(`errno` :math:`-9`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`t = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq t\leq m`.
(`errno` :math:`-3`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 0`.
(`errno` :math:`-2`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m\geq 0`.
(`errno` :math:`-1`)
On entry, :math:`\mathrm{irevcm} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{irevcm} = 0`, :math:`1` or :math:`2`.
(`errno` :math:`-1`)
On initial entry, :math:`\mathrm{irevcm} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{irevcm} = 0`.
(`errno` :math:`1`)
Internal error; please contact `NAG <https://www.nag.com>`__.
.. _f04zd-py2-py-notes:
**Notes**
``complex_gen_norm_rcomm`` computes an estimate (a lower bound) for the :math:`1`-norm
.. math::
\left\lVert A\right\rVert_1 = \mathrm{max}_{{1\leq j\leq n}}\sum_{{i = 1}}^m\left\lvert a_{{ij}}\right\rvert
of an :math:`m\times n` complex matrix :math:`A = \left(a_{{ij}}\right)`.
The function regards the matrix :math:`A` as being defined by a user-supplied 'Black Box' which, given an :math:`n\times t` matrix :math:`X` (with :math:`t≪n`) or an :math:`m\times t` matrix :math:`Y`, can return :math:`AX` or :math:`A^\mathrm{H}Y`, where :math:`A^\mathrm{H}` is the complex conjugate transpose.
A reverse communication interface is used; thus control is returned to the calling program whenever a matrix product is required.
**Note:** this function is **not** **recommended** for use when the elements of :math:`A` are known explicitly; it is then more efficient to compute the :math:`1`-norm directly from the formula `(1) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/f04/f04zdf.html#eqn1>`__ above.
The **main use** of the function is for estimating :math:`\left\lVert B^{-1}\right\rVert_1` for a square matrix :math:`B`, and hence the **condition number** :math:`\kappa_1\left(B\right) = \left\lVert B\right\rVert_1\left\lVert B^{-1}\right\rVert_1`, without forming :math:`B^{-1}` explicitly (:math:`A = B^{-1}` above).
If, for example, an :math:`LU` factorization of :math:`B` is available, the matrix products :math:`B^{-1}X` and :math:`B^{-\mathrm{H}}Y` required by ``complex_gen_norm_rcomm`` may be computed by back - and forward-substitutions, without computing :math:`B^{-1}`.
The function can also be used to estimate :math:`1`-norms of matrix products such as :math:`A^{-1}B` and :math:`ABC`, without forming the products explicitly.
Further applications are described in Higham (1988).
Since :math:`\left\lVert A\right\rVert_\infty = \left\lVert A^\mathrm{H}\right\rVert_1`, ``complex_gen_norm_rcomm`` can be used to estimate the :math:`\infty`-norm of :math:`A` by working with :math:`A^\mathrm{H}` instead of :math:`A`.
The algorithm used is described in Higham and Tisseur (2000).
.. _f04zd-py2-py-references:
**References**
Higham, N J, 1988, `FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation`, ACM Trans. Math. Software (14), 381--396
Higham, N J and Tisseur, F, 2000, `A block algorithm for matrix` :math:`1` `-norm estimation, with an application to` :math:`1` `-norm pseudospectra`, SIAM J. Matrix. Anal. Appl. (21), 1185--1201
"""
raise NotImplementedError
