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__ = ['lu_factor', 'lu_solve', 'rcond_from_lu'] 

10 

11import numpy as np 

12from . import lapack 

13 

14 

15def lu_factor(a, overwrite_a=False): 

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

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

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

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

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

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

22 

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

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

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

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

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

28 is standard, this version is better to use. 

29 

30 Parameters 

31 ---------- 

32 a : array, shape (M, M) 

33 Matrix to factorize 

34 overwrite_a : boolean 

35 Whether to overwrite data in a (may increase performance) 

36 

37 Returns 

38 ------- 

39 lu : array, shape (N, N) 

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

41 The unit diagonal elements of L are not stored. 

42 piv : array, shape (N,) 

43 Pivot indices representing the permutation matrix P: 

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

45 singular : boolean 

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

47 """ 

48 a = lapack.prepare_for_lapack(overwrite_a, a) 

49 return lapack.getrf(a) 

50 

51 

52def lu_solve(matrix_factorization, b): 

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

54 factorization of a 

55 

56 Parameters 

57 ---------- 

58 matrix_factorization 

59 Factorization of the coefficient matrix a, as given by lu_factor 

60 b : array (vector or matrix) 

61 Right-hand side 

62 

63 Returns 

64 ------- 

65 x : array (vector or matrix) 

66 Solution to the system 

67 """ 

68 (lu, ipiv, singular) = matrix_factorization 

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

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

71 "a singular matrix. Result of solve step " 

72 "are probably unreliable") 

73 

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

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

76 return lapack.getrs(lu, ipiv, b) 

77 

78 

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

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

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

82 norm_a. 

83 

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

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

86 

87 Parameters 

88 ---------- 

89 matrix_factorization 

90 Factorization of the matrix a, as given by lu_factor 

91 norm_a : float or complex 

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

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

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

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

96 

97 Returns 

98 ------- 

99 rcond : float or complex 

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

101 norm specified in norm 

102 """ 

103 (lu, ipiv, singular) = matrix_factorization 

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

105 lu = lapack.prepare_for_lapack(False, lu) 

106 return lapack.gecon(lu, norm_a, norm)