Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_decomp_schur.py: 17%
102 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Schur decomposition functions."""
2import numpy
3from numpy import asarray_chkfinite, single, asarray, array
4from numpy.linalg import norm
7# Local imports.
8from ._misc import LinAlgError, _datacopied
9from .lapack import get_lapack_funcs
10from ._decomp import eigvals
12__all__ = ['schur', 'rsf2csf']
14_double_precision = ['i', 'l', 'd']
17def schur(a, output='real', lwork=None, overwrite_a=False, sort=None,
18 check_finite=True):
19 """
20 Compute Schur decomposition of a matrix.
22 The Schur decomposition is::
24 A = Z T Z^H
26 where Z is unitary and T is either upper-triangular, or for real
27 Schur decomposition (output='real'), quasi-upper triangular. In
28 the quasi-triangular form, 2x2 blocks describing complex-valued
29 eigenvalue pairs may extrude from the diagonal.
31 Parameters
32 ----------
33 a : (M, M) array_like
34 Matrix to decompose
35 output : {'real', 'complex'}, optional
36 Construct the real or complex Schur decomposition (for real matrices).
37 lwork : int, optional
38 Work array size. If None or -1, it is automatically computed.
39 overwrite_a : bool, optional
40 Whether to overwrite data in a (may improve performance).
41 sort : {None, callable, 'lhp', 'rhp', 'iuc', 'ouc'}, optional
42 Specifies whether the upper eigenvalues should be sorted. A callable
43 may be passed that, given a eigenvalue, returns a boolean denoting
44 whether the eigenvalue should be sorted to the top-left (True).
45 Alternatively, string parameters may be used::
47 'lhp' Left-hand plane (x.real < 0.0)
48 'rhp' Right-hand plane (x.real > 0.0)
49 'iuc' Inside the unit circle (x*x.conjugate() <= 1.0)
50 'ouc' Outside the unit circle (x*x.conjugate() > 1.0)
52 Defaults to None (no sorting).
53 check_finite : bool, optional
54 Whether to check that the input matrix contains only finite numbers.
55 Disabling may give a performance gain, but may result in problems
56 (crashes, non-termination) if the inputs do contain infinities or NaNs.
58 Returns
59 -------
60 T : (M, M) ndarray
61 Schur form of A. It is real-valued for the real Schur decomposition.
62 Z : (M, M) ndarray
63 An unitary Schur transformation matrix for A.
64 It is real-valued for the real Schur decomposition.
65 sdim : int
66 If and only if sorting was requested, a third return value will
67 contain the number of eigenvalues satisfying the sort condition.
69 Raises
70 ------
71 LinAlgError
72 Error raised under three conditions:
74 1. The algorithm failed due to a failure of the QR algorithm to
75 compute all eigenvalues.
76 2. If eigenvalue sorting was requested, the eigenvalues could not be
77 reordered due to a failure to separate eigenvalues, usually because
78 of poor conditioning.
79 3. If eigenvalue sorting was requested, roundoff errors caused the
80 leading eigenvalues to no longer satisfy the sorting condition.
82 See Also
83 --------
84 rsf2csf : Convert real Schur form to complex Schur form
86 Examples
87 --------
88 >>> import numpy as np
89 >>> from scipy.linalg import schur, eigvals
90 >>> A = np.array([[0, 2, 2], [0, 1, 2], [1, 0, 1]])
91 >>> T, Z = schur(A)
92 >>> T
93 array([[ 2.65896708, 1.42440458, -1.92933439],
94 [ 0. , -0.32948354, -0.49063704],
95 [ 0. , 1.31178921, -0.32948354]])
96 >>> Z
97 array([[0.72711591, -0.60156188, 0.33079564],
98 [0.52839428, 0.79801892, 0.28976765],
99 [0.43829436, 0.03590414, -0.89811411]])
101 >>> T2, Z2 = schur(A, output='complex')
102 >>> T2
103 array([[ 2.65896708, -1.22839825+1.32378589j, 0.42590089+1.51937378j],
104 [ 0. , -0.32948354+0.80225456j, -0.59877807+0.56192146j],
105 [ 0. , 0. , -0.32948354-0.80225456j]])
106 >>> eigvals(T2)
107 array([2.65896708, -0.32948354+0.80225456j, -0.32948354-0.80225456j])
109 An arbitrary custom eig-sorting condition, having positive imaginary part,
110 which is satisfied by only one eigenvalue
112 >>> T3, Z3, sdim = schur(A, output='complex', sort=lambda x: x.imag > 0)
113 >>> sdim
114 1
116 """
117 if output not in ['real', 'complex', 'r', 'c']:
118 raise ValueError("argument must be 'real', or 'complex'")
119 if check_finite:
120 a1 = asarray_chkfinite(a)
121 else:
122 a1 = asarray(a)
123 if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
124 raise ValueError('expected square matrix')
125 typ = a1.dtype.char
126 if output in ['complex', 'c'] and typ not in ['F', 'D']:
127 if typ in _double_precision:
128 a1 = a1.astype('D')
129 typ = 'D'
130 else:
131 a1 = a1.astype('F')
132 typ = 'F'
133 overwrite_a = overwrite_a or (_datacopied(a1, a))
134 gees, = get_lapack_funcs(('gees',), (a1,))
135 if lwork is None or lwork == -1:
136 # get optimal work array
137 result = gees(lambda x: None, a1, lwork=-1)
138 lwork = result[-2][0].real.astype(numpy.int_)
140 if sort is None:
141 sort_t = 0
142 sfunction = lambda x: None
143 else:
144 sort_t = 1
145 if callable(sort):
146 sfunction = sort
147 elif sort == 'lhp':
148 sfunction = lambda x: (x.real < 0.0)
149 elif sort == 'rhp':
150 sfunction = lambda x: (x.real >= 0.0)
151 elif sort == 'iuc':
152 sfunction = lambda x: (abs(x) <= 1.0)
153 elif sort == 'ouc':
154 sfunction = lambda x: (abs(x) > 1.0)
155 else:
156 raise ValueError("'sort' parameter must either be 'None', or a "
157 "callable, or one of ('lhp','rhp','iuc','ouc')")
159 result = gees(sfunction, a1, lwork=lwork, overwrite_a=overwrite_a,
160 sort_t=sort_t)
162 info = result[-1]
163 if info < 0:
164 raise ValueError('illegal value in {}-th argument of internal gees'
165 ''.format(-info))
166 elif info == a1.shape[0] + 1:
167 raise LinAlgError('Eigenvalues could not be separated for reordering.')
168 elif info == a1.shape[0] + 2:
169 raise LinAlgError('Leading eigenvalues do not satisfy sort condition.')
170 elif info > 0:
171 raise LinAlgError("Schur form not found. Possibly ill-conditioned.")
173 if sort_t == 0:
174 return result[0], result[-3]
175 else:
176 return result[0], result[-3], result[1]
179eps = numpy.finfo(float).eps
180feps = numpy.finfo(single).eps
182_array_kind = {'b': 0, 'h': 0, 'B': 0, 'i': 0, 'l': 0,
183 'f': 0, 'd': 0, 'F': 1, 'D': 1}
184_array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
185_array_type = [['f', 'd'], ['F', 'D']]
188def _commonType(*arrays):
189 kind = 0
190 precision = 0
191 for a in arrays:
192 t = a.dtype.char
193 kind = max(kind, _array_kind[t])
194 precision = max(precision, _array_precision[t])
195 return _array_type[kind][precision]
198def _castCopy(type, *arrays):
199 cast_arrays = ()
200 for a in arrays:
201 if a.dtype.char == type:
202 cast_arrays = cast_arrays + (a.copy(),)
203 else:
204 cast_arrays = cast_arrays + (a.astype(type),)
205 if len(cast_arrays) == 1:
206 return cast_arrays[0]
207 else:
208 return cast_arrays
211def rsf2csf(T, Z, check_finite=True):
212 """
213 Convert real Schur form to complex Schur form.
215 Convert a quasi-diagonal real-valued Schur form to the upper-triangular
216 complex-valued Schur form.
218 Parameters
219 ----------
220 T : (M, M) array_like
221 Real Schur form of the original array
222 Z : (M, M) array_like
223 Schur transformation matrix
224 check_finite : bool, optional
225 Whether to check that the input arrays contain only finite numbers.
226 Disabling may give a performance gain, but may result in problems
227 (crashes, non-termination) if the inputs do contain infinities or NaNs.
229 Returns
230 -------
231 T : (M, M) ndarray
232 Complex Schur form of the original array
233 Z : (M, M) ndarray
234 Schur transformation matrix corresponding to the complex form
236 See Also
237 --------
238 schur : Schur decomposition of an array
240 Examples
241 --------
242 >>> import numpy as np
243 >>> from scipy.linalg import schur, rsf2csf
244 >>> A = np.array([[0, 2, 2], [0, 1, 2], [1, 0, 1]])
245 >>> T, Z = schur(A)
246 >>> T
247 array([[ 2.65896708, 1.42440458, -1.92933439],
248 [ 0. , -0.32948354, -0.49063704],
249 [ 0. , 1.31178921, -0.32948354]])
250 >>> Z
251 array([[0.72711591, -0.60156188, 0.33079564],
252 [0.52839428, 0.79801892, 0.28976765],
253 [0.43829436, 0.03590414, -0.89811411]])
254 >>> T2 , Z2 = rsf2csf(T, Z)
255 >>> T2
256 array([[2.65896708+0.j, -1.64592781+0.743164187j, -1.21516887+1.00660462j],
257 [0.+0.j , -0.32948354+8.02254558e-01j, -0.82115218-2.77555756e-17j],
258 [0.+0.j , 0.+0.j, -0.32948354-0.802254558j]])
259 >>> Z2
260 array([[0.72711591+0.j, 0.28220393-0.31385693j, 0.51319638-0.17258824j],
261 [0.52839428+0.j, 0.24720268+0.41635578j, -0.68079517-0.15118243j],
262 [0.43829436+0.j, -0.76618703+0.01873251j, -0.03063006+0.46857912j]])
264 """
265 if check_finite:
266 Z, T = map(asarray_chkfinite, (Z, T))
267 else:
268 Z, T = map(asarray, (Z, T))
270 for ind, X in enumerate([Z, T]):
271 if X.ndim != 2 or X.shape[0] != X.shape[1]:
272 raise ValueError("Input '{}' must be square.".format('ZT'[ind]))
274 if T.shape[0] != Z.shape[0]:
275 raise ValueError("Input array shapes must match: Z: {} vs. T: {}"
276 "".format(Z.shape, T.shape))
277 N = T.shape[0]
278 t = _commonType(Z, T, array([3.0], 'F'))
279 Z, T = _castCopy(t, Z, T)
281 for m in range(N-1, 0, -1):
282 if abs(T[m, m-1]) > eps*(abs(T[m-1, m-1]) + abs(T[m, m])):
283 mu = eigvals(T[m-1:m+1, m-1:m+1]) - T[m, m]
284 r = norm([mu[0], T[m, m-1]])
285 c = mu[0] / r
286 s = T[m, m-1] / r
287 G = array([[c.conj(), s], [-s, c]], dtype=t)
289 T[m-1:m+1, m-1:] = G.dot(T[m-1:m+1, m-1:])
290 T[:m+1, m-1:m+1] = T[:m+1, m-1:m+1].dot(G.conj().T)
291 Z[:, m-1:m+1] = Z[:, m-1:m+1].dot(G.conj().T)
293 T[m, m-1] = 0.0
294 return T, Z