Coverage for kwant/linalg/lapack.pyx : 75%

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
# # 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.
"""Low-level access to LAPACK functions. """
'sgetrs', 'dgetrs', 'cgetrs', 'zgetrs', 'sgecon', 'dgecon', 'cgecon', 'zgecon', 'sggev', 'dggev', 'cggev', 'zggev', 'sgees', 'dgees', 'cgees', 'zgees', 'strsen', 'dtrsen', 'ctrsen', 'ztrsen', 'strevc', 'dtrevc', 'ctrevc', 'ztrevc', 'sgges', 'dgges', 'cgges', 'zgges', 'stgsen', 'dtgsen', 'ctgsen', 'ztgsen', 'stgevc', 'dtgevc', 'ctgevc', 'ztgevc', 'prepare_for_lapack']
cimport numpy as np
# Strangely absolute imports from kwant.linalg.f_lapack don't work here. So # continue using traditional implicit relative imports. from . cimport f_lapack from .f_lapack cimport l_int, l_logical
# exceptions
pass
# some helper functions
# This is a workaround for a bug in NumPy version < 2.0, # where 1x1 matrices do not have the F_Contiguous flag set correctly. raise ValueError("Input matrix must be Fortran contiguous")
# Wrappers for xGETRF cdef l_int M, N, info cdef np.ndarray[l_int, ndim=1] ipiv
<l_int *>ipiv.data, &info)
cdef l_int M, N, info cdef np.ndarray[l_int, ndim=1] ipiv
<l_int *>ipiv.data, &info)
cdef l_int M, N, info cdef np.ndarray[l_int, ndim=1] ipiv
<l_int *>ipiv.data, &info)
cdef l_int M, N, info cdef np.ndarray[l_int, ndim=1] ipiv
<l_int *>ipiv.data, &info)
# Wrappers for xGETRS
np.ndarray[l_int, ndim=1] IPIV, B): cdef l_int N, NRHS, info cdef np.ndarray b
# again: workaround for 1x1-Fortran bug in NumPy < v2.0 raise ValueError("In dgetrs: B must be a Fortran ordered NumPy array")
elif b.ndim == 2: else: raise ValueError("In sgetrs: B must be a vector or matrix")
<l_int *>IPIV.data, <float *>b.data, &N, &info)
np.ndarray[l_int, ndim=1] IPIV, B): cdef l_int N, NRHS, info cdef np.ndarray b
# again: workaround for 1x1-Fortran bug in NumPy < v2.0 raise ValueError("In dgetrs: B must be a Fortran ordered NumPy array")
elif b.ndim == 2: else: raise ValueError("In dgetrs: B must be a vector or matrix")
<l_int *>IPIV.data, <double *>b.data, &N, &info)
np.ndarray[l_int, ndim=1] IPIV, B): cdef l_int N, NRHS, info cdef np.ndarray b
# again: workaround for 1x1-Fortran bug in NumPy < v2.0 raise ValueError("In dgetrs: B must be a Fortran ordered NumPy array")
elif b.ndim == 2: else: raise ValueError("In cgetrs: B must be a vector or matrix")
<l_int *>IPIV.data, <float complex *>b.data, &N, &info)
np.ndarray[l_int, ndim=1] IPIV, B): cdef l_int N, NRHS, info cdef np.ndarray b
# again: workaround for 1x1-Fortran bug in NumPy < v2.0 raise ValueError("In dgetrs: B must be a Fortran ordered NumPy array")
elif b.ndim == 2: else: raise ValueError("In zgetrs: B must be a vector or matrix")
<l_int *>IPIV.data, <double complex *>b.data, &N, &info)
# Wrappers for xGECON
float normA, char *norm = b"1"): cdef l_int N, info cdef float rcond cdef np.ndarray[np.float32_t, ndim=1] work cdef np.ndarray[l_int, ndim=1] iwork
&rcond, <float *>work.data, <l_int *>iwork.data, &info)
double normA, char *norm = b"1"): cdef l_int N, info cdef double rcond cdef np.ndarray[np.float64_t, ndim=1] work cdef np.ndarray[l_int, ndim=1] iwork
&rcond, <double *>work.data, <l_int *>iwork.data, &info)
float normA, char *norm = b"1"): cdef l_int N, info cdef float rcond cdef np.ndarray[np.complex64_t, ndim=1] work cdef np.ndarray[np.float32_t, ndim=1] rwork
&rcond, <float complex *>work.data, <float *>rwork.data, &info)
double normA, char *norm = b"1"): cdef l_int N, info cdef double rcond cdef np.ndarray[np.complex128_t, ndim=1] work cdef np.ndarray[np.float64_t, ndim=1] rwork
&rcond, <double complex *>work.data, <double *>rwork.data, &info)
# Wrappers for xGGEV
# Helper function for xGGEV # depending on whether the eigenvalues are purely real or complex, # some post-processing of the eigenvalues and -vectors is necessary
else: vl = None
else: vr = None else: alpha = alphar vl = vl_r vr = vr_r
np.ndarray[np.float32_t, ndim=2] B, left=False, right=True): cdef l_int N, info, lwork cdef char *jobvl cdef char *jobvr cdef np.ndarray[np.float32_t, ndim=2] vl_r, vr_r cdef float *vl_ptr cdef float *vr_ptr cdef float qwork cdef np.ndarray[np.float32_t, ndim=1] work, alphar, alphai, beta
else: vl_r = None vl_ptr = NULL jobvl = "N"
else: vr_r = None vr_ptr = NULL jobvr = "N"
# workspace query
<float *>B.data, &N, <float *>alphar.data, <float *> alphai.data, <float *>beta.data, vl_ptr, &N, vr_ptr, &N, &qwork, &lwork, &info)
# Now the real calculation <float *>B.data, &N, <float *>alphar.data, <float *> alphai.data, <float *>beta.data, vl_ptr, &N, vr_ptr, &N, <float *>work.data, &lwork, &info)
raise LinAlgError("QZ iteration failed to converge in sggev")
np.ndarray[np.float64_t, ndim=2] B, left=False, right=True): cdef l_int N, info, lwork cdef char *jobvl cdef char *jobvr cdef np.ndarray[np.float64_t, ndim=2] vl_r, vr_r cdef double *vl_ptr cdef double *vr_ptr cdef double qwork cdef np.ndarray[np.float64_t, ndim=1] work, alphar, alphai, beta
else: vl_r = None vl_ptr = NULL jobvl = "N"
else: vr_r = None vr_ptr = NULL jobvr = "N"
# workspace query
<double *>B.data, &N, <double *>alphar.data, <double *> alphai.data, <double *>beta.data, vl_ptr, &N, vr_ptr, &N, &qwork, &lwork, &info)
# Now the real calculation <double *>B.data, &N, <double *>alphar.data, <double *> alphai.data, <double *>beta.data, vl_ptr, &N, vr_ptr, &N, <double *>work.data, &lwork, &info)
raise LinAlgError("QZ iteration failed to converge in dggev")
np.ndarray[np.complex64_t, ndim=2] B, left=False, right=True): cdef l_int N, info, lwork cdef char *jobvl cdef char *jobvr cdef np.ndarray[np.complex64_t, ndim=2] vl, vr cdef float complex *vl_ptr cdef float complex *vr_ptr cdef float complex qwork cdef np.ndarray[np.complex64_t, ndim=1] work, alpha, beta cdef np.ndarray[np.float32_t, ndim=1] rwork
else: vl_ptr = NULL jobvl = "N"
else: vr_ptr = NULL jobvr = "N"
# workspace query
<float complex *>B.data, &N, <float complex *>alpha.data, <float complex *>beta.data, vl_ptr, &N, vr_ptr, &N, &qwork, &lwork, <float *>rwork.data, &info)
# Now the real calculation <float complex *>B.data, &N, <float complex *>alpha.data, <float complex *>beta.data, vl_ptr, &N, vr_ptr, &N, <float complex *>work.data, &lwork, <float *>rwork.data, &info)
raise LinAlgError("QZ iteration failed to converge in cggev")
np.ndarray[np.complex128_t, ndim=2] B, left=False, right=True): cdef l_int N, info, lwork cdef char *jobvl cdef char *jobvr cdef np.ndarray[np.complex128_t, ndim=2] vl, vr cdef double complex *vl_ptr cdef double complex *vr_ptr cdef double complex qwork cdef np.ndarray[np.complex128_t, ndim=1] work, alpha, beta cdef np.ndarray[np.float64_t, ndim=1] rwork
else: vl_ptr = NULL jobvl = "N"
else: vr_ptr = NULL jobvr = "N"
# workspace query
<double complex *>B.data, &N, <double complex *>alpha.data, <double complex *>beta.data, vl_ptr, &N, vr_ptr, &N, &qwork, &lwork, <double *>rwork.data, &info)
# Now the real calculation <double complex *>B.data, &N, <double complex *>alpha.data, <double complex *>beta.data, vl_ptr, &N, vr_ptr, &N, <double complex *>work.data, &lwork, <double *>rwork.data, &info)
raise LinAlgError("QZ iteration failed to converge in zggev")
# Wrapper for xGEES calc_q=True, calc_ev=True): cdef l_int N, lwork, sdim, info cdef char *jobvs cdef float *vs_ptr cdef float qwork cdef np.ndarray[np.float32_t, ndim=2] vs cdef np.ndarray[np.float32_t] wr, wi, work
else: vs_ptr = NULL jobvs = "N"
# workspace query &sdim, <float *>wr.data, <float *>wi.data, vs_ptr, &N, &qwork, &lwork, NULL, &info)
# Now the real calculation &sdim, <float *>wr.data, <float *>wi.data, vs_ptr, &N, <float *>work.data, &lwork, NULL, &info)
raise LinAlgError("QR iteration failed to converge in sgees")
else: w = wr
calc_q=True, calc_ev=True): cdef l_int N, lwork, sdim, info cdef char *jobvs cdef double *vs_ptr cdef double qwork cdef np.ndarray[np.float64_t, ndim=2] vs cdef np.ndarray[np.float64_t] wr, wi, work
else: vs_ptr = NULL jobvs = "N"
# workspace query &sdim, <double *>wr.data, <double *>wi.data, vs_ptr, &N, &qwork, &lwork, NULL, &info)
# Now the real calculation &sdim, <double *>wr.data, <double *>wi.data, vs_ptr, &N, <double *>work.data, &lwork, NULL, &info)
raise LinAlgError("QR iteration failed to converge in dgees")
else:
calc_q=True, calc_ev=True): cdef l_int N, lwork, sdim, info cdef char *jobvs cdef float complex *vs_ptr cdef float complex qwork cdef np.ndarray[np.complex64_t, ndim=2] vs cdef np.ndarray[np.complex64_t] w, work cdef np.ndarray[np.float32_t] rwork
else: vs_ptr = NULL jobvs = "N"
# workspace query &sdim, <float complex *>w.data, vs_ptr, &N, &qwork, &lwork, <float *>rwork.data, NULL, &info)
# Now the real calculation &sdim, <float complex *>w.data, vs_ptr, &N, <float complex *>work.data, &lwork, <float *>rwork.data, NULL, &info)
raise LinAlgError("QR iteration failed to converge in cgees")
calc_q=True, calc_ev=True): cdef l_int N, lwork, sdim, info cdef char *jobvs cdef double complex *vs_ptr cdef double complex qwork cdef np.ndarray[np.complex128_t, ndim=2] vs cdef np.ndarray[np.complex128_t] w, work cdef np.ndarray[np.float64_t] rwork
else: vs_ptr = NULL jobvs = "N"
# workspace query &sdim, <double complex *>w.data, vs_ptr, &N, &qwork, &lwork, <double *>rwork.data, NULL, &info)
# Now the real calculation &sdim, <double complex *>w.data, vs_ptr, &N, <double complex *>work.data, &lwork, <double *>rwork.data, NULL, &info)
raise LinAlgError("QR iteration failed to converge in zgees")
# Wrapper for xTRSEN np.ndarray[np.float32_t, ndim=2] T, np.ndarray[np.float32_t, ndim=2] Q=None, calc_ev=True): cdef l_int N, M, lwork, liwork, qiwork, info cdef char *compq cdef float qwork cdef float *q_ptr cdef np.ndarray[np.float32_t] wr, wi, work cdef np.ndarray[l_int] iwork
assert_fortran_mat(T, Q)
N = T.shape[0] wr = np.empty(N, dtype = np.float32) wi = np.empty(N, dtype = np.float32)
if Q is not None: compq = "V" q_ptr = <float *>Q.data else: compq = "N" q_ptr = NULL
# workspace query lwork = liwork = -1 f_lapack.strsen_("N", compq, <l_logical *>select.data, &N, <float *>T.data, &N, q_ptr, &N, <float *>wr.data, <float *>wi.data, &M, NULL, NULL, &qwork, &lwork, &qiwork, &liwork, &info)
assert info == 0, "Argument error in strsen"
lwork = <int>qwork work = np.empty(lwork, dtype = np.float32) liwork = qiwork iwork = np.empty(liwork, dtype = f_lapack.l_int_dtype)
# Now the real calculation f_lapack.strsen_("N", compq, <l_logical *>select.data, &N, <float *>T.data, &N, q_ptr, &N, <float *>wr.data, <float *>wi.data, &M, NULL, NULL, <float *>work.data, &lwork, <int *>iwork.data, &liwork, &info)
if info > 0: raise LinAlgError("Reordering failed; problem is very ill-conditioned")
assert info == 0, "Argument error in strsen"
if wi.nonzero()[0].size: w = wr + 1j * wi else: w = wr
return filter_args((True, Q is not None, calc_ev), (T, Q, w))
np.ndarray[np.float64_t, ndim=2] T, np.ndarray[np.float64_t, ndim=2] Q=None, calc_ev=True): cdef l_int N, M, lwork, liwork, qiwork, info cdef char *compq cdef double qwork cdef double *q_ptr cdef np.ndarray[np.float64_t] wr, wi, work cdef np.ndarray[l_int] iwork
else: compq = "N" q_ptr = NULL
# workspace query &N, <double *>T.data, &N, q_ptr, &N, <double *>wr.data, <double *>wi.data, &M, NULL, NULL, &qwork, &lwork, &qiwork, &liwork, &info)
# Now the real calculation &N, <double *>T.data, &N, q_ptr, &N, <double *>wr.data, <double *>wi.data, &M, NULL, NULL, <double *>work.data, &lwork, <int *>iwork.data, &liwork, &info)
raise LinAlgError("Reordering failed; problem is very ill-conditioned")
else:
np.ndarray[np.complex64_t, ndim=2] T, np.ndarray[np.complex64_t, ndim=2] Q=None, calc_ev=True): cdef l_int N, M, lwork, info cdef char *compq cdef float complex qwork cdef float complex *q_ptr cdef np.ndarray[np.complex64_t] w, work
else: compq = "N" q_ptr = NULL
# workspace query &N, <float complex *>T.data, &N, q_ptr, &N, <float complex *>w.data, &M, NULL, NULL, &qwork, &lwork, &info)
# Now the real calculation &N, <float complex *>T.data, &N, q_ptr, &N, <float complex *>w.data, &M, NULL, NULL, <float complex *>work.data, &lwork, &info)
raise LinAlgError("Reordering failed; problem is very ill-conditioned")
np.ndarray[np.complex128_t, ndim=2] T, np.ndarray[np.complex128_t, ndim=2] Q=None, calc_ev=True): cdef l_int N, MM, M, lwork, info cdef char *compq cdef double complex qwork cdef double complex *q_ptr cdef np.ndarray[np.complex128_t] w, work
else: compq = "N" q_ptr = NULL
# workspace query &N, <double complex *>T.data, &N, q_ptr, &N, <double complex *>w.data, &M, NULL, NULL, &qwork, &lwork, &info)
# Now the real calculation &N, <double complex *>T.data, &N, q_ptr, &N, <double complex *>w.data, &M, NULL, NULL, <double complex *>work.data, &lwork, &info)
raise LinAlgError("Reordering failed; problem is very ill-conditioned")
# Helper function for xTREVC and xTGEVC cdef int N, M, i, m, indx
# we have the situation of a 2x2 block, and # the eigenvalue with the positive imaginary part desired
# Check if the eigenvalue with negative real part is also # selected, if it is, we need the same entries in vr # we have the situation of a 2x2 block, and # the eigenvalue with the negative imaginary part desired
else: # real eigenvalue
# Wrappers for xTREVC np.ndarray[np.float32_t, ndim=2] Q=None, np.ndarray[l_logical] select=None, left=False, right=True): cdef l_int N, info, M, MM cdef char *side cdef char *howmny cdef np.ndarray[np.float32_t, ndim=2] vl_r, vr_r cdef float *vl_r_ptr cdef float *vr_r_ptr cdef np.ndarray[l_logical] select_cpy cdef l_logical *select_ptr cdef np.ndarray[np.float32_t] work
elif left: side = "L" elif right: side = "R" else: return
# Correct for possible additional storage if a single complex # eigenvalue is selected. # For that: Figure out the positions of the 2x2 blocks.
# Select is overwritten in strevc. order = 'F') else: else: howmny = "A"
else: else: vl_r_ptr = NULL
else: else: vr_r_ptr = NULL
&N, <float *>T.data, &N, vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M, <float *>work.data, &info)
# If there are complex eigenvalues, we need to postprocess the # eigenvectors. else: if left: vl = vl_r if right: vr = vr_r
elif left: return vl else: return vr
np.ndarray[np.float64_t, ndim=2] Q=None, np.ndarray[l_logical] select=None, left=False, right=True): cdef l_int N, info, M, MM cdef char *side cdef char *howmny cdef np.ndarray[np.float64_t, ndim=2] vl_r, vr_r cdef double *vl_r_ptr cdef double *vr_r_ptr cdef np.ndarray[l_logical] select_cpy cdef l_logical *select_ptr cdef np.ndarray[np.float64_t] work
side = "L" else: return
# Correct for possible additional storage if a single complex # eigenvalue is selected. # For that: Figure out the positions of the 2x2 blocks.
# Select is overwritten in dtrevc. order = 'F') else: else: howmny = "A"
else: else:
else: else: vr_r_ptr = NULL
&N, <double *>T.data, &N, vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M, <double *>work.data, &info)
# If there are complex eigenvalues, we need to postprocess the eigenvectors else: vl = vl_r
return vl else:
np.ndarray[np.complex64_t, ndim=2] Q=None, np.ndarray[l_logical] select=None, left=False, right=True): cdef l_int N, info, M, MM cdef char *side cdef char *howmny cdef np.ndarray[np.complex64_t, ndim=2] vl, vr cdef float complex *vl_ptr cdef float complex *vr_ptr cdef l_logical *select_ptr cdef np.ndarray[np.complex64_t] work cdef np.ndarray[np.float32_t] rwork
elif left: side = "L" elif right: side = "R" else: return
else: else: howmny = "A"
else: else: vl_ptr = NULL
else: else: vr_ptr = NULL
&N, <float complex *>T.data, &N, vl_ptr, &N, vr_ptr, &N, &MM, &M, <float complex *>work.data, <float *>rwork.data, &info)
elif left: return vl else: return vr
np.ndarray[np.complex128_t, ndim=2] Q=None, np.ndarray[l_logical] select=None, left=False, right=True): cdef l_int N, info, M, MM cdef char *side cdef char *howmny cdef np.ndarray[np.complex128_t, ndim=2] vl, vr cdef double complex *vl_ptr cdef double complex *vr_ptr cdef l_logical *select_ptr cdef np.ndarray[np.complex128_t] work cdef np.ndarray[np.float64_t] rwork
side = "L" else: return
else: else: howmny = "A"
else: else:
else: else: vr_ptr = NULL
&N, <double complex *>T.data, &N, vl_ptr, &N, vr_ptr, &N, &MM, &M, <double complex *>work.data, <double *>rwork.data, &info)
return vl else:
# wrappers for xGGES np.ndarray[np.float32_t, ndim=2] B, calc_q=True, calc_z=True, calc_ev=True): cdef l_int N, lwork, sdim, info cdef char *jobvsl cdef char *jobvsr cdef float *vsl_ptr cdef float *vsr_ptr cdef float qwork cdef np.ndarray[np.float32_t, ndim=2] vsl, vsr cdef np.ndarray[np.float32_t] alphar, alphai, beta, work
else: vsl = None vsl_ptr = NULL jobvsl = "N"
else: vsr = None vsr_ptr = NULL jobvsr = "N"
# workspace query &N, <float *>A.data, &N, <float *>B.data, &N, &sdim, <float *>alphar.data, <float *>alphai.data, <float *>beta.data, vsl_ptr, &N, vsr_ptr, &N, &qwork, &lwork, NULL, &info)
# Now the real calculation &N, <float *>A.data, &N, <float *>B.data, &N, &sdim, <float *>alphar.data, <float *>alphai.data, <float *>beta.data, vsl_ptr, &N, vsr_ptr, &N, <float *>work.data, &lwork, NULL, &info)
raise LinAlgError("QZ iteration failed to converge in sgges")
else: alpha = alphar
np.ndarray[np.float64_t, ndim=2] B, calc_q=True, calc_z=True, calc_ev=True): cdef l_int N, lwork, sdim, info cdef char *jobvsl cdef char *jobvsr cdef double *vsl_ptr cdef double *vsr_ptr cdef double qwork cdef np.ndarray[np.float64_t, ndim=2] vsl, vsr cdef np.ndarray[np.float64_t] alphar, alphai, beta, work
else:
else: vsr = None vsr_ptr = NULL jobvsr = "N"
# workspace query &N, <double *>A.data, &N, <double *>B.data, &N, &sdim, <double *>alphar.data, <double *>alphai.data, <double *>beta.data, vsl_ptr, &N, vsr_ptr, &N, &qwork, &lwork, NULL, &info)
# Now the real calculation &N, <double *>A.data, &N, <double *>B.data, &N, &sdim, <double *>alphar.data, <double *>alphai.data, <double *>beta.data, vsl_ptr, &N, vsr_ptr, &N, <double *>work.data, &lwork, NULL, &info)
raise LinAlgError("QZ iteration failed to converge in dgges")
else:
np.ndarray[np.complex64_t, ndim=2] B, calc_q=True, calc_z=True, calc_ev=True): cdef l_int N, lwork, sdim, info cdef char *jobvsl cdef char *jobvsr cdef float complex *vsl_ptr cdef float complex *vsr_ptr cdef float complex qwork cdef np.ndarray[np.complex64_t, ndim=2] vsl, vsr cdef np.ndarray[np.complex64_t] alpha, beta, work cdef np.ndarray[np.float32_t] rwork
else: vsl = None vsl_ptr = NULL jobvsl = "N"
else: vsr = None vsr_ptr = NULL jobvsr = "N"
# workspace query &N, <float complex *>A.data, &N, <float complex *>B.data, &N, &sdim, <float complex *>alpha.data, <float complex *>beta.data, vsl_ptr, &N, vsr_ptr, &N, &qwork, &lwork, <float *>rwork.data, NULL, &info)
# Now the real calculation &N, <float complex *>A.data, &N, <float complex *>B.data, &N, &sdim, <float complex *>alpha.data, <float complex *>beta.data, vsl_ptr, &N, vsr_ptr, &N, <float complex *>work.data, &lwork, <float *>rwork.data, NULL, &info)
raise LinAlgError("QZ iteration failed to converge in cgges")
np.ndarray[np.complex128_t, ndim=2] B, calc_q=True, calc_z=True, calc_ev=True): cdef l_int N, lwork, sdim, info cdef char *jobvsl cdef char *jobvsr cdef double complex *vsl_ptr cdef double complex *vsr_ptr cdef double complex qwork cdef np.ndarray[np.complex128_t, ndim=2] vsl, vsr cdef np.ndarray[np.complex128_t] alpha, beta, work cdef np.ndarray[np.float64_t] rwork
else: vsl = None vsl_ptr = NULL jobvsl = "N"
else: vsr = None vsr_ptr = NULL jobvsr = "N"
# workspace query &N, <double complex *>A.data, &N, <double complex *>B.data, &N, &sdim, <double complex *>alpha.data, <double complex *>beta.data, vsl_ptr, &N, vsr_ptr, &N, &qwork, &lwork, <double *>rwork.data, NULL, &info)
# Now the real calculation &N, <double complex *>A.data, &N, <double complex *>B.data, &N, &sdim, <double complex *>alpha.data, <double complex *>beta.data, vsl_ptr, &N, vsr_ptr, &N, <double complex *>work.data, &lwork, <double *>rwork.data, NULL, &info)
raise LinAlgError("QZ iteration failed to converge in zgges")
# wrappers for xTGSEN np.ndarray[np.float32_t, ndim=2] S, np.ndarray[np.float32_t, ndim=2] T, np.ndarray[np.float32_t, ndim=2] Q=None, np.ndarray[np.float32_t, ndim=2] Z=None, calc_ev=True): cdef l_int N, M, lwork, liwork, qiwork, info, ijob cdef l_logical wantq, wantz cdef float qwork cdef float *q_ptr cdef float *z_ptr cdef np.ndarray[np.float32_t] alphar, alphai, beta, work cdef np.ndarray[l_int] iwork
assert_fortran_mat(S, T, Q, Z)
N = S.shape[0] alphar = np.empty(N, dtype = np.float32) alphai = np.empty(N, dtype = np.float32) beta = np.empty(N, dtype = np.float32) ijob = 0
if Q is not None: wantq = 1 q_ptr = <float *>Q.data else: wantq = 0 q_ptr = NULL
if Z is not None: wantz = 1 z_ptr = <float *>Z.data else: wantz = 0 z_ptr = NULL
# workspace query lwork = -1 liwork = -1 f_lapack.stgsen_(&ijob, &wantq, &wantz, <l_logical *>select.data, &N, <float *>S.data, &N, <float *>T.data, &N, <float *>alphar.data, <float *>alphai.data, <float *>beta.data, q_ptr, &N, z_ptr, &N, &M, NULL, NULL, NULL, &qwork, &lwork, &qiwork, &liwork, &info)
assert info == 0, "Argument error in stgsen"
lwork = <int>qwork work = np.empty(lwork, dtype = np.float32) liwork = qiwork iwork = np.empty(liwork, dtype = int_dtype)
# Now the real calculation f_lapack.stgsen_(&ijob, &wantq, &wantz, <l_logical *>select.data, &N, <float *>S.data, &N, <float *>T.data, &N, <float *>alphar.data, <float *>alphai.data, <float *>beta.data, q_ptr, &N, z_ptr, &N, &M, NULL, NULL, NULL, <float *>work.data, &lwork, <l_int *>iwork.data, &liwork, &info)
if info > 0: raise LinAlgError("Reordering failed; problem is very ill-conditioned")
assert info == 0, "Argument error in stgsen"
if alphai.nonzero()[0].size: alpha = alphar + 1j * alphai else: alpha = alphar
return filter_args((True, True, Q is not None, Z is not None, calc_ev, calc_ev), (S, T, Q, Z, alpha, beta))
np.ndarray[np.float64_t, ndim=2] S, np.ndarray[np.float64_t, ndim=2] T, np.ndarray[np.float64_t, ndim=2] Q=None, np.ndarray[np.float64_t, ndim=2] Z=None, calc_ev=True): cdef l_int N, M, lwork, liwork, qiwork, info, ijob cdef l_logical wantq, wantz cdef double qwork cdef double *q_ptr cdef double *z_ptr cdef np.ndarray[np.float64_t] alphar, alphai, beta, work cdef np.ndarray[l_int] iwork
wantq = 1 q_ptr = <double *>Q.data else:
else: wantz = 0 z_ptr = NULL
# workspace query &N, <double *>S.data, &N, <double *>T.data, &N, <double *>alphar.data, <double *>alphai.data, <double *>beta.data, q_ptr, &N, z_ptr, &N, &M, NULL, NULL, NULL, &qwork, &lwork, &qiwork, &liwork, &info)
# Now the real calculation &N, <double *>S.data, &N, <double *>T.data, &N, <double *>alphar.data, <double *>alphai.data, <double *>beta.data, q_ptr, &N, z_ptr, &N, &M, NULL, NULL, NULL, <double *>work.data, &lwork, <l_int *>iwork.data, &liwork, &info)
raise LinAlgError("Reordering failed; problem is very ill-conditioned")
else:
np.ndarray[np.complex64_t, ndim=2] S, np.ndarray[np.complex64_t, ndim=2] T, np.ndarray[np.complex64_t, ndim=2] Q=None, np.ndarray[np.complex64_t, ndim=2] Z=None, calc_ev=True): cdef l_int N, M, lwork, liwork, qiwork, info, ijob cdef l_logical wantq, wantz cdef float complex qwork cdef float complex *q_ptr cdef float complex *z_ptr cdef np.ndarray[np.complex64_t] alpha, beta, work cdef np.ndarray[l_int] iwork
else: wantq = 0 q_ptr = NULL
else: wantz = 0 z_ptr = NULL
# workspace query &N, <float complex *>S.data, &N, <float complex *>T.data, &N, <float complex *>alpha.data, <float complex *>beta.data, q_ptr, &N, z_ptr, &N, &M, NULL, NULL, NULL, &qwork, &lwork, &qiwork, &liwork, &info)
# Now the real calculation &N, <float complex *>S.data, &N, <float complex *>T.data, &N, <float complex *>alpha.data, <float complex *>beta.data, q_ptr, &N, z_ptr, &N, &M, NULL, NULL, NULL, <float complex *>work.data, &lwork, <l_int *>iwork.data, &liwork, &info)
raise LinAlgError("Reordering failed; problem is very ill-conditioned")
np.ndarray[np.complex128_t, ndim=2] S, np.ndarray[np.complex128_t, ndim=2] T, np.ndarray[np.complex128_t, ndim=2] Q=None, np.ndarray[np.complex128_t, ndim=2] Z=None, calc_ev=True): cdef l_int N, M, lwork, liwork, qiwork, info, ijob cdef l_logical wantq, wantz cdef double complex qwork cdef double complex *q_ptr cdef double complex *z_ptr cdef np.ndarray[np.complex128_t] alpha, beta, work cdef np.ndarray[l_int] iwork
else: wantq = 0 q_ptr = NULL
else: wantz = 0 z_ptr = NULL
# workspace query &N, <double complex *>S.data, &N, <double complex *>T.data, &N, <double complex *>alpha.data, <double complex *>beta.data, q_ptr, &N, z_ptr, &N, &M, NULL, NULL, NULL, &qwork, &lwork, &qiwork, &liwork, &info)
# Now the real calculation &N, <double complex *>S.data, &N, <double complex *>T.data, &N, <double complex *>alpha.data, <double complex *>beta.data, q_ptr, &N, z_ptr, &N, &M, NULL, NULL, NULL, <double complex *>work.data, &lwork, <l_int *>iwork.data, &liwork, &info)
raise LinAlgError("Reordering failed; problem is very ill-conditioned")
# xTGEVC np.ndarray[np.float32_t, ndim=2] T, np.ndarray[np.float32_t, ndim=2] Q=None, np.ndarray[np.float32_t, ndim=2] Z=None, np.ndarray[l_logical] select=None, left=False, right=True): cdef l_int N, info, M, MM cdef char *side cdef char *howmny cdef np.ndarray[np.float32_t, ndim=2] vl_r, vr_r cdef float *vl_r_ptr cdef float *vr_r_ptr cdef np.ndarray[l_logical] select_cpy cdef l_logical *select_ptr cdef np.ndarray[np.float32_t] work
elif left: side = "L" elif right: side = "R" else: return
# Correct for possible additional storage if a single complex # eigenvalue is selected. # For that: Figure out the positions of the 2x2 blocks. MM += 1
# select is overwritten in stgevc order = 'F') else: (left and not right and Q is not None) or (right and not left and Z is not None)): else: howmny = "A"
else: else: vl_r_ptr = NULL
else: else: vr_r_ptr = NULL
&N, <float *>S.data, &N, <float *>T.data, &N, vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M, <float *>work.data, &info)
# If there are complex eigenvalues, we need to postprocess the eigenvectors else: if left: vl = vl_r if right: vr = vr_r
elif left: return vl else: return vr
np.ndarray[np.float64_t, ndim=2] T, np.ndarray[np.float64_t, ndim=2] Q=None, np.ndarray[np.float64_t, ndim=2] Z=None, np.ndarray[l_logical] select=None, left=False, right=True): cdef l_int N, info, M, MM cdef char *side cdef char *howmny cdef np.ndarray[np.float64_t, ndim=2] vl_r, vr_r cdef double *vl_r_ptr cdef double *vr_r_ptr cdef np.ndarray[l_logical] select_cpy cdef l_logical *select_ptr cdef np.ndarray[np.float64_t] work
side = "L" else: return
# Correct for possible additional storage if a single complex # eigenvalue is selected. # For that: Figure out the positions of the 2x2 blocks. MM += 1
# select is overwritten in dtgevc order = 'F') else: (left and not right and Q is not None) or (right and not left and Z is not None)): else: howmny = "A"
else: else:
else: else: vr_r_ptr = NULL
&N, <double *>S.data, &N, <double *>T.data, &N, vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M, <double *>work.data, &info)
# If there are complex eigenvalues, we need to postprocess the # eigenvectors. else:
return vl else:
np.ndarray[np.complex64_t, ndim=2] T, np.ndarray[np.complex64_t, ndim=2] Q=None, np.ndarray[np.complex64_t, ndim=2] Z=None, np.ndarray[l_logical] select=None, left=False, right=True): cdef l_int N, info, M, MM cdef char *side cdef char *howmny cdef np.ndarray[np.complex64_t, ndim=2] vl, vr cdef float complex *vl_ptr cdef float complex *vr_ptr cdef l_logical *select_ptr cdef np.ndarray[np.complex64_t] work cdef np.ndarray[np.float32_t] rwork
elif left: side = "L" elif right: side = "R" else: return
else: (left and not right and Q is not None) or (right and not left and Z is not None)): else: howmny = "A"
else: else: vl_ptr = NULL
else: else: vr_ptr = NULL
&N, <float complex *>S.data, &N, <float complex *>T.data, &N, vl_ptr, &N, vr_ptr, &N, &MM, &M, <float complex *>work.data, <float *>rwork.data, &info)
elif left: return vl else: return vr
np.ndarray[np.complex128_t, ndim=2] T, np.ndarray[np.complex128_t, ndim=2] Q=None, np.ndarray[np.complex128_t, ndim=2] Z=None, np.ndarray[l_logical] select=None, left=False, right=True): cdef l_int N, info, M, MM cdef char *side cdef char *howmny cdef np.ndarray[np.complex128_t, ndim=2] vl, vr cdef double complex *vl_ptr cdef double complex *vr_ptr cdef l_logical *select_ptr cdef np.ndarray[np.complex128_t] work cdef np.ndarray[np.float64_t] rwork
elif left: side = "L" elif right: side = "R" else: return
else: (left and not right and Q is not None) or (right and not left and Z is not None)): else: howmny = "A"
else: else: vl_ptr = NULL
else: else: vr_ptr = NULL
&N, <double complex *>S.data, &N, <double complex *>T.data, &N, vl_ptr, &N, vr_ptr, &N, &MM, &M, <double complex *>work.data, <double *>rwork.data, &info)
elif left: return vl else: return vr
"""Convert arrays to Fortran format.
This function takes a number of array objects in `args` and converts them to a format that can be directly passed to a Fortran function (Fortran contiguous NumPy array). If the arrays have different data type, they converted arrays are cast to a common compatible data type (one of NumPy's `float32`, `float64`, `complex64`, `complex128` data types).
If `overwrite` is ``False``, an NumPy array that would already be in the correct format (Fortran contiguous, right data type) is neverthelessed copied. (Hence, overwrite = True does not imply that acting on the converted array in the return values will overwrite the original array in all cases -- it does only so if the original array was already in the correct format. The conversions require copying. In fact, that's the same behavior as in SciPy, it's just not explicitly stated there)
If an argument is ``None``, it is just passed through and not used to determine the proper LAPACK type.
`prepare_for_lapack` returns a character indicating the proper LAPACK data type ('s', 'd', 'c', 'z') and a list of properly converted arrays. """
# Make sure we have NumPy arrays raise ValueError("Argument cannot be interpreted " "as a numeric array")
else:
# First figure out common dtype # Note: The return type of common_type is guaranteed to be a floating point # kind.
else: raise AssertionError("Unexpected data type from common_type")
# Now make sure that the array is contiguous, and copy if necessary. # ugly here: copy makes always C-array, no way to tell it # to make a Fortran array. npmat = np.ascontiguousarray(npmat, dtype = dtype) else: raise ValueError("Dimensionality of array is not 1 or 2")
|