Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

# Copyright 2011-2013 Kwant authors. 

# 

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

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

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

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

# http://kwant-project.org/authors. 

 

__all__ = ['lu_factor', 'lu_solve', 'rcond_from_lu'] 

 

import numpy as np 

from . import lapack 

 

 

def lu_factor(a, overwrite_a=False): 

"""Compute the LU factorization of a matrix A = P * L * U. The function 

returns a tuple (lu, p, singular), where lu contains the LU factorization 

storing the unit lower triangular matrix L in the strictly lower triangle 

(the unit diagonal is not stored) and the upper triangular matrix U in the 

upper triangle. p is a vector of pivot indices, and singular a Boolean 

value indicating whether the matrix A is singular up to machine precision. 

 

NOTE: This function mimics the behavior of scipy.linalg.lu_factor (except 

that it has in addition the flag singular). The main reason is that 

lu_factor in SciPy has a bug that depending on the type of NumPy matrix 

passed to it, it would not return what was descirbed in the 

documentation. This bug will be (probably) fixed in 0.10.0 but until this 

is standard, this version is better to use. 

 

Parameters 

---------- 

a : array, shape (M, M) 

Matrix to factorize 

overwrite_a : boolean 

Whether to overwrite data in a (may increase performance) 

 

Returns 

------- 

lu : array, shape (N, N) 

Matrix containing U in its upper triangle, and L in its lower triangle. 

The unit diagonal elements of L are not stored. 

piv : array, shape (N,) 

Pivot indices representing the permutation matrix P: 

row i of matrix was interchanged with row piv[i]. 

singular : boolean 

Whether the matrix a is singular (up to machine precision) 

""" 

 

ltype, a = lapack.prepare_for_lapack(overwrite_a, a) 

 

51 ↛ 52line 51 didn't jump to line 52, because the condition on line 51 was never true if a.ndim != 2: 

raise ValueError("lu_factor expects a matrix") 

 

if ltype == 'd': 

return lapack.dgetrf(a) 

elif ltype == 'z': 

return lapack.zgetrf(a) 

elif ltype == 's': 

return lapack.sgetrf(a) 

else: 

return lapack.cgetrf(a) 

 

 

def lu_solve(matrix_factorization, b): 

"""Solve a linear system of equations, a x = b, given the LU 

factorization of a 

 

Parameters 

---------- 

matrix_factorization 

Factorization of the coefficient matrix a, as given by lu_factor 

b : array (vector or matrix) 

Right-hand side 

 

Returns 

------- 

x : array (vector or matrix) 

Solution to the system 

""" 

(lu, ipiv, singular) = matrix_factorization 

81 ↛ 82line 81 didn't jump to line 82, because the condition on line 81 was never true if singular: 

raise RuntimeWarning("In lu_solve: the flag singular indicates " 

"a singular matrix. Result of solve step " 

"are probably unreliable") 

 

ltype, lu, b = lapack.prepare_for_lapack(False, lu, b) 

ipiv = np.ascontiguousarray(np.asanyarray(ipiv), dtype=lapack.int_dtype) 

 

89 ↛ 90line 89 didn't jump to line 90, because the condition on line 89 was never true if b.ndim > 2: 

raise ValueError("lu_solve: b must be a vector or matrix") 

 

92 ↛ 93line 92 didn't jump to line 93, because the condition on line 92 was never true if lu.shape[0] != b.shape[0]: 

raise ValueError("lu_solve: incompatible dimensions of b") 

 

if ltype == 'd': 

return lapack.dgetrs(lu, ipiv, b) 

elif ltype == 'z': 

return lapack.zgetrs(lu, ipiv, b) 

elif ltype == 's': 

return lapack.sgetrs(lu, ipiv, b) 

else: 

return lapack.cgetrs(lu, ipiv, b) 

 

 

def rcond_from_lu(matrix_factorization, norm_a, norm="1"): 

"""Compute the reciprocal condition number from the LU decomposition as 

returned from lu_factor(), given additionally the norm of the matrix a in 

norm_a. 

 

The reciprocal condition number is given as 1/(||A||*||A^-1||), where 

||...|| is a matrix norm. 

 

Parameters 

---------- 

matrix_factorization 

Factorization of the matrix a, as given by lu_factor 

norm_a : float or complex 

norm of the original matrix a (type of norm is specified in norm) 

norm : {'1', 'I'}, optional 

type of matrix norm which should be used to compute the condition 

number ("1": 1-norm, "I": infinity norm). Default: '1'. 

 

Returns 

------- 

rcond : float or complex 

reciprocal condition number of a with respect to the type of matrix 

norm specified in norm 

""" 

(lu, ipiv, singular) = matrix_factorization 

130 ↛ 131line 130 didn't jump to line 131, because the condition on line 130 was never true if not norm in ("1", "I"): 

raise ValueError("norm in rcond_from_lu must be either '1' or 'I'") 

norm = norm.encode('utf8') # lapack expects bytes 

 

ltype, lu = lapack.prepare_for_lapack(False, lu) 

 

if ltype == 'd': 

return lapack.dgecon(lu, norm_a, norm) 

elif ltype == 'z': 

return lapack.zgecon(lu, norm_a, norm) 

elif ltype == 's': 

return lapack.sgecon(lu, norm_a, norm) 

else: 

return lapack.cgecon(lu, norm_a, norm)