Skip to content
GitLab
Verified Commit 5832d028 authored by Anton Akhmerov's avatar Anton Akhmerov
Browse files

remove custom generalized eigenproblem wrapper

parent 4a0a7962
Pipeline #74397 passed with stages
in 10 minutes and 28 seconds
Hide whitespace changes
Inline Side-by-side
kwant/linalg/__init__.py
View file @ 5832d028
# Copyright 2011-2013 Kwant authors.
# Copyright 2011-2021 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
......@@ -11,7 +11,5 @@ from . import lapack
# Merge the public interface of the other submodules.
from .decomp_schur import *
from .decomp_ev import *
__all__.extend([decomp_ev.__all__,
decomp_schur.__all__])
__all__.extend([decomp_schur.__all__])
kwant/linalg/decomp_ev.py deleted 100644 → 0
View file @ 4a0a7962
# 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
# https://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
# https://kwant-project.org/authors.
__all__ = ['gen_eig']
from . import lapack
def gen_eig(a, b, left=False, right=True, overwrite_ab=False):
"""Compute the eigenvalues and -vectors of the matrix pencil (a,b), i.e. of
the generalized (unsymmetric) eigenproblem a v = lambda b v where a and b
are square (unsymmetric) matrices, v the eigenvector and lambda the
eigenvalues.
The eigenvalues are returned as numerator alpha and denominator beta,
i.e. lambda = alpha/beta. This is advantageous, as lambda can be infinity
which is well-defined in this case as beta = 0.
Parameters
----------
a : array, shape (M, M)
b : array, shape (M, M)
`a` and `b` are the two matrices defining the generalized eigenproblem
left : boolean
Whether to calculate and return left eigenvectors
right : boolean
Whether to calculate and return right eigenvectors
overwrite_ab : boolean
Whether to overwrite data in `a` and `b` (may improve performance)
Returns
-------
alpha : complex array, shape (M,)
beta : real or complex array, shape (M,)
The eigenvalues in the form ``alpha/beta``
(if left == True)
vl : double or complex array, shape (M, M)
The left eigenvector corresponding to the eigenvalue
``alpha[i]/beta[i]`` is the column ``vl[:,i]``.
(if right == True)
vr : double or complex array, shape (M, M)
The right eigenvector corresponding to the eigenvalue
``alpha[i]/beta[i]`` is the column ``vr[:,i]``.
"""
a, b = lapack.prepare_for_lapack(overwrite_ab, a, b)
return lapack.ggev(a, b, left, right)
kwant/linalg/decomp_schur.py
View file @ 5832d028
......@@ -302,8 +302,7 @@ def gen_schur(a, b, calc_q=True, calc_z=True, calc_ev=True,
problem (the entries of the diagonal of the complex Schur form are the
eigenvalues of the matrix, for example), and the routine can optionally
also return the generalized eigenvalues in the form (alpha, beta), such
that alpha/beta is a generalized eigenvalue of the pencil (a, b) (see also
gen_eig()).
that alpha/beta is a generalized eigenvalue of the pencil (a, b).
Parameters
----------
......
kwant/linalg/lapack.pyx
View file @ 5832d028
......@@ -183,155 +183,6 @@ def ggev_postprocess(dtype, alphar, alphai, vl_r=None, vr_r=None):
return (alpha, vl, vr)
def ggev(np.ndarray[scalar, ndim=2] A, np.ndarray[scalar, ndim=2] B,
left=False, right=True):
cdef l_int N, info, lwork
# Parameter checks
assert_fortran_mat(A, B)
if A.ndim != 2 or A.ndim != 2:
raise ValueError("gen_eig requires both a and be to be matrices")
if A.shape[0] != A.shape[1]:
raise ValueError("gen_eig requires square matrix input")
if A.shape[0] != B.shape[0] or A.shape[1] != B.shape[1]:
raise ValueError("A and B do not have the same shape")
# Allocate workspaces
N = A.shape[0]
cdef np.ndarray[scalar] alphar, alphai
if scalar in cmplx:
alphar = np.empty(N, dtype=A.dtype)
alphai = None
else:
alphar = np.empty(N, dtype=A.dtype)
alphai = np.empty(N, dtype=A.dtype)
cdef np.ndarray[scalar] beta = np.empty(N, dtype=A.dtype)
cdef np.ndarray rwork = None
if scalar is float_complex:
rwork = np.empty(8 * N, dtype=np.float32)
elif scalar is double_complex:
rwork = np.empty(8 * N, dtype=np.float64)
cdef np.ndarray vl
cdef scalar *vl_ptr
cdef char *jobvl
if left:
vl = np.empty((N,N), dtype=A.dtype, order='F')
vl_ptr = <scalar *>vl.data
jobvl = "V"
else:
vl = None
vl_ptr = NULL
jobvl = "N"
cdef np.ndarray vr
cdef scalar *vr_ptr
cdef char *jobvr
if right:
vr = np.empty((N,N), dtype=A.dtype, order='F')
vr_ptr = <scalar *>vr.data
jobvr = "V"
else:
vr = None
vr_ptr = NULL
jobvr = "N"
# Workspace query
# Xggev expects &qwork as a <scalar *> (even though it's an integer)
lwork = -1
cdef scalar qwork
if scalar is float:
lapack.sggev(jobvl, jobvr, &N, <float *>A.data, &N,
<float *>B.data, &N,
<float *>alphar.data, <float *> alphai.data,
<float *>beta.data,
vl_ptr, &N, vr_ptr, &N,
&qwork, &lwork, &info)
elif scalar is double:
lapack.dggev(jobvl, jobvr, &N, <double *>A.data, &N,
<double *>B.data, &N,
<double *>alphar.data, <double *> alphai.data,
<double *>beta.data,
vl_ptr, &N, vr_ptr, &N,
&qwork, &lwork, &info)
elif scalar is float_complex:
lapack.cggev(jobvl, jobvr, &N, <float complex *>A.data, &N,
<float complex *>B.data, &N,
<float complex *>alphar.data, <float complex *>beta.data,
vl_ptr, &N, vr_ptr, &N,
&qwork, &lwork,
<float *>rwork.data, &info)
elif scalar is double_complex:
lapack.zggev(jobvl, jobvr, &N, <double complex *>A.data, &N,
<double complex *>B.data, &N,
<double complex *>alphar.data, <double complex *>beta.data,
vl_ptr, &N, vr_ptr, &N,
&qwork, &lwork,
<double *>rwork.data, &info)
assert info == 0, "Argument error in ggev"
lwork = lwork_from_qwork(qwork)
cdef np.ndarray[scalar] work = np.empty(lwork, dtype=A.dtype)
# The actual calculation
if scalar is float:
lapack.sggev(jobvl, jobvr, &N, <float *>A.data, &N,
<float *>B.data, &N,
<float *>alphar.data, <float *> alphai.data,
<float *>beta.data,
vl_ptr, &N, vr_ptr, &N,
<float *>work.data, &lwork, &info)
elif scalar is double:
lapack.dggev(jobvl, jobvr, &N, <double *>A.data, &N,
<double *>B.data, &N,
<double *>alphar.data, <double *> alphai.data,
<double *>beta.data,
vl_ptr, &N, vr_ptr, &N,
<double *>work.data, &lwork, &info)
elif scalar is float_complex:
lapack.cggev(jobvl, jobvr, &N, <float complex *>A.data, &N,
<float complex *>B.data, &N,
<float complex *>alphar.data, <float complex *>beta.data,
vl_ptr, &N, vr_ptr, &N,
<float complex *>work.data, &lwork,
<float *>rwork.data, &info)
elif scalar is double_complex:
lapack.zggev(jobvl, jobvr, &N, <double complex *>A.data, &N,
<double complex *>B.data, &N,
<double complex *>alphar.data, <double complex *>beta.data,
vl_ptr, &N, vr_ptr, &N,
<double complex *>work.data, &lwork,
<double *>rwork.data, &info)
if info > 0:
raise LinAlgError("QZ iteration failed to converge in sggev")
assert info == 0, "Argument error in ggev"
if scalar is float:
post_dtype = np.complex64
elif scalar is double:
post_dtype = np.complex128
cdef np.ndarray alpha
alpha = alphar
if scalar in floating:
alpha, vl, vr = ggev_postprocess(post_dtype, alphar, alphai, vl, vr)
return filter_args((True, True, left, right), (alpha, beta, vl, vr))
def gees(np.ndarray[scalar, ndim=2] A, calc_q=True, calc_ev=True):
cdef l_int N, lwork, sdim, info
......
kwant/linalg/tests/test_linalg.py
View file @ 5832d028
......@@ -8,12 +8,11 @@
import pytest
import numpy as np
from scipy import linalg
from kwant.linalg import (
gen_eig, schur,
convert_r2c_schur, order_schur, evecs_from_schur, gen_schur,
convert_r2c_gen_schur, order_gen_schur, evecs_from_gen_schur)
schur, convert_r2c_schur, order_schur, evecs_from_schur, gen_schur,
convert_r2c_gen_schur, order_gen_schur, evecs_from_gen_schur
)
from ._test_utils import _Random, assert_array_almost_equal
......@@ -25,18 +24,6 @@ def dtype(request):
return request.param
def test_gen_eig(dtype):
rand = _Random()
a = rand.randmat(4, 4, dtype)
b = rand.randmat(4, 4, dtype)
(alpha, beta, vl, vr) = gen_eig(a, b, True, True)
assert_array_almost_equal(dtype, a @ vr @ beta, b @ vr @ alpha)
assert_array_almost_equal(dtype, beta @ vl.T.conj() @ a,
alpha @ vl.T.conj() @ b)
def test_schur(dtype):
rand = _Random()
a = rand.randmat(5, 5, dtype)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment