Coverage for /usr/lib/python3/dist-packages/mpmath/matrices/eigen.py: 6%
352 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1#!/usr/bin/python
2# -*- coding: utf-8 -*-
4##################################################################################################
5# module for the eigenvalue problem
6# Copyright 2013 Timo Hartmann (thartmann15 at gmail.com)
7#
8# todo:
9# - implement balancing
10# - agressive early deflation
11#
12##################################################################################################
14"""
15The eigenvalue problem
16----------------------
18This file contains routines for the eigenvalue problem.
20high level routines:
22 hessenberg : reduction of a real or complex square matrix to upper Hessenberg form
23 schur : reduction of a real or complex square matrix to upper Schur form
24 eig : eigenvalues and eigenvectors of a real or complex square matrix
26low level routines:
28 hessenberg_reduce_0 : reduction of a real or complex square matrix to upper Hessenberg form
29 hessenberg_reduce_1 : auxiliary routine to hessenberg_reduce_0
30 qr_step : a single implicitly shifted QR step for an upper Hessenberg matrix
31 hessenberg_qr : Schur decomposition of an upper Hessenberg matrix
32 eig_tr_r : right eigenvectors of an upper triangular matrix
33 eig_tr_l : left eigenvectors of an upper triangular matrix
34"""
36from ..libmp.backend import xrange
38class Eigen(object):
39 pass
41def defun(f):
42 setattr(Eigen, f.__name__, f)
43 return f
45def hessenberg_reduce_0(ctx, A, T):
46 """
47 This routine computes the (upper) Hessenberg decomposition of a square matrix A.
48 Given A, an unitary matrix Q is calculated such that
50 Q' A Q = H and Q' Q = Q Q' = 1
52 where H is an upper Hessenberg matrix, meaning that it only contains zeros
53 below the first subdiagonal. Here ' denotes the hermitian transpose (i.e.
54 transposition and conjugation).
56 parameters:
57 A (input/output) On input, A contains the square matrix A of
58 dimension (n,n). On output, A contains a compressed representation
59 of Q and H.
60 T (output) An array of length n containing the first elements of
61 the Householder reflectors.
62 """
64 # internally we work with householder reflections from the right.
65 # let u be a row vector (i.e. u[i]=A[i,:i]). then
66 # Q is build up by reflectors of the type (1-v'v) where v is a suitable
67 # modification of u. these reflectors are applyed to A from the right.
68 # because we work with reflectors from the right we have to start with
69 # the bottom row of A and work then upwards (this corresponds to
70 # some kind of RQ decomposition).
71 # the first part of the vectors v (i.e. A[i,:(i-1)]) are stored as row vectors
72 # in the lower left part of A (excluding the diagonal and subdiagonal).
73 # the last entry of v is stored in T.
74 # the upper right part of A (including diagonal and subdiagonal) becomes H.
77 n = A.rows
78 if n <= 2: return
80 for i in xrange(n-1, 1, -1):
82 # scale the vector
84 scale = 0
85 for k in xrange(0, i):
86 scale += abs(ctx.re(A[i,k])) + abs(ctx.im(A[i,k]))
88 scale_inv = 0
89 if scale != 0:
90 scale_inv = 1 / scale
92 if scale == 0 or ctx.isinf(scale_inv):
93 # sadly there are floating point numbers not equal to zero whose reciprocal is infinity
94 T[i] = 0
95 A[i,i-1] = 0
96 continue
98 # calculate parameters for housholder transformation
100 H = 0
101 for k in xrange(0, i):
102 A[i,k] *= scale_inv
103 rr = ctx.re(A[i,k])
104 ii = ctx.im(A[i,k])
105 H += rr * rr + ii * ii
107 F = A[i,i-1]
108 f = abs(F)
109 G = ctx.sqrt(H)
110 A[i,i-1] = - G * scale
112 if f == 0:
113 T[i] = G
114 else:
115 ff = F / f
116 T[i] = F + G * ff
117 A[i,i-1] *= ff
119 H += G * f
120 H = 1 / ctx.sqrt(H)
122 T[i] *= H
123 for k in xrange(0, i - 1):
124 A[i,k] *= H
126 for j in xrange(0, i):
127 # apply housholder transformation (from right)
129 G = ctx.conj(T[i]) * A[j,i-1]
130 for k in xrange(0, i-1):
131 G += ctx.conj(A[i,k]) * A[j,k]
133 A[j,i-1] -= G * T[i]
134 for k in xrange(0, i-1):
135 A[j,k] -= G * A[i,k]
137 for j in xrange(0, n):
138 # apply housholder transformation (from left)
140 G = T[i] * A[i-1,j]
141 for k in xrange(0, i-1):
142 G += A[i,k] * A[k,j]
144 A[i-1,j] -= G * ctx.conj(T[i])
145 for k in xrange(0, i-1):
146 A[k,j] -= G * ctx.conj(A[i,k])
150def hessenberg_reduce_1(ctx, A, T):
151 """
152 This routine forms the unitary matrix Q described in hessenberg_reduce_0.
154 parameters:
155 A (input/output) On input, A is the same matrix as delivered by
156 hessenberg_reduce_0. On output, A is set to Q.
158 T (input) On input, T is the same array as delivered by hessenberg_reduce_0.
159 """
161 n = A.rows
163 if n == 1:
164 A[0,0] = 1
165 return
167 A[0,0] = A[1,1] = 1
168 A[0,1] = A[1,0] = 0
170 for i in xrange(2, n):
171 if T[i] != 0:
173 for j in xrange(0, i):
174 G = T[i] * A[i-1,j]
175 for k in xrange(0, i-1):
176 G += A[i,k] * A[k,j]
178 A[i-1,j] -= G * ctx.conj(T[i])
179 for k in xrange(0, i-1):
180 A[k,j] -= G * ctx.conj(A[i,k])
182 A[i,i] = 1
183 for j in xrange(0, i):
184 A[j,i] = A[i,j] = 0
188@defun
189def hessenberg(ctx, A, overwrite_a = False):
190 """
191 This routine computes the Hessenberg decomposition of a square matrix A.
192 Given A, an unitary matrix Q is determined such that
194 Q' A Q = H and Q' Q = Q Q' = 1
196 where H is an upper right Hessenberg matrix. Here ' denotes the hermitian
197 transpose (i.e. transposition and conjugation).
199 input:
200 A : a real or complex square matrix
201 overwrite_a : if true, allows modification of A which may improve
202 performance. if false, A is not modified.
204 output:
205 Q : an unitary matrix
206 H : an upper right Hessenberg matrix
208 example:
209 >>> from mpmath import mp
210 >>> A = mp.matrix([[3, -1, 2], [2, 5, -5], [-2, -3, 7]])
211 >>> Q, H = mp.hessenberg(A)
212 >>> mp.nprint(H, 3) # doctest:+SKIP
213 [ 3.15 2.23 4.44]
214 [-0.769 4.85 3.05]
215 [ 0.0 3.61 7.0]
216 >>> print(mp.chop(A - Q * H * Q.transpose_conj()))
217 [0.0 0.0 0.0]
218 [0.0 0.0 0.0]
219 [0.0 0.0 0.0]
221 return value: (Q, H)
222 """
224 n = A.rows
226 if n == 1:
227 return (ctx.matrix([[1]]), A)
229 if not overwrite_a:
230 A = A.copy()
232 T = ctx.matrix(n, 1)
234 hessenberg_reduce_0(ctx, A, T)
235 Q = A.copy()
236 hessenberg_reduce_1(ctx, Q, T)
238 for x in xrange(n):
239 for y in xrange(x+2, n):
240 A[y,x] = 0
242 return Q, A
245###########################################################################
248def qr_step(ctx, n0, n1, A, Q, shift):
249 """
250 This subroutine executes a single implicitly shifted QR step applied to an
251 upper Hessenberg matrix A. Given A and shift as input, first an QR
252 decomposition is calculated:
254 Q R = A - shift * 1 .
256 The output is then following matrix:
258 R Q + shift * 1
260 parameters:
261 n0, n1 (input) Two integers which specify the submatrix A[n0:n1,n0:n1]
262 on which this subroutine operators. The subdiagonal elements
263 to the left and below this submatrix must be deflated (i.e. zero).
264 following restriction is imposed: n1>=n0+2
265 A (input/output) On input, A is an upper Hessenberg matrix.
266 On output, A is replaced by "R Q + shift * 1"
267 Q (input/output) The parameter Q is multiplied by the unitary matrix
268 Q arising from the QR decomposition. Q can also be false, in which
269 case the unitary matrix Q is not computated.
270 shift (input) a complex number specifying the shift. idealy close to an
271 eigenvalue of the bottemmost part of the submatrix A[n0:n1,n0:n1].
273 references:
274 Stoer, Bulirsch - Introduction to Numerical Analysis.
275 Kresser : Numerical Methods for General and Structured Eigenvalue Problems
276 """
278 # implicitly shifted and bulge chasing is explained at p.398/399 in "Stoer, Bulirsch - Introduction to Numerical Analysis"
279 # for bulge chasing see also "Watkins - The Matrix Eigenvalue Problem" sec.4.5,p.173
281 # the Givens rotation we used is determined as follows: let c,s be two complex
282 # numbers. then we have following relation:
283 #
284 # v = sqrt(|c|^2 + |s|^2)
285 #
286 # 1/v [ c~ s~] [c] = [v]
287 # [-s c ] [s] [0]
288 #
289 # the matrix on the left is our Givens rotation.
291 n = A.rows
293 # first step
295 # calculate givens rotation
296 c = A[n0 ,n0] - shift
297 s = A[n0+1,n0]
299 v = ctx.hypot(ctx.hypot(ctx.re(c), ctx.im(c)), ctx.hypot(ctx.re(s), ctx.im(s)))
301 if v == 0:
302 v = 1
303 c = 1
304 s = 0
305 else:
306 c /= v
307 s /= v
309 cc = ctx.conj(c)
310 cs = ctx.conj(s)
312 for k in xrange(n0, n):
313 # apply givens rotation from the left
314 x = A[n0 ,k]
315 y = A[n0+1,k]
316 A[n0 ,k] = cc * x + cs * y
317 A[n0+1,k] = c * y - s * x
319 for k in xrange(min(n1, n0+3)):
320 # apply givens rotation from the right
321 x = A[k,n0 ]
322 y = A[k,n0+1]
323 A[k,n0 ] = c * x + s * y
324 A[k,n0+1] = cc * y - cs * x
326 if not isinstance(Q, bool):
327 for k in xrange(n):
328 # eigenvectors
329 x = Q[k,n0 ]
330 y = Q[k,n0+1]
331 Q[k,n0 ] = c * x + s * y
332 Q[k,n0+1] = cc * y - cs * x
334 # chase the bulge
336 for j in xrange(n0, n1 - 2):
337 # calculate givens rotation
339 c = A[j+1,j]
340 s = A[j+2,j]
342 v = ctx.hypot(ctx.hypot(ctx.re(c), ctx.im(c)), ctx.hypot(ctx.re(s), ctx.im(s)))
344 if v == 0:
345 A[j+1,j] = 0
346 v = 1
347 c = 1
348 s = 0
349 else:
350 A[j+1,j] = v
351 c /= v
352 s /= v
354 A[j+2,j] = 0
356 cc = ctx.conj(c)
357 cs = ctx.conj(s)
359 for k in xrange(j+1, n):
360 # apply givens rotation from the left
361 x = A[j+1,k]
362 y = A[j+2,k]
363 A[j+1,k] = cc * x + cs * y
364 A[j+2,k] = c * y - s * x
366 for k in xrange(0, min(n1, j+4)):
367 # apply givens rotation from the right
368 x = A[k,j+1]
369 y = A[k,j+2]
370 A[k,j+1] = c * x + s * y
371 A[k,j+2] = cc * y - cs * x
373 if not isinstance(Q, bool):
374 for k in xrange(0, n):
375 # eigenvectors
376 x = Q[k,j+1]
377 y = Q[k,j+2]
378 Q[k,j+1] = c * x + s * y
379 Q[k,j+2] = cc * y - cs * x
383def hessenberg_qr(ctx, A, Q):
384 """
385 This routine computes the Schur decomposition of an upper Hessenberg matrix A.
386 Given A, an unitary matrix Q is determined such that
388 Q' A Q = R and Q' Q = Q Q' = 1
390 where R is an upper right triangular matrix. Here ' denotes the hermitian
391 transpose (i.e. transposition and conjugation).
393 parameters:
394 A (input/output) On input, A contains an upper Hessenberg matrix.
395 On output, A is replace by the upper right triangluar matrix R.
397 Q (input/output) The parameter Q is multiplied by the unitary
398 matrix Q arising from the Schur decomposition. Q can also be
399 false, in which case the unitary matrix Q is not computated.
400 """
402 n = A.rows
404 norm = 0
405 for x in xrange(n):
406 for y in xrange(min(x+2, n)):
407 norm += ctx.re(A[y,x]) ** 2 + ctx.im(A[y,x]) ** 2
408 norm = ctx.sqrt(norm) / n
410 if norm == 0:
411 return
413 n0 = 0
414 n1 = n
416 eps = ctx.eps / (100 * n)
417 maxits = ctx.dps * 4
419 its = totalits = 0
421 while 1:
422 # kressner p.32 algo 3
423 # the active submatrix is A[n0:n1,n0:n1]
425 k = n0
427 while k + 1 < n1:
428 s = abs(ctx.re(A[k,k])) + abs(ctx.im(A[k,k])) + abs(ctx.re(A[k+1,k+1])) + abs(ctx.im(A[k+1,k+1]))
429 if s < eps * norm:
430 s = norm
431 if abs(A[k+1,k]) < eps * s:
432 break
433 k += 1
435 if k + 1 < n1:
436 # deflation found at position (k+1, k)
438 A[k+1,k] = 0
439 n0 = k + 1
441 its = 0
443 if n0 + 1 >= n1:
444 # block of size at most two has converged
445 n0 = 0
446 n1 = k + 1
447 if n1 < 2:
448 # QR algorithm has converged
449 return
450 else:
451 if (its % 30) == 10:
452 # exceptional shift
453 shift = A[n1-1,n1-2]
454 elif (its % 30) == 20:
455 # exceptional shift
456 shift = abs(A[n1-1,n1-2])
457 elif (its % 30) == 29:
458 # exceptional shift
459 shift = norm
460 else:
461 # A = [ a b ] det(x-A)=x*x-x*tr(A)+det(A)
462 # [ c d ]
463 #
464 # eigenvalues bad: (tr(A)+sqrt((tr(A))**2-4*det(A)))/2
465 # bad because of cancellation if |c| is small and |a-d| is small, too.
466 #
467 # eigenvalues good: (a+d+sqrt((a-d)**2+4*b*c))/2
469 t = A[n1-2,n1-2] + A[n1-1,n1-1]
470 s = (A[n1-1,n1-1] - A[n1-2,n1-2]) ** 2 + 4 * A[n1-1,n1-2] * A[n1-2,n1-1]
471 if ctx.re(s) > 0:
472 s = ctx.sqrt(s)
473 else:
474 s = ctx.sqrt(-s) * 1j
475 a = (t + s) / 2
476 b = (t - s) / 2
477 if abs(A[n1-1,n1-1] - a) > abs(A[n1-1,n1-1] - b):
478 shift = b
479 else:
480 shift = a
482 its += 1
483 totalits += 1
485 qr_step(ctx, n0, n1, A, Q, shift)
487 if its > maxits:
488 raise RuntimeError("qr: failed to converge after %d steps" % its)
491@defun
492def schur(ctx, A, overwrite_a = False):
493 """
494 This routine computes the Schur decomposition of a square matrix A.
495 Given A, an unitary matrix Q is determined such that
497 Q' A Q = R and Q' Q = Q Q' = 1
499 where R is an upper right triangular matrix. Here ' denotes the
500 hermitian transpose (i.e. transposition and conjugation).
502 input:
503 A : a real or complex square matrix
504 overwrite_a : if true, allows modification of A which may improve
505 performance. if false, A is not modified.
507 output:
508 Q : an unitary matrix
509 R : an upper right triangular matrix
511 return value: (Q, R)
513 example:
514 >>> from mpmath import mp
515 >>> A = mp.matrix([[3, -1, 2], [2, 5, -5], [-2, -3, 7]])
516 >>> Q, R = mp.schur(A)
517 >>> mp.nprint(R, 3) # doctest:+SKIP
518 [2.0 0.417 -2.53]
519 [0.0 4.0 -4.74]
520 [0.0 0.0 9.0]
521 >>> print(mp.chop(A - Q * R * Q.transpose_conj()))
522 [0.0 0.0 0.0]
523 [0.0 0.0 0.0]
524 [0.0 0.0 0.0]
526 warning: The Schur decomposition is not unique.
527 """
529 n = A.rows
531 if n == 1:
532 return (ctx.matrix([[1]]), A)
534 if not overwrite_a:
535 A = A.copy()
537 T = ctx.matrix(n, 1)
539 hessenberg_reduce_0(ctx, A, T)
540 Q = A.copy()
541 hessenberg_reduce_1(ctx, Q, T)
543 for x in xrange(n):
544 for y in xrange(x + 2, n):
545 A[y,x] = 0
547 hessenberg_qr(ctx, A, Q)
549 return Q, A
552def eig_tr_r(ctx, A):
553 """
554 This routine calculates the right eigenvectors of an upper right triangular matrix.
556 input:
557 A an upper right triangular matrix
559 output:
560 ER a matrix whose columns form the right eigenvectors of A
562 return value: ER
563 """
565 # this subroutine is inspired by the lapack routines ctrevc.f,clatrs.f
567 n = A.rows
569 ER = ctx.eye(n)
571 eps = ctx.eps
573 unfl = ctx.ldexp(ctx.one, -ctx.prec * 30)
574 # since mpmath effectively has no limits on the exponent, we simply scale doubles up
575 # original double has prec*20
577 smlnum = unfl * (n / eps)
578 simin = 1 / ctx.sqrt(eps)
580 rmax = 1
582 for i in xrange(1, n):
583 s = A[i,i]
585 smin = max(eps * abs(s), smlnum)
587 for j in xrange(i - 1, -1, -1):
589 r = 0
590 for k in xrange(j + 1, i + 1):
591 r += A[j,k] * ER[k,i]
593 t = A[j,j] - s
594 if abs(t) < smin:
595 t = smin
597 r = -r / t
598 ER[j,i] = r
600 rmax = max(rmax, abs(r))
601 if rmax > simin:
602 for k in xrange(j, i+1):
603 ER[k,i] /= rmax
604 rmax = 1
606 if rmax != 1:
607 for k in xrange(0, i + 1):
608 ER[k,i] /= rmax
610 return ER
612def eig_tr_l(ctx, A):
613 """
614 This routine calculates the left eigenvectors of an upper right triangular matrix.
616 input:
617 A an upper right triangular matrix
619 output:
620 EL a matrix whose rows form the left eigenvectors of A
622 return value: EL
623 """
625 n = A.rows
627 EL = ctx.eye(n)
629 eps = ctx.eps
631 unfl = ctx.ldexp(ctx.one, -ctx.prec * 30)
632 # since mpmath effectively has no limits on the exponent, we simply scale doubles up
633 # original double has prec*20
635 smlnum = unfl * (n / eps)
636 simin = 1 / ctx.sqrt(eps)
638 rmax = 1
640 for i in xrange(0, n - 1):
641 s = A[i,i]
643 smin = max(eps * abs(s), smlnum)
645 for j in xrange(i + 1, n):
647 r = 0
648 for k in xrange(i, j):
649 r += EL[i,k] * A[k,j]
651 t = A[j,j] - s
652 if abs(t) < smin:
653 t = smin
655 r = -r / t
656 EL[i,j] = r
658 rmax = max(rmax, abs(r))
659 if rmax > simin:
660 for k in xrange(i, j + 1):
661 EL[i,k] /= rmax
662 rmax = 1
664 if rmax != 1:
665 for k in xrange(i, n):
666 EL[i,k] /= rmax
668 return EL
670@defun
671def eig(ctx, A, left = False, right = True, overwrite_a = False):
672 """
673 This routine computes the eigenvalues and optionally the left and right
674 eigenvectors of a square matrix A. Given A, a vector E and matrices ER
675 and EL are calculated such that
677 A ER[:,i] = E[i] ER[:,i]
678 EL[i,:] A = EL[i,:] E[i]
680 E contains the eigenvalues of A. The columns of ER contain the right eigenvectors
681 of A whereas the rows of EL contain the left eigenvectors.
684 input:
685 A : a real or complex square matrix of shape (n, n)
686 left : if true, the left eigenvectors are calculated.
687 right : if true, the right eigenvectors are calculated.
688 overwrite_a : if true, allows modification of A which may improve
689 performance. if false, A is not modified.
691 output:
692 E : a list of length n containing the eigenvalues of A.
693 ER : a matrix whose columns contain the right eigenvectors of A.
694 EL : a matrix whose rows contain the left eigenvectors of A.
696 return values:
697 E if left and right are both false.
698 (E, ER) if right is true and left is false.
699 (E, EL) if left is true and right is false.
700 (E, EL, ER) if left and right are true.
703 examples:
704 >>> from mpmath import mp
705 >>> A = mp.matrix([[3, -1, 2], [2, 5, -5], [-2, -3, 7]])
706 >>> E, ER = mp.eig(A)
707 >>> print(mp.chop(A * ER[:,0] - E[0] * ER[:,0]))
708 [0.0]
709 [0.0]
710 [0.0]
712 >>> E, EL, ER = mp.eig(A,left = True, right = True)
713 >>> E, EL, ER = mp.eig_sort(E, EL, ER)
714 >>> mp.nprint(E)
715 [2.0, 4.0, 9.0]
716 >>> print(mp.chop(A * ER[:,0] - E[0] * ER[:,0]))
717 [0.0]
718 [0.0]
719 [0.0]
720 >>> print(mp.chop( EL[0,:] * A - EL[0,:] * E[0]))
721 [0.0 0.0 0.0]
723 warning:
724 - If there are multiple eigenvalues, the eigenvectors do not necessarily
725 span the whole vectorspace, i.e. ER and EL may have not full rank.
726 Furthermore in that case the eigenvectors are numerical ill-conditioned.
727 - In the general case the eigenvalues have no natural order.
729 see also:
730 - eigh (or eigsy, eighe) for the symmetric eigenvalue problem.
731 - eig_sort for sorting of eigenvalues and eigenvectors
732 """
734 n = A.rows
736 if n == 1:
737 if left and (not right):
738 return ([A[0]], ctx.matrix([[1]]))
740 if right and (not left):
741 return ([A[0]], ctx.matrix([[1]]))
743 return ([A[0]], ctx.matrix([[1]]), ctx.matrix([[1]]))
745 if not overwrite_a:
746 A = A.copy()
748 T = ctx.zeros(n, 1)
750 hessenberg_reduce_0(ctx, A, T)
752 if left or right:
753 Q = A.copy()
754 hessenberg_reduce_1(ctx, Q, T)
755 else:
756 Q = False
758 for x in xrange(n):
759 for y in xrange(x + 2, n):
760 A[y,x] = 0
762 hessenberg_qr(ctx, A, Q)
764 E = [0 for i in xrange(n)]
765 for i in xrange(n):
766 E[i] = A[i,i]
768 if not (left or right):
769 return E
771 if left:
772 EL = eig_tr_l(ctx, A)
773 EL = EL * Q.transpose_conj()
775 if right:
776 ER = eig_tr_r(ctx, A)
777 ER = Q * ER
779 if left and (not right):
780 return (E, EL)
782 if right and (not left):
783 return (E, ER)
785 return (E, EL, ER)
787@defun
788def eig_sort(ctx, E, EL = False, ER = False, f = "real"):
789 """
790 This routine sorts the eigenvalues and eigenvectors delivered by ``eig``.
792 parameters:
793 E : the eigenvalues as delivered by eig
794 EL : the left eigenvectors as delivered by eig, or false
795 ER : the right eigenvectors as delivered by eig, or false
796 f : either a string ("real" sort by increasing real part, "imag" sort by
797 increasing imag part, "abs" sort by absolute value) or a function
798 mapping complexs to the reals, i.e. ``f = lambda x: -mp.re(x) ``
799 would sort the eigenvalues by decreasing real part.
801 return values:
802 E if EL and ER are both false.
803 (E, ER) if ER is not false and left is false.
804 (E, EL) if EL is not false and right is false.
805 (E, EL, ER) if EL and ER are not false.
807 example:
808 >>> from mpmath import mp
809 >>> A = mp.matrix([[3, -1, 2], [2, 5, -5], [-2, -3, 7]])
810 >>> E, EL, ER = mp.eig(A,left = True, right = True)
811 >>> E, EL, ER = mp.eig_sort(E, EL, ER)
812 >>> mp.nprint(E)
813 [2.0, 4.0, 9.0]
814 >>> E, EL, ER = mp.eig_sort(E, EL, ER,f = lambda x: -mp.re(x))
815 >>> mp.nprint(E)
816 [9.0, 4.0, 2.0]
817 >>> print(mp.chop(A * ER[:,0] - E[0] * ER[:,0]))
818 [0.0]
819 [0.0]
820 [0.0]
821 >>> print(mp.chop( EL[0,:] * A - EL[0,:] * E[0]))
822 [0.0 0.0 0.0]
823 """
825 if isinstance(f, str):
826 if f == "real":
827 f = ctx.re
828 elif f == "imag":
829 f = ctx.im
830 elif f == "abs":
831 f = abs
832 else:
833 raise RuntimeError("unknown function %s" % f)
835 n = len(E)
837 # Sort eigenvalues (bubble-sort)
839 for i in xrange(n):
840 imax = i
841 s = f(E[i]) # s is the current maximal element
843 for j in xrange(i + 1, n):
844 c = f(E[j])
845 if c < s:
846 s = c
847 imax = j
849 if imax != i:
850 # swap eigenvalues
852 z = E[i]
853 E[i] = E[imax]
854 E[imax] = z
856 if not isinstance(EL, bool):
857 for j in xrange(n):
858 z = EL[i,j]
859 EL[i,j] = EL[imax,j]
860 EL[imax,j] = z
862 if not isinstance(ER, bool):
863 for j in xrange(n):
864 z = ER[j,i]
865 ER[j,i] = ER[j,imax]
866 ER[j,imax] = z
868 if isinstance(EL, bool) and isinstance(ER, bool):
869 return E
871 if isinstance(EL, bool) and not(isinstance(ER, bool)):
872 return (E, ER)
874 if isinstance(ER, bool) and not(isinstance(EL, bool)):
875 return (E, EL)
877 return (E, EL, ER)