Source Code for Module linear_operators.iterative.utils

 1  """ 
 2  Useful functions for LinearOperators 
 3  """ 
 4  import numpy as np 
 5   
 6  try: 
 7      from scipy.sparse.linalg import arpack 
 8  except(ImportError): 
 9      import arpack 
10   
11  # to handle new names of eigen and eigen_symmetric 
12  if 'eigen' in arpack.__dict__.keys(): 
13      eigs = arpack.eigen 
14      eigsh = arpack.eigen_symmetric 
15  elif 'eigs' in arpack.__dict__.keys(): 
16      eigs = arpack.eigs 
17      eigsh = arpack.eigsh 
18  else: 
19      raise ImportError("Unable to find eigen or eigs in arpack.") 
20   
21  from ..operators import SymmetricOperator 
22   
23 -def eigendecomposition(A, **kwargs): 
24      """ 
25      A wrapper to arpack eigen_symmetric which output an approximation 
26      of the A input matrix, as a LinearOperator storing eigenvalues and 
27      eigenvectors. 
28   
29      It passes its arguments arguments as arpack.eigsh but 
30      forces return_eigenvectors to True. 
31      """ 
32      from ..operators import EigendecompositionOperator 
33      return EigendecompositionOperator(A, **kwargs) 
34   
35 -def cond(A, k=2, kl=None, ks=None, symmetric=True, M=None, maxiter=None, 
36           tol=1e-6, verbose=False, prune_zeros=False, loweig=None): 
37      """ 
38      Find the condition number of a LinearOperator using arpack eigen 
39      function. 
40      """ 
41      if kl is None: 
42          kl = k 
43      if ks is None: 
44          ks = k 
45   
46      if symmetric: 
47          eigen = eigsh 
48      else: 
49          eigen = eigs 
50   
51      vmax = eigen(A, which='LM', k=kl, M=M, maxiter=maxiter, tol=tol, 
52                   return_eigenvectors=False) 
53      vmin = eigen(A, which='SM', k=ks, M=M, maxiter=maxiter, tol=tol, 
54                   return_eigenvectors=False) 
55      nb_null = np.sum(vmin == 0) 
56      if verbose: 
57          print("Found %d zero-valued eigenvalues" % nb_null) 
58          if nb_null == len(vmin): 
59              print("All small eigenvalues were zeros.") 
60      if prune_zeros: 
61          vmin = vmin[vmin !=0] 
62      if vmin.size == 0: 
63          return np.inf 
64      smin = np.min(vmin) 
65      smax = np.max(vmax) 
66      if smin == 0: 
67          return np.inf 
68      if loweig is not None: 
69          if smin < loweig: 
70              return np.inf 
71      if verbose: 
72          print vmin, vmax 
73          print smin, smax 
74   
75      return smax / smin 
76