decomp_lu.py 3.93 KB
 Christoph Groth committed Jul 31, 2013 1 ``````# Copyright 2011-2013 Kwant authors. `````` Christoph Groth committed Feb 06, 2013 2 ``````# `````` Christoph Groth committed Oct 20, 2015 3 4 ``````# 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 `````` Anton Akhmerov committed Dec 22, 2020 5 ``````# https://kwant-project.org/license. A list of Kwant authors can be found in `````` Christoph Groth committed Oct 20, 2015 6 ``````# the file AUTHORS.rst at the top-level directory of this distribution and at `````` Anton Akhmerov committed Dec 22, 2020 7 ``````# https://kwant-project.org/authors. `````` Christoph Groth committed Feb 06, 2013 8 `````` `````` Kwant authors committed Jan 10, 2011 9 10 11 12 13 ``````__all__ = ['lu_factor', 'lu_solve', 'rcond_from_lu'] import numpy as np from . import lapack `````` Anton Akhmerov committed Aug 19, 2012 14 15 `````` def lu_factor(a, overwrite_a=False): `````` Kwant authors committed Jan 10, 2011 16 17 18 19 20 21 22 23 24 `````` """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 `````` Christoph Groth committed Nov 21, 2012 25 `````` lu_factor in SciPy has a bug that depending on the type of NumPy matrix `````` Kwant authors committed Jan 10, 2011 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 `````` 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) """ `````` Joseph Weston committed Sep 13, 2017 48 `````` a = lapack.prepare_for_lapack(overwrite_a, a) `````` Joseph Weston committed Sep 13, 2017 49 `````` return lapack.getrf(a) `````` Kwant authors committed Jan 10, 2011 50 `````` `````` Anton Akhmerov committed Aug 19, 2012 51 `````` `````` Joseph Weston committed Nov 07, 2015 52 ``````def lu_solve(matrix_factorization, b): `````` Kwant authors committed Jan 10, 2011 53 54 55 56 57 `````` """Solve a linear system of equations, a x = b, given the LU factorization of a Parameters ---------- `````` Joseph Weston committed Nov 07, 2015 58 `````` matrix_factorization `````` Kwant authors committed Jan 10, 2011 59 60 61 62 63 64 65 66 67 `````` 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 """ `````` Joseph Weston committed Nov 07, 2015 68 `````` (lu, ipiv, singular) = matrix_factorization `````` Kwant authors committed Jan 10, 2011 69 70 71 72 73 `````` if singular: raise RuntimeWarning("In lu_solve: the flag singular indicates " "a singular matrix. Result of solve step " "are probably unreliable") `````` Joseph Weston committed Sep 13, 2017 74 `````` lu, b = lapack.prepare_for_lapack(False, lu, b) `````` Anton Akhmerov committed Aug 19, 2012 75 `````` ipiv = np.ascontiguousarray(np.asanyarray(ipiv), dtype=lapack.int_dtype) `````` Joseph Weston committed Sep 13, 2017 76 `````` return lapack.getrs(lu, ipiv, b) `````` Kwant authors committed Jan 10, 2011 77 `````` `````` Anton Akhmerov committed Aug 19, 2012 78 `````` `````` Joseph Weston committed Nov 07, 2015 79 ``````def rcond_from_lu(matrix_factorization, norm_a, norm="1"): `````` Kwant authors committed Jan 10, 2011 80 81 82 83 84 85 86 87 88 `````` """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 ---------- `````` Joseph Weston committed Nov 07, 2015 89 `````` matrix_factorization `````` Kwant authors committed Jan 10, 2011 90 91 92 93 94 95 96 97 98 99 100 101 102 `````` 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 """ `````` Joseph Weston committed Nov 07, 2015 103 `````` (lu, ipiv, singular) = matrix_factorization `````` Joseph Weston committed Nov 05, 2015 104 `````` norm = norm.encode('utf8') # lapack expects bytes `````` Joseph Weston committed Sep 13, 2017 105 `````` lu = lapack.prepare_for_lapack(False, lu) `````` Joseph Weston committed Sep 13, 2017 106 `` return lapack.gecon(lu, norm_a, norm)``