1# Copyright 2011-2013 Kwant authors. 

2# 

3# This file is part of Kwant. It is subject to the license terms in the file 

4# LICENSE.rst found in the top-level directory of this distribution and at 

5# https://kwant-project.org/license. A list of Kwant authors can be found in 

6# the file AUTHORS.rst at the top-level directory of this distribution and at 

7# https://kwant-project.org/authors. 

8 

9__all__ = ['smatrix', 'ldos', 'wave_function', 'greens_function', 'options', 

10 'Solver'] 

11 

12import numpy as np 

13from . import common 

14from ..linalg import mumps 

15 

16 

17class Solver(common.SparseSolver): 

18 """Sparse Solver class based on the sparse direct solver MUMPS.""" 

19 

20 lhsformat = 'coo' 

21 rhsformat = 'csc' 

22 

23 def __init__(self): 

24 self.nrhs = self.ordering = self.sparse_rhs = None 

25 self.reset_options() 

26 

27 def reset_options(self): 

28 """Set the options to default values. Return the old options.""" 

29 return self.options(nrhs=6, ordering='kwant_decides', sparse_rhs=False) 

30 

31 def options(self, nrhs=None, ordering=None, sparse_rhs=None): 

32 """ 

33 Modify some options. Return the old options. 

34 

35 Parameters 

36 ---------- 

37 nrhs : number 

38 number of right hand sides that should be solved simultaneously. A 

39 value around 5-10 gives optimal performance on many machines. If 

40 memory is an issue, it can be set to 1, to minimize memory usage 

41 (at the cost of slower performance). Default value is 6. 

42 ordering : string 

43 one of the ordering methods supported by the MUMPS solver (see 

44 ``kwant.linalg.mumps``. The availability of certain orderings 

45 depends on the MUMPS installation.), or 'kwant_decides'. If 

46 ``ordering=='kwant_decides'``, the ordering that typically gives 

47 the best performance is chosen from the available ones. One can 

48 also defer the choice of ordering to MUMPS by specifying 'auto', in 

49 some cases MUMPS however chooses poorly. 

50 

51 The choice of ordering can significantly influence the performance 

52 and memory impact of the solve phase. Typically the nested 

53 dissection orderings 'metis' and 'scotch' are most suited for 

54 physical systems. Default is 'kwant_decides' 

55 sparse_rhs : True or False 

56 whether to use a sparse right hand side in the solve phase of 

57 MUMPS. Preliminary tests have not shown a significant performance 

58 increase when this feature is used, but this needs more looking 

59 into. Default value is False. 

60 

61 Returns 

62 ------- 

63 old_options: dict 

64 dictionary containing the previous options. 

65 

66 Notes 

67 ----- 

68 Thanks to this method returning the old options as a dictionary it is 

69 easy to change some options temporarily: 

70 

71 >>> saved_options = kwant.solvers.mumps.options(nrhs=12) 

72 >>> some_code() 

73 >>> kwant.solvers.mumps.options(**saved_options) 

74 """ 

75 

76 old_opts = {'nrhs': self.nrhs, 

77 'ordering': self.ordering, 

78 'sparse_rhs': self.sparse_rhs} 

79 

80 if nrhs is not None: 

81 if nrhs < 1 and int(nrhs) != nrhs: 81 ↛ 82line 81 didn't jump to line 82, because the condition on line 81 was never true

82 raise ValueError("nrhs must be an integer bigger than zero") 

83 nrhs = int(nrhs) 

84 self.nrhs = nrhs 

85 

86 if ordering is not None: 

87 if ordering == 'kwant_decides': 

88 # Choose what is considered to be the best ordering. 

89 sorted_orderings = [order 

90 for order in ['metis', 'scotch', 'auto'] 

91 if order in mumps.possible_orderings()] 

92 ordering = sorted_orderings[0] 

93 elif ordering not in mumps.orderings: 93 ↛ 94line 93 didn't jump to line 94, because the condition on line 93 was never true

94 raise ValueError("Invalid ordering: " + ordering) 

95 self.ordering = ordering 

96 

97 if sparse_rhs is not None: 

98 self.sparse_rhs = bool(sparse_rhs) 

99 

100 return old_opts 

101 

102 def _factorized(self, a): 

103 inst = mumps.MUMPSContext() 

104 inst.factor(a, ordering=self.ordering) 

105 return inst 

106 

107 def _solve_linear_sys(self, factorized_a, b, kept_vars): 

108 if b.shape[1] == 0: 

109 return b[kept_vars] 

110 

111 solve = factorized_a.solve 

112 sols = [] 

113 

114 for j in range(0, b.shape[1], self.nrhs): 

115 tmprhs = b[:, j:min(j + self.nrhs, b.shape[1])] 

116 

117 if not self.sparse_rhs: 

118 tmprhs = tmprhs.toarray() 

119 sols.append(solve(tmprhs)[kept_vars, :]) 

120 

121 return np.concatenate(sols, axis=1) 

122 

123 

124default_solver = Solver() 

125 

126smatrix = default_solver.smatrix 

127greens_function = default_solver.greens_function 

128ldos = default_solver.ldos 

129wave_function = default_solver.wave_function 

130options = default_solver.options 

131reset_options = default_solver.reset_options