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