decomp_lu.py 3.93 KB
Newer Older
1
# Copyright 2011-2013 Kwant authors.
2
#
Christoph Groth's avatar
Christoph Groth committed
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
5
# https://kwant-project.org/license.  A list of Kwant authors can be found in
Christoph Groth's avatar
Christoph Groth committed
6
# the file AUTHORS.rst at the top-level directory of this distribution and at
7
# https://kwant-project.org/authors.
8

9
10
11
12
13
__all__ = ['lu_factor', 'lu_solve', 'rcond_from_lu']

import numpy as np
from . import lapack

Anton Akhmerov's avatar
Anton Akhmerov committed
14
15

def lu_factor(a, overwrite_a=False):
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
25
    lu_factor in SciPy has a bug that depending on the type of NumPy matrix
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)
    """
48
    a = lapack.prepare_for_lapack(overwrite_a, a)
Joseph Weston's avatar
Joseph Weston committed
49
    return lapack.getrf(a)
50

Anton Akhmerov's avatar
Anton Akhmerov committed
51

52
def lu_solve(matrix_factorization, b):
53
54
55
56
57
    """Solve a linear system of equations, a x = b, given the LU
    factorization of a

    Parameters
    ----------
58
    matrix_factorization
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
    """
68
    (lu, ipiv, singular) = matrix_factorization
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")

74
    lu, b = lapack.prepare_for_lapack(False, lu, b)
Anton Akhmerov's avatar
Anton Akhmerov committed
75
    ipiv = np.ascontiguousarray(np.asanyarray(ipiv), dtype=lapack.int_dtype)
Joseph Weston's avatar
Joseph Weston committed
76
    return lapack.getrs(lu, ipiv, b)
77

Anton Akhmerov's avatar
Anton Akhmerov committed
78

79
def rcond_from_lu(matrix_factorization, norm_a, norm="1"):
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
    ----------
89
    matrix_factorization
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
    """
103
    (lu, ipiv, singular) = matrix_factorization
104
    norm = norm.encode('utf8')  # lapack expects bytes
105
    lu = lapack.prepare_for_lapack(False, lu)
Joseph Weston's avatar
Joseph Weston committed
106
    return lapack.gecon(lu, norm_a, norm)