Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/linalg/_decomp_schur.py: 15%
110 statements
« prev ^ index » next coverage.py v7.4.1, created at 2024-02-14 06:37 +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], # may vary
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 numpy.issubdtype(a1.dtype, numpy.integer):
124 a1 = asarray(a, dtype=numpy.dtype("long"))
125 if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
126 raise ValueError('expected square matrix')
127 typ = a1.dtype.char
128 if output in ['complex', 'c'] and typ not in ['F', 'D']:
129 if typ in _double_precision:
130 a1 = a1.astype('D')
131 typ = 'D'
132 else:
133 a1 = a1.astype('F')
134 typ = 'F'
135 overwrite_a = overwrite_a or (_datacopied(a1, a))
136 gees, = get_lapack_funcs(('gees',), (a1,))
137 if lwork is None or lwork == -1:
138 # get optimal work array
139 result = gees(lambda x: None, a1, lwork=-1)
140 lwork = result[-2][0].real.astype(numpy.int_)
142 if sort is None:
143 sort_t = 0
144 def sfunction(x):
145 return None
146 else:
147 sort_t = 1
148 if callable(sort):
149 sfunction = sort
150 elif sort == 'lhp':
151 def sfunction(x):
152 return x.real < 0.0
153 elif sort == 'rhp':
154 def sfunction(x):
155 return x.real >= 0.0
156 elif sort == 'iuc':
157 def sfunction(x):
158 return abs(x) <= 1.0
159 elif sort == 'ouc':
160 def sfunction(x):
161 return abs(x) > 1.0
162 else:
163 raise ValueError("'sort' parameter must either be 'None', or a "
164 "callable, or one of ('lhp','rhp','iuc','ouc')")
166 result = gees(sfunction, a1, lwork=lwork, overwrite_a=overwrite_a,
167 sort_t=sort_t)
169 info = result[-1]
170 if info < 0:
171 raise ValueError(f'illegal value in {-info}-th argument of internal gees')
172 elif info == a1.shape[0] + 1:
173 raise LinAlgError('Eigenvalues could not be separated for reordering.')
174 elif info == a1.shape[0] + 2:
175 raise LinAlgError('Leading eigenvalues do not satisfy sort condition.')
176 elif info > 0:
177 raise LinAlgError("Schur form not found. Possibly ill-conditioned.")
179 if sort_t == 0:
180 return result[0], result[-3]
181 else:
182 return result[0], result[-3], result[1]
185eps = numpy.finfo(float).eps
186feps = numpy.finfo(single).eps
188_array_kind = {'b': 0, 'h': 0, 'B': 0, 'i': 0, 'l': 0,
189 'f': 0, 'd': 0, 'F': 1, 'D': 1}
190_array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
191_array_type = [['f', 'd'], ['F', 'D']]
194def _commonType(*arrays):
195 kind = 0
196 precision = 0
197 for a in arrays:
198 t = a.dtype.char
199 kind = max(kind, _array_kind[t])
200 precision = max(precision, _array_precision[t])
201 return _array_type[kind][precision]
204def _castCopy(type, *arrays):
205 cast_arrays = ()
206 for a in arrays:
207 if a.dtype.char == type:
208 cast_arrays = cast_arrays + (a.copy(),)
209 else:
210 cast_arrays = cast_arrays + (a.astype(type),)
211 if len(cast_arrays) == 1:
212 return cast_arrays[0]
213 else:
214 return cast_arrays
217def rsf2csf(T, Z, check_finite=True):
218 """
219 Convert real Schur form to complex Schur form.
221 Convert a quasi-diagonal real-valued Schur form to the upper-triangular
222 complex-valued Schur form.
224 Parameters
225 ----------
226 T : (M, M) array_like
227 Real Schur form of the original array
228 Z : (M, M) array_like
229 Schur transformation matrix
230 check_finite : bool, optional
231 Whether to check that the input arrays contain only finite numbers.
232 Disabling may give a performance gain, but may result in problems
233 (crashes, non-termination) if the inputs do contain infinities or NaNs.
235 Returns
236 -------
237 T : (M, M) ndarray
238 Complex Schur form of the original array
239 Z : (M, M) ndarray
240 Schur transformation matrix corresponding to the complex form
242 See Also
243 --------
244 schur : Schur decomposition of an array
246 Examples
247 --------
248 >>> import numpy as np
249 >>> from scipy.linalg import schur, rsf2csf
250 >>> A = np.array([[0, 2, 2], [0, 1, 2], [1, 0, 1]])
251 >>> T, Z = schur(A)
252 >>> T
253 array([[ 2.65896708, 1.42440458, -1.92933439],
254 [ 0. , -0.32948354, -0.49063704],
255 [ 0. , 1.31178921, -0.32948354]])
256 >>> Z
257 array([[0.72711591, -0.60156188, 0.33079564],
258 [0.52839428, 0.79801892, 0.28976765],
259 [0.43829436, 0.03590414, -0.89811411]])
260 >>> T2 , Z2 = rsf2csf(T, Z)
261 >>> T2
262 array([[2.65896708+0.j, -1.64592781+0.743164187j, -1.21516887+1.00660462j],
263 [0.+0.j , -0.32948354+8.02254558e-01j, -0.82115218-2.77555756e-17j],
264 [0.+0.j , 0.+0.j, -0.32948354-0.802254558j]])
265 >>> Z2
266 array([[0.72711591+0.j, 0.28220393-0.31385693j, 0.51319638-0.17258824j],
267 [0.52839428+0.j, 0.24720268+0.41635578j, -0.68079517-0.15118243j],
268 [0.43829436+0.j, -0.76618703+0.01873251j, -0.03063006+0.46857912j]])
270 """
271 if check_finite:
272 Z, T = map(asarray_chkfinite, (Z, T))
273 else:
274 Z, T = map(asarray, (Z, T))
276 for ind, X in enumerate([Z, T]):
277 if X.ndim != 2 or X.shape[0] != X.shape[1]:
278 raise ValueError("Input '{}' must be square.".format('ZT'[ind]))
280 if T.shape[0] != Z.shape[0]:
281 message = f"Input array shapes must match: Z: {Z.shape} vs. T: {T.shape}"
282 raise ValueError(message)
283 N = T.shape[0]
284 t = _commonType(Z, T, array([3.0], 'F'))
285 Z, T = _castCopy(t, Z, T)
287 for m in range(N-1, 0, -1):
288 if abs(T[m, m-1]) > eps*(abs(T[m-1, m-1]) + abs(T[m, m])):
289 mu = eigvals(T[m-1:m+1, m-1:m+1]) - T[m, m]
290 r = norm([mu[0], T[m, m-1]])
291 c = mu[0] / r
292 s = T[m, m-1] / r
293 G = array([[c.conj(), s], [-s, c]], dtype=t)
295 T[m-1:m+1, m-1:] = G.dot(T[m-1:m+1, m-1:])
296 T[:m+1, m-1:m+1] = T[:m+1, m-1:m+1].dot(G.conj().T)
297 Z[:, m-1:m+1] = Z[:, m-1:m+1].dot(G.conj().T)
299 T[m, m-1] = 0.0
300 return T, Z